diff --git a/include/AdePT/core/AsyncAdePTTransport.cuh b/include/AdePT/core/AsyncAdePTTransport.cuh
index 77fbd854..f904c2b4 100644
--- a/include/AdePT/core/AsyncAdePTTransport.cuh
+++ b/include/AdePT/core/AsyncAdePTTransport.cuh
@@ -172,37 +172,37 @@ __device__ inline uint64_t GenerateSeedFromTrackInfo(const AsyncAdePT::TrackData
 }
 
 // Kernel function to initialize tracks comming from a Geant4 buffer
-__global__ void InitTracks(AsyncAdePT::TrackDataWithIDs *trackinfo, int ntracks, Secondaries secondaries,
+__global__ void InitTracks(AsyncAdePT::TrackDataWithIDs *trackinfo, int ntracks, ParticleManager particleManager,
                            const vecgeom::VPlacedVolume *world, adept::MParrayT<QueueIndexPair> *toBeEnqueued,
                            uint64_t initialSeed)
 {
   // constexpr double tolerance = 10. * vecgeom::kTolerance;
   for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < ntracks; i += blockDim.x * gridDim.x) {
-    ParticleGenerator *generator = nullptr;
-    const auto &trackInfo        = trackinfo[i];
-    short queueIndex             = -1;
+    SpeciesParticleManager *speciesTM = nullptr;
+    const auto &trackInfo             = trackinfo[i];
+    short queueIndex                  = -1;
     switch (trackInfo.pdg) {
     case 11:
-      generator  = &secondaries.electrons;
+      speciesTM  = &particleManager.electrons;
       queueIndex = ParticleType::Electron;
       break;
     case -11:
-      generator  = &secondaries.positrons;
+      speciesTM  = &particleManager.positrons;
       queueIndex = ParticleType::Positron;
       break;
     case 22:
-      generator  = &secondaries.gammas;
+      speciesTM  = &particleManager.gammas;
       queueIndex = ParticleType::Gamma;
     };
-    assert(generator != nullptr && "Unsupported pdg type");
+    assert(speciesTM != nullptr && "Unsupported pdg type");
 
     // TODO: Delay when not enough slots?
-    const auto slot = generator->NextSlot();
+    const auto slot = speciesTM->NextSlot();
     // we need to scramble the initial seed with some more trackinfo to generate a unique seed.
     // otherwise, if a particle returns from the device and is injected again (i.e., via lepton nuclear), it would have
     // the same random number state, causing collisions in the track IDs
     auto seed    = GenerateSeedFromTrackInfo(trackInfo, initialSeed);
-    Track &track = generator->InitTrack(
+    Track &track = speciesTM->InitTrack(
         slot, seed, trackInfo.eKin, trackInfo.globalTime, static_cast<float>(trackInfo.localTime),
         static_cast<float>(trackInfo.properTime), trackInfo.weight, trackInfo.position, trackInfo.direction,
         trackInfo.eventId, trackInfo.trackId, trackInfo.parentId, trackInfo.threadId, trackInfo.stepCounter);
@@ -952,13 +952,15 @@ void TransportLoop(int trackCapacity, int leakCapacity, int scoringCapacity, int
       positrons.queues.SwapActive();
       gammas.queues.SwapActive();
 
-      const Secondaries secondaries = {
-          .electrons = {electrons.tracks, electrons.slotManager, electrons.slotManagerLeaks,
-                        electrons.queues.nextActive},
-          .positrons = {positrons.tracks, positrons.slotManager, positrons.slotManagerLeaks,
-                        positrons.queues.nextActive},
-          .gammas    = {gammas.tracks, gammas.slotManager, gammas.slotManagerLeaks, gammas.queues.nextActive},
-      };
+      const ParticleManager particleManager = {
+          .electrons = {electrons.tracks, electrons.leaks, electrons.slotManager, electrons.slotManagerLeaks,
+                        electrons.queues.initiallyActive, electrons.queues.nextActive,
+                        electrons.queues.leakedTracksCurrent},
+          .positrons = {positrons.tracks, positrons.leaks, positrons.slotManager, positrons.slotManagerLeaks,
+                        positrons.queues.initiallyActive, positrons.queues.nextActive,
+                        positrons.queues.leakedTracksCurrent},
+          .gammas    = {gammas.tracks, gammas.leaks, gammas.slotManager, gammas.slotManagerLeaks,
+                        gammas.queues.initiallyActive, gammas.queues.nextActive, gammas.queues.leakedTracksCurrent}};
       const AllParticleQueues allParticleQueues = {{electrons.queues, positrons.queues, gammas.queues}};
 #ifdef USE_SPLIT_KERNELS
       const AllInteractionQueues allGammaInteractionQueues = {
@@ -1009,7 +1011,7 @@ void TransportLoop(int trackCapacity, int leakCapacity, int scoringCapacity, int
           constexpr auto injectThreads = 128u;
           const auto injectBlocks      = (nInject + injectThreads - 1) / injectThreads;
           InitTracks<<<injectBlocks, injectThreads, 0, injectStream>>>(
-              trackBuffer.toDevice_dev.get(), nInject, secondaries, world_dev, gpuState.injectionQueue, adeptSeed);
+              trackBuffer.toDevice_dev.get(), nInject, particleManager, world_dev, gpuState.injectionQueue, adeptSeed);
           COPCORE_CUDA_CHECK(cudaLaunchHostFunc(
               injectStream,
               [](void *arg) { (*static_cast<decltype(GPUstate::injectState) *>(arg)) = InjectState::ReadyToEnqueue; },
@@ -1057,33 +1059,28 @@ void TransportLoop(int trackCapacity, int leakCapacity, int scoringCapacity, int
         const auto [threads, blocks] = computeThreadsAndBlocks(particlesInFlight[ParticleType::Electron]);
 #ifdef USE_SPLIT_KERNELS
         ElectronHowFar<true><<<blocks, threads, 0, electrons.stream>>>(
-            electrons.tracks, gpuState.hepEmBuffers_d.electronsHepEm, electrons.queues.initiallyActive,
-            electrons.queues.nextActive, electrons.queues.propagation, electrons.queues.leakedTracksCurrent,
-            gpuState.stats_dev, allowFinishOffEvent);
-        ElectronPropagation<true><<<blocks, threads, 0, electrons.stream>>>(
-            electrons.tracks, gpuState.hepEmBuffers_d.electronsHepEm, electrons.queues.propagation,
-            electrons.queues.leakedTracksCurrent);
+            particleManager, gpuState.hepEmBuffers_d.electronsHepEm, electrons.queues.propagation, gpuState.stats_dev,
+            allowFinishOffEvent);
+        ElectronPropagation<true>
+            <<<blocks, threads, 0, electrons.stream>>>(particleManager, gpuState.hepEmBuffers_d.electronsHepEm);
         ElectronMSC<true><<<blocks, threads, 0, electrons.stream>>>(
             electrons.tracks, gpuState.hepEmBuffers_d.electronsHepEm, electrons.queues.propagation);
         ElectronSetupInteractions<true, PerEventScoring><<<blocks, threads, 0, electrons.stream>>>(
-            electrons.tracks, electrons.leaks, gpuState.hepEmBuffers_d.electronsHepEm, electrons.queues.propagation,
-            secondaries, electrons.queues.nextActive, allElectronInteractionQueues,
-            electrons.queues.leakedTracksCurrent, gpuState.fScoring_dev, returnAllSteps, returnLastStep);
+            gpuState.hepEmBuffers_d.electronsHepEm, electrons.queues.propagation, particleManager,
+            allElectronInteractionQueues, gpuState.fScoring_dev, returnAllSteps, returnLastStep);
         ElectronRelocation<true, PerEventScoring><<<blocks, threads, 0, electrons.stream>>>(
-            electrons.tracks, electrons.leaks, gpuState.hepEmBuffers_d.electronsHepEm, secondaries,
-            electrons.queues.nextActive, electrons.queues.interactionQueues[4], electrons.queues.leakedTracksCurrent,
+            gpuState.hepEmBuffers_d.electronsHepEm, particleManager, electrons.queues.interactionQueues[4],
             gpuState.fScoring_dev, returnAllSteps, returnLastStep);
         ElectronIonization<true, PerEventScoring><<<blocks, threads, 0, electrons.stream>>>(
-            electrons.tracks, gpuState.hepEmBuffers_d.electronsHepEm, secondaries, electrons.queues.nextActive,
-            electrons.queues.interactionQueues[0], gpuState.fScoring_dev, returnAllSteps, returnLastStep);
+            gpuState.hepEmBuffers_d.electronsHepEm, particleManager, electrons.queues.interactionQueues[0],
+            gpuState.fScoring_dev, returnAllSteps, returnLastStep);
         ElectronBremsstrahlung<true, PerEventScoring><<<blocks, threads, 0, electrons.stream>>>(
-            electrons.tracks, gpuState.hepEmBuffers_d.electronsHepEm, secondaries, electrons.queues.nextActive,
-            electrons.queues.interactionQueues[1], gpuState.fScoring_dev, returnAllSteps, returnLastStep);
+            gpuState.hepEmBuffers_d.electronsHepEm, particleManager, electrons.queues.interactionQueues[1],
+            gpuState.fScoring_dev, returnAllSteps, returnLastStep);
 #else
         TransportElectrons<PerEventScoring, SteppingAction><<<blocks, threads, 0, electrons.stream>>>(
-            electrons.tracks, electrons.leaks, electrons.queues.initiallyActive, secondaries,
-            electrons.queues.nextActive, electrons.queues.leakedTracksCurrent, gpuState.fScoring_dev,
-            gpuState.stats_dev, steppingActionParams, allowFinishOffEvent, returnAllSteps, returnLastStep);
+            particleManager, gpuState.fScoring_dev, gpuState.stats_dev, steppingActionParams, allowFinishOffEvent,
+            returnAllSteps, returnLastStep);
 #endif
         COPCORE_CUDA_CHECK(cudaEventRecord(electrons.event, electrons.stream));
         COPCORE_CUDA_CHECK(cudaStreamWaitEvent(gpuState.stream, electrons.event, 0));
@@ -1098,39 +1095,34 @@ void TransportLoop(int trackCapacity, int leakCapacity, int scoringCapacity, int
         const auto [threads, blocks] = computeThreadsAndBlocks(particlesInFlight[ParticleType::Positron]);
 #ifdef USE_SPLIT_KERNELS
         ElectronHowFar<false><<<blocks, threads, 0, positrons.stream>>>(
-            positrons.tracks, gpuState.hepEmBuffers_d.positronsHepEm, positrons.queues.initiallyActive,
-            positrons.queues.nextActive, positrons.queues.propagation, positrons.queues.leakedTracksCurrent,
-            gpuState.stats_dev, allowFinishOffEvent);
-        ElectronPropagation<false><<<blocks, threads, 0, positrons.stream>>>(
-            positrons.tracks, gpuState.hepEmBuffers_d.positronsHepEm, positrons.queues.propagation,
-            positrons.queues.leakedTracksCurrent);
+            particleManager, gpuState.hepEmBuffers_d.positronsHepEm, positrons.queues.propagation, gpuState.stats_dev,
+            allowFinishOffEvent);
+        ElectronPropagation<false>
+            <<<blocks, threads, 0, positrons.stream>>>(particleManager, gpuState.hepEmBuffers_d.positronsHepEm);
         ElectronMSC<false><<<blocks, threads, 0, positrons.stream>>>(
             positrons.tracks, gpuState.hepEmBuffers_d.positronsHepEm, positrons.queues.propagation);
         ElectronSetupInteractions<false, PerEventScoring><<<blocks, threads, 0, positrons.stream>>>(
-            positrons.tracks, positrons.leaks, gpuState.hepEmBuffers_d.positronsHepEm, positrons.queues.propagation,
-            secondaries, positrons.queues.nextActive, allPositronInteractionQueues,
-            positrons.queues.leakedTracksCurrent, gpuState.fScoring_dev, returnAllSteps, returnLastStep);
+            gpuState.hepEmBuffers_d.positronsHepEm, positrons.queues.propagation, particleManager,
+            allPositronInteractionQueues, gpuState.fScoring_dev, returnAllSteps, returnLastStep);
         ElectronRelocation<false, PerEventScoring><<<blocks, threads, 0, positrons.stream>>>(
-            positrons.tracks, positrons.leaks, gpuState.hepEmBuffers_d.positronsHepEm, secondaries,
-            positrons.queues.nextActive, positrons.queues.interactionQueues[4], positrons.queues.leakedTracksCurrent,
+            gpuState.hepEmBuffers_d.positronsHepEm, particleManager, positrons.queues.interactionQueues[4],
             gpuState.fScoring_dev, returnAllSteps, returnLastStep);
         ElectronIonization<false, PerEventScoring><<<blocks, threads, 0, positrons.stream>>>(
-            positrons.tracks, gpuState.hepEmBuffers_d.positronsHepEm, secondaries, positrons.queues.nextActive,
-            positrons.queues.interactionQueues[0], gpuState.fScoring_dev, returnAllSteps, returnLastStep);
+            gpuState.hepEmBuffers_d.positronsHepEm, particleManager, positrons.queues.interactionQueues[0],
+            gpuState.fScoring_dev, returnAllSteps, returnLastStep);
         ElectronBremsstrahlung<false, PerEventScoring><<<blocks, threads, 0, positrons.stream>>>(
-            positrons.tracks, gpuState.hepEmBuffers_d.positronsHepEm, secondaries, positrons.queues.nextActive,
-            positrons.queues.interactionQueues[1], gpuState.fScoring_dev, returnAllSteps, returnLastStep);
+            gpuState.hepEmBuffers_d.positronsHepEm, particleManager, positrons.queues.interactionQueues[1],
+            gpuState.fScoring_dev, returnAllSteps, returnLastStep);
         PositronAnnihilation<PerEventScoring><<<blocks, threads, 0, positrons.stream>>>(
-            positrons.tracks, gpuState.hepEmBuffers_d.positronsHepEm, secondaries, positrons.queues.nextActive,
-            positrons.queues.interactionQueues[2], gpuState.fScoring_dev, returnAllSteps, returnLastStep);
+            gpuState.hepEmBuffers_d.positronsHepEm, particleManager, positrons.queues.interactionQueues[2],
+            gpuState.fScoring_dev, returnAllSteps, returnLastStep);
         PositronStoppedAnnihilation<PerEventScoring><<<blocks, threads, 0, positrons.stream>>>(
-            positrons.tracks, gpuState.hepEmBuffers_d.positronsHepEm, secondaries, positrons.queues.nextActive,
-            positrons.queues.interactionQueues[3], gpuState.fScoring_dev, returnAllSteps, returnLastStep);
+            gpuState.hepEmBuffers_d.positronsHepEm, particleManager, positrons.queues.interactionQueues[3],
+            gpuState.fScoring_dev, returnAllSteps, returnLastStep);
 #else
         TransportPositrons<PerEventScoring, SteppingAction><<<blocks, threads, 0, positrons.stream>>>(
-            positrons.tracks, positrons.leaks, positrons.queues.initiallyActive, secondaries,
-            positrons.queues.nextActive, positrons.queues.leakedTracksCurrent, gpuState.fScoring_dev,
-            gpuState.stats_dev, steppingActionParams, allowFinishOffEvent, returnAllSteps, returnLastStep);
+            particleManager, gpuState.fScoring_dev, gpuState.stats_dev, steppingActionParams, allowFinishOffEvent,
+            returnAllSteps, returnLastStep);
 #endif
 
         COPCORE_CUDA_CHECK(cudaEventRecord(positrons.event, positrons.stream));
@@ -1145,43 +1137,39 @@ void TransportLoop(int trackCapacity, int leakCapacity, int scoringCapacity, int
 
         const auto [threads, blocks] = computeThreadsAndBlocks(particlesInFlight[ParticleType::Gamma]);
 #ifdef USE_SPLIT_KERNELS
-        GammaHowFar<<<blocks, threads, 0, gammas.stream>>>(
-            gammas.tracks, gammas.leaks, gpuState.hepEmBuffers_d.gammasHepEm, gammas.queues.initiallyActive,
-            secondaries, gammas.queues.propagation, gammas.queues.leakedTracksCurrent, gpuState.stats_dev,
-            allowFinishOffEvent);
+        GammaHowFar<<<blocks, threads, 0, gammas.stream>>>(gpuState.hepEmBuffers_d.gammasHepEm, particleManager,
+                                                           gammas.queues.propagation, gpuState.stats_dev,
+                                                           allowFinishOffEvent);
         GammaPropagation<<<blocks, threads, 0, gammas.stream>>>(gammas.tracks, gpuState.hepEmBuffers_d.gammasHepEm,
                                                                 gammas.queues.propagation);
-        GammaSetupInteractions<PerEventScoring><<<blocks, threads, 0, gammas.stream>>>(
-            gammas.tracks, gammas.leaks, gpuState.hepEmBuffers_d.gammasHepEm, gammas.queues.propagation, secondaries,
-            gammas.queues.nextActive, allGammaInteractionQueues, gammas.queues.leakedTracksCurrent,
-            gpuState.fScoring_dev);
+        GammaSetupInteractions<PerEventScoring>
+            <<<blocks, threads, 0, gammas.stream>>>(gpuState.hepEmBuffers_d.gammasHepEm, gammas.queues.propagation,
+                                                    particleManager, allGammaInteractionQueues, gpuState.fScoring_dev);
         GammaRelocation<PerEventScoring><<<blocks, threads, 0, gammas.stream>>>(
-            gammas.tracks, gammas.leaks, gpuState.hepEmBuffers_d.gammasHepEm, secondaries, gammas.queues.nextActive,
-            gammas.queues.interactionQueues[4], gammas.queues.leakedTracksCurrent, gpuState.fScoring_dev,
-            returnAllSteps, returnLastStep);
+            gpuState.hepEmBuffers_d.gammasHepEm, particleManager, gammas.queues.interactionQueues[4],
+            gpuState.fScoring_dev, returnAllSteps, returnLastStep);
         // Copying the number of interacting tracks back to host and using this information to adjust
         // the launch parameters of the interactions kernel is complicated and expensive due to a
         // required additional kernel launch and copy. Instead, launch the kernels with the same
         // parameters, the unneeded threads inmediately return and become free again.
         GammaConversion<PerEventScoring><<<blocks, threads, 0, gammas.stream>>>(
-            gammas.tracks, gpuState.hepEmBuffers_d.gammasHepEm, secondaries, gammas.queues.nextActive,
-            gammas.queues.interactionQueues[0], gpuState.fScoring_dev, returnAllSteps, returnLastStep);
+            gpuState.hepEmBuffers_d.gammasHepEm, particleManager, gammas.queues.interactionQueues[0],
+            gpuState.fScoring_dev, returnAllSteps, returnLastStep);
         GammaCompton<PerEventScoring><<<blocks, threads, 0, gammas.stream>>>(
-            gammas.tracks, gpuState.hepEmBuffers_d.gammasHepEm, secondaries, gammas.queues.nextActive,
-            gammas.queues.interactionQueues[1], gpuState.fScoring_dev, returnAllSteps, returnLastStep);
+            gpuState.hepEmBuffers_d.gammasHepEm, particleManager, gammas.queues.interactionQueues[1],
+            gpuState.fScoring_dev, returnAllSteps, returnLastStep);
         GammaPhotoelectric<PerEventScoring><<<blocks, threads, 0, gammas.stream>>>(
-            gammas.tracks, gpuState.hepEmBuffers_d.gammasHepEm, secondaries, gammas.queues.nextActive,
-            gammas.queues.interactionQueues[2], gpuState.fScoring_dev, returnAllSteps, returnLastStep);
+            gpuState.hepEmBuffers_d.gammasHepEm, particleManager, gammas.queues.interactionQueues[2],
+            gpuState.fScoring_dev, returnAllSteps, returnLastStep);
 #else
         TransportGammas<PerEventScoring, SteppingAction><<<blocks, threads, 0, gammas.stream>>>(
-            gammas.tracks, gammas.leaks, gammas.queues.initiallyActive, secondaries, gammas.queues.nextActive,
-            gammas.queues.leakedTracksCurrent, gpuState.fScoring_dev, gpuState.stats_dev, steppingActionParams,
-            allowFinishOffEvent, returnAllSteps, returnLastStep); //, gpuState.gammaInteractions);
+            particleManager, gpuState.fScoring_dev, gpuState.stats_dev, steppingActionParams, allowFinishOffEvent,
+            returnAllSteps, returnLastStep); //, gpuState.gammaInteractions);
 #endif
 
         // constexpr unsigned int intThreads = 128;
         // ApplyGammaInteractions<PerEventScoring><<<dim3(20, 3, 1), intThreads, 0, gammas.stream>>>(
-        //     gammas.tracks, secondaries, gammas.queues.nextActive, gpuState.fScoring_dev,
+        //     gammas.tracks, particleManager, gammas.queues.nextActive, gpuState.fScoring_dev,
         //     gpuState.gammaInteractions);
 
         COPCORE_CUDA_CHECK(cudaEventRecord(gammas.event, gammas.stream));
diff --git a/include/AdePT/core/AsyncAdePTTransportStruct.cuh b/include/AdePT/core/AsyncAdePTTransportStruct.cuh
index f3711c43..9f3c3a54 100644
--- a/include/AdePT/core/AsyncAdePTTransportStruct.cuh
+++ b/include/AdePT/core/AsyncAdePTTransportStruct.cuh
@@ -26,24 +26,44 @@
 namespace AsyncAdePT {
 
 // A bundle of pointers to generate particles of an implicit type.
-struct ParticleGenerator {
+struct SpeciesParticleManager {
   Track *fTracks;
+  Track *fLeakedTracks;
   SlotManager *fSlotManager;
   SlotManager *fSlotManagerLeaks;
   adept::MParray *fActiveQueue;
+  adept::MParray *fNextActiveQueue;
+  adept::MParray *fActiveLeaksQueue;
 
 public:
-  __host__ __device__ ParticleGenerator(Track *tracks, SlotManager *slotManager, SlotManager *slotManagerLeaks,
-                                        adept::MParray *activeQueue)
-      : fTracks(tracks), fSlotManager(slotManager), fSlotManagerLeaks(slotManagerLeaks), fActiveQueue(activeQueue)
+  __host__ __device__ SpeciesParticleManager(Track *tracks, Track *leakedTracks, SlotManager *slotManager,
+                                             SlotManager *slotManagerLeaks, adept::MParray *activeQueue,
+                                             adept::MParray *nextActiveQueue, adept::MParray *activeLeaksQueue)
+      : fTracks(tracks), fLeakedTracks(leakedTracks), fSlotManager(slotManager), fSlotManagerLeaks(slotManagerLeaks),
+        fActiveQueue(activeQueue), fNextActiveQueue(nextActiveQueue), fActiveLeaksQueue(activeLeaksQueue)
   {
   }
 
+  /// Obtain track and leaked track at given slot position
+  __device__ __forceinline__ Track &TrackAt(SlotManager::value_type slot) { return fTracks[slot]; }
+  __device__ __forceinline__ Track &LeakTrackAt(SlotManager::value_type slot) { return fLeakedTracks[slot]; }
+
   /// Obtain a slot for a track, but don't enqueue.
   __device__ auto NextSlot() { return fSlotManager->NextSlot(); }
-
   __device__ auto NextLeakSlot() { return fSlotManagerLeaks->NextSlot(); }
 
+  // enqueue into next-active queue
+  __device__ __forceinline__ bool EnqueueNext(SlotManager::value_type slot)
+  {
+    return fNextActiveQueue->push_back(slot);
+  }
+
+  // size of the active queue
+  __device__ __forceinline__ int ActiveSize() const { return fActiveQueue->size(); }
+
+  // read slot from active queue by index
+  __device__ __forceinline__ SlotManager::value_type ActiveAt(int i) const { return (*fActiveQueue)[i]; }
+
   /// Construct a track at the given location, forwarding all arguments to the constructor.
   template <typename... Ts>
   __device__ Track &InitTrack(SlotManager::value_type slot, Ts &&...args)
@@ -56,12 +76,36 @@ public:
   __device__ Track &NextTrack(Ts &&...args)
   {
     const auto slot = NextSlot();
-    fActiveQueue->push_back(slot);
+    // next track is only visible in next GPU iteration, therefore pushed in the NextActiveQueue
+    fNextActiveQueue->push_back(slot);
     auto &track = InitTrack(slot, std::forward<Ts>(args)...);
     return track;
   }
 
-  void SetActiveQueue(adept::MParray *queue) { fActiveQueue = queue; }
+  /// Obtains a leak slot, copies the track from the source slot to the leaks, and marks the slot in the active queue
+  /// for freeing
+  __device__ __forceinline__ void CopyTrackToLeaked(SlotManager::value_type srcSlot)
+  {
+    // get a leak slot
+    const auto leakSlot = NextLeakSlot();
+
+    // Create and construct track from other track
+    new (fLeakedTracks + leakSlot) Track{TrackAt(srcSlot)};
+
+    // enqueue into leak queue
+    const bool success = fActiveLeaksQueue->push_back(leakSlot);
+    if (!success) {
+      printf("ERROR: No space left in leaks queue.\n"
+             "\tThe threshold for flushing the leak buffer may be too high\n"
+             "\tThe space allocated to the leak buffer may be too small\n");
+      asm("trap;");
+    }
+
+    // free the source slot
+    fSlotManager->MarkSlotForFreeing(srcSlot);
+
+    return;
+  }
 };
 
 struct LeakedTracks {
@@ -86,10 +130,10 @@ struct LeakedTracks {
 // };
 
 // A bundle of generators for the three particle types.
-struct Secondaries {
-  ParticleGenerator electrons;
-  ParticleGenerator positrons;
-  ParticleGenerator gammas;
+struct ParticleManager {
+  SpeciesParticleManager electrons;
+  SpeciesParticleManager positrons;
+  SpeciesParticleManager gammas;
 };
 
 // Holds the leaked track structs for all three particle types
diff --git a/include/AdePT/core/Track.cuh b/include/AdePT/core/Track.cuh
index 993dc79d..cc20590b 100644
--- a/include/AdePT/core/Track.cuh
+++ b/include/AdePT/core/Track.cuh
@@ -77,6 +77,9 @@ struct Track {
 
   LeakStatus leakStatus{LeakStatus::NoLeak};
 
+  __host__ __device__ Track(const Track &)            = default;
+  __host__ __device__ Track &operator=(const Track &) = default;
+
   /// Construct a new track for GPU transport.
   /// NB: The navState remains uninitialised.
   __device__ Track(uint64_t rngSeed, double eKin, double globalTime, float localTime, float properTime, float weight,
diff --git a/include/AdePT/kernels/electrons.cuh b/include/AdePT/kernels/electrons.cuh
index fd1a494f..5f6c23c8 100644
--- a/include/AdePT/kernels/electrons.cuh
+++ b/include/AdePT/kernels/electrons.cuh
@@ -48,9 +48,7 @@ namespace AsyncAdePT {
 // applying the continuous effects and maybe a discrete process that could
 // generate secondaries.
 template <bool IsElectron, typename Scoring, class SteppingActionT>
-static __device__ __forceinline__ void TransportElectrons(Track *electrons, Track *leaks, const adept::MParray *active,
-                                                          Secondaries &secondaries, adept::MParray *nextActiveQueue,
-                                                          adept::MParray *leakedQueue, Scoring *userScoring,
+static __device__ __forceinline__ void TransportElectrons(ParticleManager &particleManager, Scoring *userScoring,
                                                           Stats *InFlightStats, const StepActionParam params,
                                                           AllowFinishOffEventArray allowFinishOffEvent,
                                                           const bool returnAllSteps, const bool returnLastStep)
@@ -73,12 +71,13 @@ static __device__ __forceinline__ void TransportElectrons(Track *electrons, Trac
 
   auto &magneticField = *gMagneticField;
 
-  int activeSize = active->size();
-  for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < activeSize; i += blockDim.x * gridDim.x) {
-    const int slot           = (*active)[i];
-    SlotManager &slotManager = IsElectron ? *secondaries.electrons.fSlotManager : *secondaries.positrons.fSlotManager;
+  auto &electronsOrPositrons = (IsElectron ? particleManager.electrons : particleManager.positrons);
+  SlotManager &slotManager   = *electronsOrPositrons.fSlotManager;
 
-    Track &currentTrack = electrons[slot];
+  const int activeSize = electronsOrPositrons.ActiveSize();
+  for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < activeSize; i += blockDim.x * gridDim.x) {
+    const auto slot     = electronsOrPositrons.ActiveAt(i);
+    Track &currentTrack = electronsOrPositrons.TrackAt(slot);
 #else
 
 template <bool IsElectron, typename Scoring, class SteppingActionT>
@@ -164,29 +163,11 @@ static __device__ __forceinline__ void TransportElectrons(adept::TrackManager<Tr
       currentTrack.navState   = nextState;
       currentTrack.leakStatus = leakReason;
 #ifdef ASYNC_MODE
-      // NOTE: When adapting the split kernels for async mode this won't
-      // work if we want to re-use slots on the fly. Directly copying to
-      // a trackdata struct would be better
       if (leakReason != LeakStatus::NoLeak) {
-        // Get a slot in the leaks array
-        int leakSlot;
-        if (IsElectron)
-          leakSlot = secondaries.electrons.NextLeakSlot();
-        else
-          leakSlot = secondaries.positrons.NextLeakSlot();
-        // Copy the track to the leaks array and store the index in the leak queue
-        leaks[leakSlot] = electrons[slot];
-        auto success    = leakedQueue->push_back(leakSlot);
-        if (!success) {
-          printf("ERROR: No space left in e-/+ leaks queue.\n\
-\tThe threshold for flushing the leak buffer may be too high\n\
-\tThe space allocated to the leak buffer may be too small\n");
-          asm("trap;");
-        }
-        // Free the slot in the tracks slot manager
-        slotManager.MarkSlotForFreeing(slot);
+        // Copy track at slot to the leaked tracks
+        electronsOrPositrons.CopyTrackToLeaked(slot);
       } else {
-        nextActiveQueue->push_back(slot);
+        electronsOrPositrons.EnqueueNext(slot);
       }
 #else
       currentTrack.CopyTo(trackdata, Pdg);
@@ -589,7 +570,7 @@ static __device__ __forceinline__ void TransportElectrons(adept::TrackManager<Tr
 
           } else {
 #ifdef ASYNC_MODE
-            Track &secondary = secondaries.electrons.NextTrack(
+            Track &secondary = particleManager.electrons.NextTrack(
                 newRNG, deltaEkin, pos, vecgeom::Vector3D<double>{dirSecondary[0], dirSecondary[1], dirSecondary[2]},
                 navState, currentTrack, globalTime);
 #else
@@ -671,7 +652,7 @@ static __device__ __forceinline__ void TransportElectrons(adept::TrackManager<Tr
 #endif
           } else {
 #ifdef ASYNC_MODE
-            Track &gamma = secondaries.gammas.NextTrack(
+            Track &gamma = particleManager.gammas.NextTrack(
                 newRNG, deltaEkin, pos, vecgeom::Vector3D<double>{dirSecondary[0], dirSecondary[1], dirSecondary[2]},
                 navState, currentTrack, globalTime);
 #else
@@ -750,7 +731,7 @@ static __device__ __forceinline__ void TransportElectrons(adept::TrackManager<Tr
 
           } else {
 #ifdef ASYNC_MODE
-            Track &gamma1 = secondaries.gammas.NextTrack(
+            Track &gamma1 = particleManager.gammas.NextTrack(
                 newRNG, theGamma1Ekin, pos,
                 vecgeom::Vector3D<double>{theGamma1Dir[0], theGamma1Dir[1], theGamma1Dir[2]}, navState, currentTrack,
                 globalTime);
@@ -794,7 +775,7 @@ static __device__ __forceinline__ void TransportElectrons(adept::TrackManager<Tr
 
           } else {
 #ifdef ASYNC_MODE
-            Track &gamma2 = secondaries.gammas.NextTrack(
+            Track &gamma2 = particleManager.gammas.NextTrack(
                 currentTrack.rngState, theGamma2Ekin, pos,
                 vecgeom::Vector3D<double>{theGamma2Dir[0], theGamma2Dir[1], theGamma2Dir[2]}, navState, currentTrack,
                 globalTime);
@@ -872,13 +853,14 @@ static __device__ __forceinline__ void TransportElectrons(adept::TrackManager<Tr
           RanluxppDouble newRNG2(currentTrack.rngState.Branch());
 
 #ifdef ASYNC_MODE
-          Track &gamma1 = secondaries.gammas.NextTrack(newRNG2, double{copcore::units::kElectronMassC2}, pos,
-                                                       vecgeom::Vector3D<double>{sint * cosPhi, sint * sinPhi, cost},
-                                                       navState, currentTrack, globalTime);
+          Track &gamma1 = particleManager.gammas.NextTrack(
+              newRNG2, double{copcore::units::kElectronMassC2}, pos,
+              vecgeom::Vector3D<double>{sint * cosPhi, sint * sinPhi, cost}, navState, currentTrack, globalTime);
 
           // Reuse the RNG state of the dying track.
-          Track &gamma2 = secondaries.gammas.NextTrack(currentTrack.rngState, double{copcore::units::kElectronMassC2},
-                                                       pos, -gamma1.dir, navState, currentTrack, globalTime);
+          Track &gamma2 =
+              particleManager.gammas.NextTrack(currentTrack.rngState, double{copcore::units::kElectronMassC2}, pos,
+                                               -gamma1.dir, navState, currentTrack, globalTime);
 #else
           Track &gamma1 = secondaries.gammas->NextTrack();
           Track &gamma2 = secondaries.gammas->NextTrack();
@@ -1024,26 +1006,20 @@ static __device__ __forceinline__ void TransportElectrons(adept::TrackManager<Tr
 // Instantiate kernels for electrons and positrons.
 #ifdef ASYNC_MODE
 template <typename Scoring, class SteppingActionT>
-__global__ void TransportElectrons(Track *electrons, Track *leaks, const adept::MParray *active,
-                                   Secondaries secondaries, adept::MParray *nextActiveQueue,
-                                   adept::MParray *leakedQueue, Scoring *userScoring, Stats *InFlightStats,
+__global__ void TransportElectrons(ParticleManager particleManager, Scoring *userScoring, Stats *InFlightStats,
                                    const StepActionParam params, AllowFinishOffEventArray allowFinishOffEvent,
                                    const bool returnAllSteps, const bool returnLastStep)
 {
   TransportElectrons</*IsElectron*/ true, Scoring, SteppingActionT>(
-      electrons, leaks, active, secondaries, nextActiveQueue, leakedQueue, userScoring, InFlightStats, params,
-      allowFinishOffEvent, returnAllSteps, returnLastStep);
+      particleManager, userScoring, InFlightStats, params, allowFinishOffEvent, returnAllSteps, returnLastStep);
 }
 template <typename Scoring, class SteppingActionT>
-__global__ void TransportPositrons(Track *positrons, Track *leaks, const adept::MParray *active,
-                                   Secondaries secondaries, adept::MParray *nextActiveQueue,
-                                   adept::MParray *leakedQueue, Scoring *userScoring, Stats *InFlightStats,
+__global__ void TransportPositrons(ParticleManager particleManager, Scoring *userScoring, Stats *InFlightStats,
                                    const StepActionParam params, AllowFinishOffEventArray allowFinishOffEvent,
                                    const bool returnAllSteps, const bool returnLastStep)
 {
   TransportElectrons</*IsElectron*/ false, Scoring, SteppingActionT>(
-      positrons, leaks, active, secondaries, nextActiveQueue, leakedQueue, userScoring, InFlightStats, params,
-      allowFinishOffEvent, returnAllSteps, returnLastStep);
+      particleManager, userScoring, InFlightStats, params, allowFinishOffEvent, returnAllSteps, returnLastStep);
 }
 #else
 template <typename Scoring, class SteppingActionT>
diff --git a/include/AdePT/kernels/electrons_split.cuh b/include/AdePT/kernels/electrons_split.cuh
index 2a7e6cb9..a13956bd 100644
--- a/include/AdePT/kernels/electrons_split.cuh
+++ b/include/AdePT/kernels/electrons_split.cuh
@@ -41,9 +41,8 @@ __device__ double GetVelocity(double eKin)
 namespace AsyncAdePT {
 
 template <bool IsElectron>
-__global__ void ElectronHowFar(Track *electrons, G4HepEmElectronTrack *hepEMTracks, const adept::MParray *active,
-                               adept::MParray *nextActiveQueue, adept::MParray *propagationQueue,
-                               adept::MParray *leakedQueue, Stats *InFlightStats,
+__global__ void ElectronHowFar(ParticleManager particleManager, G4HepEmElectronTrack *hepEMTracks,
+                               adept::MParray *propagationQueue, Stats *InFlightStats,
                                AllowFinishOffEventArray allowFinishOffEvent)
 {
   constexpr unsigned short maxSteps        = 10'000;
@@ -63,14 +62,16 @@ __global__ void ElectronHowFar(Track *electrons, G4HepEmElectronTrack *hepEMTrac
 
   auto &magneticField = *gMagneticField;
 
-  int activeSize = active->size();
+  auto &electronsOrPositrons = (IsElectron ? particleManager.electrons : particleManager.positrons);
+
+  const int activeSize = electronsOrPositrons.ActiveSize();
   for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < activeSize; i += blockDim.x * gridDim.x) {
-    const int slot      = (*active)[i];
-    Track &currentTrack = electrons[slot];
+    const auto slot     = electronsOrPositrons.ActiveAt(i);
+    Track &currentTrack = electronsOrPositrons.TrackAt(slot);
     // the MCC vector is indexed by the logical volume id
     const int lvolID = currentTrack.navState.GetLogicalId();
 
-    VolAuxData const &auxData = AsyncAdePT::gVolAuxData[lvolID]; // FIXME unify VolAuxData
+    VolAuxData const &auxData = AsyncAdePT::gVolAuxData[lvolID];
 
     currentTrack.preStepEKin       = currentTrack.eKin;
     currentTrack.preStepGlobalTime = currentTrack.globalTime;
@@ -87,20 +88,12 @@ __global__ void ElectronHowFar(Track *electrons, G4HepEmElectronTrack *hepEMTrac
     }
 
     auto survive = [&](LeakStatus leakReason = LeakStatus::NoLeak) {
-      // NOTE: When adapting the split kernels for async mode this won't
-      // work if we want to re-use slots on the fly. Directly copying to
-      // a trackdata struct would be better
       currentTrack.leakStatus = leakReason;
       if (leakReason != LeakStatus::NoLeak) {
-        auto success = leakedQueue->push_back(slot);
-        if (!success) {
-          printf("ERROR: No space left in e-/+ leaks queue.\n\
-\tThe threshold for flushing the leak buffer may be too high\n\
-\tThe space allocated to the leak buffer may be too small\n");
-          asm("trap;");
-        }
+        // Copy track at slot to the leaked tracks
+        electronsOrPositrons.CopyTrackToLeaked(slot);
       } else {
-        nextActiveQueue->push_back(slot);
+        electronsOrPositrons.EnqueueNext(slot);
       }
     };
 
@@ -199,8 +192,7 @@ __global__ void ElectronHowFar(Track *electrons, G4HepEmElectronTrack *hepEMTrac
 }
 
 template <bool IsElectron>
-__global__ void ElectronPropagation(Track *electrons, G4HepEmElectronTrack *hepEMTracks, const adept::MParray *active,
-                                    adept::MParray *leakedQueue)
+__global__ void ElectronPropagation(ParticleManager particleManager, G4HepEmElectronTrack *hepEMTracks)
 {
   constexpr double kPushDistance           = 1000 * vecgeom::kTolerance;
   constexpr int Charge                     = IsElectron ? -1 : 1;
@@ -221,11 +213,12 @@ __global__ void ElectronPropagation(Track *electrons, G4HepEmElectronTrack *hepE
 
   auto &magneticField = *gMagneticField;
 
-  int activeSize = active->size();
-  for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < activeSize; i += blockDim.x * gridDim.x) {
-    const int slot = (*active)[i];
+  auto &electronsOrPositrons = (IsElectron ? particleManager.electrons : particleManager.positrons);
 
-    Track &currentTrack = electrons[slot];
+  const int activeSize = electronsOrPositrons.ActiveSize();
+  for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < activeSize; i += blockDim.x * gridDim.x) {
+    const auto slot     = electronsOrPositrons.ActiveAt(i);
+    Track &currentTrack = electronsOrPositrons.TrackAt(slot);
     // the MCC vector is indexed by the logical volume id
     const int lvolID = currentTrack.navState.GetLogicalId();
 
@@ -371,22 +364,21 @@ __global__ void ElectronMSC(Track *electrons, G4HepEmElectronTrack *hepEMTracks,
  * @brief Adds tracks to interaction and relocation queues depending on their state
  */
 template <bool IsElectron, typename Scoring>
-__global__ void ElectronSetupInteractions(Track *electrons, Track *leaks, G4HepEmElectronTrack *hepEMTracks,
-                                          const adept::MParray *active, Secondaries secondaries,
-                                          adept::MParray *nextActiveQueue, AllInteractionQueues interactionQueues,
-                                          adept::MParray *leakedQueue, Scoring *userScoring, const bool returnAllSteps,
-                                          const bool returnLastStep)
+__global__ void ElectronSetupInteractions(G4HepEmElectronTrack *hepEMTracks, const adept::MParray *propagationQueue,
+                                          ParticleManager particleManager, AllInteractionQueues interactionQueues,
+                                          Scoring *userScoring, const bool returnAllSteps, const bool returnLastStep)
 {
-  int activeSize = active->size();
-  for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < activeSize; i += blockDim.x * gridDim.x) {
-    const int slot           = (*active)[i];
-    SlotManager &slotManager = IsElectron ? *secondaries.electrons.fSlotManager : *secondaries.positrons.fSlotManager;
+  auto &electronsOrPositrons = (IsElectron ? particleManager.electrons : particleManager.positrons);
+  SlotManager &slotManager   = *electronsOrPositrons.fSlotManager;
 
-    Track &currentTrack = electrons[slot];
+  int activeSize = propagationQueue->size();
+  for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < activeSize; i += blockDim.x * gridDim.x) {
+    const int slot      = (*propagationQueue)[i];
+    Track &currentTrack = electronsOrPositrons.TrackAt(slot);
     // the MCC vector is indexed by the logical volume id
     const int lvolID = currentTrack.navState.GetLogicalId();
 
-    VolAuxData const &auxData = AsyncAdePT::gVolAuxData[lvolID]; // FIXME unify VolAuxData
+    VolAuxData const &auxData = AsyncAdePT::gVolAuxData[lvolID];
 
     auto survive = [&](LeakStatus leakReason = LeakStatus::NoLeak) {
       // NOTE: When adapting the split kernels for async mode this won't
@@ -394,25 +386,10 @@ __global__ void ElectronSetupInteractions(Track *electrons, Track *leaks, G4HepE
       // a trackdata struct would be better
       currentTrack.leakStatus = leakReason;
       if (leakReason != LeakStatus::NoLeak) {
-        // Get a slot in the leaks array
-        int leakSlot;
-        if (IsElectron)
-          leakSlot = secondaries.electrons.NextLeakSlot();
-        else
-          leakSlot = secondaries.positrons.NextLeakSlot();
-        // Copy the track to the leaks array and store the index in the leak queue
-        leaks[leakSlot] = electrons[slot];
-        auto success    = leakedQueue->push_back(leakSlot);
-        if (!success) {
-          printf("ERROR: No space left in e-/+ leaks queue.\n\
-\tThe threshold for flushing the leak buffer may be too high\n\
-\tThe space allocated to the leak buffer may be too small\n");
-          asm("trap;");
-        }
-        // Free the slot in the tracks slot manager
-        slotManager.MarkSlotForFreeing(slot);
+        // Copy track at slot to the leaked tracks
+        electronsOrPositrons.CopyTrackToLeaked(slot);
       } else {
-        nextActiveQueue->push_back(slot);
+        electronsOrPositrons.EnqueueNext(slot);
       }
     };
 
@@ -540,22 +517,22 @@ __global__ void ElectronSetupInteractions(Track *electrons, Track *leaks, G4HepE
 }
 
 template <bool IsElectron, typename Scoring>
-__global__ void ElectronRelocation(Track *electrons, Track *leaks, G4HepEmElectronTrack *hepEMTracks,
-                                   Secondaries secondaries, adept::MParray *nextActiveQueue,
-                                   adept::MParray *relocatingQueue, adept::MParray *leakedQueue, Scoring *userScoring,
-                                   const bool returnAllSteps, const bool returnLastStep)
+__global__ void ElectronRelocation(G4HepEmElectronTrack *hepEMTracks, ParticleManager particleManager,
+                                   adept::MParray *relocatingQueue, Scoring *userScoring, const bool returnAllSteps,
+                                   const bool returnLastStep)
 {
   constexpr double kPushDistance = 1000 * vecgeom::kTolerance;
-  int activeSize                 = relocatingQueue->size();
-  for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < activeSize; i += blockDim.x * gridDim.x) {
-    const int slot           = (*relocatingQueue)[i];
-    SlotManager &slotManager = IsElectron ? *secondaries.electrons.fSlotManager : *secondaries.positrons.fSlotManager;
+  auto &electronsOrPositrons     = (IsElectron ? particleManager.electrons : particleManager.positrons);
 
-    Track &currentTrack = electrons[slot];
+  SlotManager &slotManager = *electronsOrPositrons.fSlotManager;
+  int activeSize           = relocatingQueue->size();
+  for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < activeSize; i += blockDim.x * gridDim.x) {
+    const int slot      = (*relocatingQueue)[i];
+    Track &currentTrack = electronsOrPositrons.TrackAt(slot);
     // the MCC vector is indexed by the logical volume id
     const int lvolID = currentTrack.navState.GetLogicalId();
 
-    VolAuxData const &auxData = AsyncAdePT::gVolAuxData[lvolID]; // FIXME unify VolAuxData
+    VolAuxData const &auxData = AsyncAdePT::gVolAuxData[lvolID];
 
     auto survive = [&](LeakStatus leakReason = LeakStatus::NoLeak) {
       // NOTE: When adapting the split kernels for async mode this won't
@@ -563,25 +540,10 @@ __global__ void ElectronRelocation(Track *electrons, Track *leaks, G4HepEmElectr
       // a trackdata struct would be better
       currentTrack.leakStatus = leakReason;
       if (leakReason != LeakStatus::NoLeak) {
-        // Get a slot in the leaks array
-        int leakSlot;
-        if (IsElectron)
-          leakSlot = secondaries.electrons.NextLeakSlot();
-        else
-          leakSlot = secondaries.positrons.NextLeakSlot();
-        // Copy the track to the leaks array and store the index in the leak queue
-        leaks[leakSlot] = electrons[slot];
-        auto success    = leakedQueue->push_back(leakSlot);
-        if (!success) {
-          printf("ERROR: No space left in e-/+ leaks queue.\n\
-\tThe threshold for flushing the leak buffer may be too high\n\
-\tThe space allocated to the leak buffer may be too small\n");
-          asm("trap;");
-        }
-        // Free the slot in the tracks slot manager
-        slotManager.MarkSlotForFreeing(slot);
+        // Copy track at slot to the leaked tracks
+        electronsOrPositrons.CopyTrackToLeaked(slot);
       } else {
-        nextActiveQueue->push_back(slot);
+        electronsOrPositrons.EnqueueNext(slot);
       }
     };
 
@@ -678,7 +640,7 @@ __global__ void ElectronRelocation(Track *electrons, Track *leaks, G4HepEmElectr
 
 template <bool IsElectron, typename Scoring>
 __device__ __forceinline__ void PerformStoppedAnnihilation(const int slot, Track &currentTrack,
-                                                           Secondaries &secondaries, double &energyDeposit,
+                                                           ParticleManager &particleManager, double &energyDeposit,
                                                            SlotManager &slotManager, const bool ApplyCuts,
                                                            const double theGammaCut, Scoring *userScoring,
                                                            const bool returnLastStep = false)
@@ -707,14 +669,15 @@ __device__ __forceinline__ void PerformStoppedAnnihilation(const int slot, Track
       currentTrack.rngState.Advance();
       RanluxppDouble newRNG(currentTrack.rngState.Branch());
 
-      Track &gamma1 = secondaries.gammas.NextTrack(newRNG, double{copcore::units::kElectronMassC2}, currentTrack.pos,
-                                                   vecgeom::Vector3D<double>{sint * cosPhi, sint * sinPhi, cost},
-                                                   currentTrack.navState, currentTrack, currentTrack.globalTime);
+      Track &gamma1 =
+          particleManager.gammas.NextTrack(newRNG, double{copcore::units::kElectronMassC2}, currentTrack.pos,
+                                           vecgeom::Vector3D<double>{sint * cosPhi, sint * sinPhi, cost},
+                                           currentTrack.navState, currentTrack, currentTrack.globalTime);
 
       // Reuse the RNG state of the dying track.
-      Track &gamma2 =
-          secondaries.gammas.NextTrack(currentTrack.rngState, double{copcore::units::kElectronMassC2}, currentTrack.pos,
-                                       -gamma1.dir, currentTrack.navState, currentTrack, currentTrack.globalTime);
+      Track &gamma2 = particleManager.gammas.NextTrack(currentTrack.rngState, double{copcore::units::kElectronMassC2},
+                                                       currentTrack.pos, -gamma1.dir, currentTrack.navState,
+                                                       currentTrack, currentTrack.globalTime);
 
       // if tracking or stepping action is called, return initial step
       if (returnLastStep) {
@@ -764,26 +727,26 @@ __device__ __forceinline__ void PerformStoppedAnnihilation(const int slot, Track
 }
 
 template <bool IsElectron, typename Scoring>
-__global__ void ElectronIonization(Track *electrons, G4HepEmElectronTrack *hepEMTracks, Secondaries secondaries,
-                                   adept::MParray *nextActiveQueue, adept::MParray *interactingQueue,
-                                   Scoring *userScoring, const bool returnAllSteps, const bool returnLastStep)
+__global__ void ElectronIonization(G4HepEmElectronTrack *hepEMTracks, ParticleManager particleManager,
+                                   adept::MParray *interactingQueue, Scoring *userScoring, const bool returnAllSteps,
+                                   const bool returnLastStep)
 {
-  int activeSize = interactingQueue->size();
+  auto &electronsOrPositrons = (IsElectron ? particleManager.electrons : particleManager.positrons);
+  SlotManager &slotManager   = *electronsOrPositrons.fSlotManager;
+  int activeSize             = interactingQueue->size();
   for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < activeSize; i += blockDim.x * gridDim.x) {
     // const int slot           = (*active)[i];
-    const int slot           = (*interactingQueue)[i];
-    SlotManager &slotManager = IsElectron ? *secondaries.electrons.fSlotManager : *secondaries.positrons.fSlotManager;
-
-    Track &currentTrack = electrons[slot];
+    const int slot      = (*interactingQueue)[i];
+    Track &currentTrack = electronsOrPositrons.TrackAt(slot);
     // the MCC vector is indexed by the logical volume id
     const int lvolID = currentTrack.navState.GetLogicalId();
 
-    VolAuxData const &auxData = AsyncAdePT::gVolAuxData[lvolID]; // FIXME unify VolAuxData
+    VolAuxData const &auxData = AsyncAdePT::gVolAuxData[lvolID];
     bool isLastStep           = true;
 
     auto survive = [&]() {
       isLastStep = false; // track survived, do not force return of step
-      nextActiveQueue->push_back(slot);
+      electronsOrPositrons.EnqueueNext(slot);
     };
 
     // Retrieve HepEM track
@@ -825,10 +788,10 @@ __global__ void ElectronIonization(Track *electrons, G4HepEmElectronTrack *hepEM
       energyDeposit += deltaEkin;
 
     } else {
-      Track &secondary =
-          secondaries.electrons.NextTrack(newRNG, deltaEkin, currentTrack.pos,
-                                          vecgeom::Vector3D<double>{dirSecondary[0], dirSecondary[1], dirSecondary[2]},
-                                          currentTrack.navState, currentTrack, currentTrack.globalTime);
+      Track &secondary = particleManager.electrons.NextTrack(
+          newRNG, deltaEkin, currentTrack.pos,
+          vecgeom::Vector3D<double>{dirSecondary[0], dirSecondary[1], dirSecondary[2]}, currentTrack.navState,
+          currentTrack, currentTrack.globalTime);
 
       // if tracking or stepping action is called, return initial step
       if (returnLastStep) {
@@ -864,7 +827,7 @@ __global__ void ElectronIonization(Track *electrons, G4HepEmElectronTrack *hepEM
         energyDeposit += currentTrack.eKin;
       }
       currentTrack.stopped = true;
-      PerformStoppedAnnihilation<IsElectron, Scoring>(slot, currentTrack, secondaries, energyDeposit, slotManager,
+      PerformStoppedAnnihilation<IsElectron, Scoring>(slot, currentTrack, particleManager, energyDeposit, slotManager,
                                                       ApplyCuts, theGammaCut, userScoring, returnLastStep);
     } else {
       currentTrack.dir.Set(dirPrimary[0], dirPrimary[1], dirPrimary[2]);
@@ -899,17 +862,17 @@ __global__ void ElectronIonization(Track *electrons, G4HepEmElectronTrack *hepEM
 }
 
 template <bool IsElectron, typename Scoring>
-__global__ void ElectronBremsstrahlung(Track *electrons, G4HepEmElectronTrack *hepEMTracks, Secondaries secondaries,
-                                       adept::MParray *nextActiveQueue, adept::MParray *interactingQueue,
-                                       Scoring *userScoring, const bool returnAllSteps, const bool returnLastStep)
+__global__ void ElectronBremsstrahlung(G4HepEmElectronTrack *hepEMTracks, ParticleManager particleManager,
+                                       adept::MParray *interactingQueue, Scoring *userScoring,
+                                       const bool returnAllSteps, const bool returnLastStep)
 {
-  int activeSize = interactingQueue->size();
+  auto &electronsOrPositrons = (IsElectron ? particleManager.electrons : particleManager.positrons);
+  SlotManager &slotManager   = *electronsOrPositrons.fSlotManager;
+  int activeSize             = interactingQueue->size();
   for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < activeSize; i += blockDim.x * gridDim.x) {
     // const int slot           = (*active)[i];
-    const int slot           = (*interactingQueue)[i];
-    SlotManager &slotManager = IsElectron ? *secondaries.electrons.fSlotManager : *secondaries.positrons.fSlotManager;
-
-    Track &currentTrack = electrons[slot];
+    const int slot      = (*interactingQueue)[i];
+    Track &currentTrack = electronsOrPositrons.TrackAt(slot);
     // the MCC vector is indexed by the logical volume id
     const int lvolID = currentTrack.navState.GetLogicalId();
 
@@ -918,7 +881,7 @@ __global__ void ElectronBremsstrahlung(Track *electrons, G4HepEmElectronTrack *h
 
     auto survive = [&]() {
       isLastStep = false; // track survived, do not force return of step
-      nextActiveQueue->push_back(slot);
+      electronsOrPositrons.EnqueueNext(slot);
     };
 
     // Retrieve HepEM track
@@ -963,9 +926,9 @@ __global__ void ElectronBremsstrahlung(Track *electrons, G4HepEmElectronTrack *h
 
     } else {
       Track &gamma =
-          secondaries.gammas.NextTrack(newRNG, deltaEkin, currentTrack.pos,
-                                       vecgeom::Vector3D<double>{dirSecondary[0], dirSecondary[1], dirSecondary[2]},
-                                       currentTrack.navState, currentTrack, currentTrack.globalTime);
+          particleManager.gammas.NextTrack(newRNG, deltaEkin, currentTrack.pos,
+                                           vecgeom::Vector3D<double>{dirSecondary[0], dirSecondary[1], dirSecondary[2]},
+                                           currentTrack.navState, currentTrack, currentTrack.globalTime);
       // if tracking or stepping action is called, return initial step
       if (returnLastStep) {
         adept_scoring::RecordHit(userScoring, gamma.trackId, gamma.parentId, /*CreatorProcessId*/ short(1),
@@ -999,7 +962,7 @@ __global__ void ElectronBremsstrahlung(Track *electrons, G4HepEmElectronTrack *h
         energyDeposit += currentTrack.eKin;
       }
       currentTrack.stopped = true;
-      PerformStoppedAnnihilation<IsElectron, Scoring>(slot, currentTrack, secondaries, energyDeposit, slotManager,
+      PerformStoppedAnnihilation<IsElectron, Scoring>(slot, currentTrack, particleManager, energyDeposit, slotManager,
                                                       ApplyCuts, theGammaCut, userScoring, returnLastStep);
     } else {
       currentTrack.dir.Set(dirPrimary[0], dirPrimary[1], dirPrimary[2]);
@@ -1034,26 +997,26 @@ __global__ void ElectronBremsstrahlung(Track *electrons, G4HepEmElectronTrack *h
 }
 
 template <typename Scoring>
-__global__ void PositronAnnihilation(Track *electrons, G4HepEmElectronTrack *hepEMTracks, Secondaries secondaries,
-                                     adept::MParray *nextActiveQueue, adept::MParray *interactingQueue,
-                                     Scoring *userScoring, const bool returnAllSteps, const bool returnLastStep)
+__global__ void PositronAnnihilation(G4HepEmElectronTrack *hepEMTracks, ParticleManager particleManager,
+                                     adept::MParray *interactingQueue, Scoring *userScoring, const bool returnAllSteps,
+                                     const bool returnLastStep)
 {
-  int activeSize = interactingQueue->size();
+  SlotManager &slotManager = *particleManager.positrons.fSlotManager;
+  int activeSize           = interactingQueue->size();
   for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < activeSize; i += blockDim.x * gridDim.x) {
     // const int slot           = (*active)[i];
-    const int slot           = (*interactingQueue)[i];
-    SlotManager &slotManager = *secondaries.positrons.fSlotManager;
+    const int slot = (*interactingQueue)[i];
 
-    Track &currentTrack = electrons[slot];
+    Track &currentTrack = particleManager.positrons.TrackAt(slot);
     // the MCC vector is indexed by the logical volume id
     const int lvolID = currentTrack.navState.GetLogicalId();
 
-    VolAuxData const &auxData = AsyncAdePT::gVolAuxData[lvolID]; // FIXME unify VolAuxData
+    VolAuxData const &auxData = AsyncAdePT::gVolAuxData[lvolID];
     bool isLastStep           = true;
 
     auto survive = [&]() {
       isLastStep = false; // track survived, do not force return of step
-      nextActiveQueue->push_back(slot);
+      particleManager.positrons.EnqueueNext(slot);
     };
 
     // Retrieve HepEM track
@@ -1095,9 +1058,9 @@ __global__ void PositronAnnihilation(Track *electrons, G4HepEmElectronTrack *hep
 
     } else {
       Track &gamma1 =
-          secondaries.gammas.NextTrack(newRNG, theGamma1Ekin, currentTrack.pos,
-                                       vecgeom::Vector3D<double>{theGamma1Dir[0], theGamma1Dir[1], theGamma1Dir[2]},
-                                       currentTrack.navState, currentTrack, currentTrack.globalTime);
+          particleManager.gammas.NextTrack(newRNG, theGamma1Ekin, currentTrack.pos,
+                                           vecgeom::Vector3D<double>{theGamma1Dir[0], theGamma1Dir[1], theGamma1Dir[2]},
+                                           currentTrack.navState, currentTrack, currentTrack.globalTime);
       // if tracking or stepping action is called, return initial step
       if (returnLastStep) {
         adept_scoring::RecordHit(userScoring, gamma1.trackId, gamma1.parentId, /*CreatorProcessId*/ short(2),
@@ -1127,9 +1090,9 @@ __global__ void PositronAnnihilation(Track *electrons, G4HepEmElectronTrack *hep
 
     } else {
       Track &gamma2 =
-          secondaries.gammas.NextTrack(currentTrack.rngState, theGamma2Ekin, currentTrack.pos,
-                                       vecgeom::Vector3D<double>{theGamma2Dir[0], theGamma2Dir[1], theGamma2Dir[2]},
-                                       currentTrack.navState, currentTrack, currentTrack.globalTime);
+          particleManager.gammas.NextTrack(currentTrack.rngState, theGamma2Ekin, currentTrack.pos,
+                                           vecgeom::Vector3D<double>{theGamma2Dir[0], theGamma2Dir[1], theGamma2Dir[2]},
+                                           currentTrack.navState, currentTrack, currentTrack.globalTime);
       // if tracking or stepping action is called, return initial step
       if (returnLastStep) {
         adept_scoring::RecordHit(userScoring, gamma2.trackId, gamma2.parentId, /*CreatorProcessId*/ short(2),
@@ -1185,18 +1148,17 @@ __global__ void PositronAnnihilation(Track *electrons, G4HepEmElectronTrack *hep
 }
 
 template <typename Scoring>
-__global__ void PositronStoppedAnnihilation(Track *electrons, G4HepEmElectronTrack *hepEMTracks,
-                                            Secondaries secondaries, adept::MParray *nextActiveQueue,
+__global__ void PositronStoppedAnnihilation(G4HepEmElectronTrack *hepEMTracks, ParticleManager particleManager,
                                             adept::MParray *interactingQueue, Scoring *userScoring,
                                             const bool returnAllSteps, const bool returnLastStep)
 {
-  int activeSize = interactingQueue->size();
+  SlotManager &slotManager = *particleManager.positrons.fSlotManager;
+  int activeSize           = interactingQueue->size();
   for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < activeSize; i += blockDim.x * gridDim.x) {
     // const int slot           = (*active)[i];
-    const int slot           = (*interactingQueue)[i];
-    SlotManager &slotManager = *secondaries.positrons.fSlotManager;
+    const int slot = (*interactingQueue)[i];
 
-    Track &currentTrack = electrons[slot];
+    Track &currentTrack = particleManager.positrons.TrackAt(slot);
     // the MCC vector is indexed by the logical volume id
     const int lvolID = currentTrack.navState.GetLogicalId();
 
@@ -1205,7 +1167,7 @@ __global__ void PositronStoppedAnnihilation(Track *electrons, G4HepEmElectronTra
 
     auto survive = [&]() {
       isLastStep = false; // track survived, do not force return of step
-      nextActiveQueue->push_back(slot);
+      particleManager.positrons.EnqueueNext(slot);
     };
 
     // Retrieve HepEM track
@@ -1225,8 +1187,8 @@ __global__ void PositronStoppedAnnihilation(Track *electrons, G4HepEmElectronTra
     // Annihilate the stopped positron into two gammas heading to opposite
     // directions (isotropic).
 
-    PerformStoppedAnnihilation<false, Scoring>(slot, currentTrack, secondaries, energyDeposit, slotManager, ApplyCuts,
-                                               theGammaCut, userScoring, returnLastStep);
+    PerformStoppedAnnihilation<false, Scoring>(slot, currentTrack, particleManager, energyDeposit, slotManager,
+                                               ApplyCuts, theGammaCut, userScoring, returnLastStep);
 
     // Record the step. Edep includes the continuous energy loss and edep from secondaries which were cut
     if ((energyDeposit > 0 && auxData.fSensIndex >= 0) || returnAllSteps || returnLastStep)
diff --git a/include/AdePT/kernels/gammas.cuh b/include/AdePT/kernels/gammas.cuh
index bddab945..6339a111 100644
--- a/include/AdePT/kernels/gammas.cuh
+++ b/include/AdePT/kernels/gammas.cuh
@@ -26,18 +26,17 @@ namespace AsyncAdePT {
 // Asynchronous TransportGammas Interface
 template <typename Scoring, class SteppingActionT>
 __global__ void __launch_bounds__(256, 1)
-    TransportGammas(Track *gammas, Track *leaks, const adept::MParray *active, Secondaries secondaries,
-                    adept::MParray *nextActiveQueue, adept::MParray *leakedQueue, Scoring *userScoring,
-                    Stats *InFlightStats, const StepActionParam params, AllowFinishOffEventArray allowFinishOffEvent,
+    TransportGammas(ParticleManager particleManager, Scoring *userScoring, Stats *InFlightStats,
+                    const StepActionParam params, AllowFinishOffEventArray allowFinishOffEvent,
                     const bool returnAllSteps, const bool returnLastStep)
 {
   constexpr double kPushDistance    = 1000 * vecgeom::kTolerance;
   constexpr unsigned short maxSteps = 10'000;
-  int activeSize                    = active->size();
+  auto &slotManager                 = *particleManager.gammas.fSlotManager;
+  const int activeSize              = particleManager.gammas.ActiveSize();
   for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < activeSize; i += blockDim.x * gridDim.x) {
-    const int slot      = (*active)[i];
-    auto &slotManager   = *secondaries.gammas.fSlotManager;
-    Track &currentTrack = gammas[slot];
+    const auto slot     = particleManager.gammas.ActiveAt(i);
+    Track &currentTrack = particleManager.gammas.TrackAt(slot);
 #else
 
 // Synchronous TransportGammas Interface
@@ -109,21 +108,10 @@ __global__ void TransportGammas(adept::TrackManager<Track> *gammas, Secondaries
       currentTrack.leakStatus = leakReason;
 #ifdef ASYNC_MODE
       if (leakReason != LeakStatus::NoLeak) {
-        // Get a slot in the leaks array
-        int leakSlot = secondaries.gammas.NextLeakSlot();
-        // Copy the track to the leaks array and store the index in the leak queue
-        leaks[leakSlot] = gammas[slot];
-        auto success    = leakedQueue->push_back(leakSlot);
-        if (!success) {
-          printf("ERROR: No space left in gammas leaks queue.\n\
-\tThe threshold for flushing the leak buffer may be too high\n\
-\tThe space allocated to the leak buffer may be too small\n");
-          asm("trap;");
-        }
-        // Free the slot in the tracks slot manager
-        slotManager.MarkSlotForFreeing(slot);
+        // Copy track at slot to the leaked tracks
+        particleManager.gammas.CopyTrackToLeaked(slot);
       } else {
-        nextActiveQueue->push_back(slot);
+        particleManager.gammas.EnqueueNext(slot);
       }
 #else
       currentTrack.CopyTo(trackdata, Pdg);
@@ -311,7 +299,7 @@ __global__ void TransportGammas(adept::TrackManager<Track> *gammas, Secondaries
           edep = elKinEnergy;
         } else {
 #ifdef ASYNC_MODE
-          Track &electron = secondaries.electrons.NextTrack(
+          Track &electron = particleManager.electrons.NextTrack(
               newRNG, elKinEnergy, pos,
               vecgeom::Vector3D<double>{dirSecondaryEl[0], dirSecondaryEl[1], dirSecondaryEl[2]}, navState,
               currentTrack, globalTime);
@@ -355,7 +343,7 @@ __global__ void TransportGammas(adept::TrackManager<Track> *gammas, Secondaries
           edep += posKinEnergy + 2 * copcore::units::kElectronMassC2;
         } else {
 #ifdef ASYNC_MODE
-          Track &positron = secondaries.positrons.NextTrack(
+          Track &positron = particleManager.positrons.NextTrack(
               currentTrack.rngState, posKinEnergy, pos,
               vecgeom::Vector3D<double>{dirSecondaryPos[0], dirSecondaryPos[1], dirSecondaryPos[2]}, navState,
               currentTrack, globalTime);
@@ -423,7 +411,7 @@ __global__ void TransportGammas(adept::TrackManager<Track> *gammas, Secondaries
         if (ApplyCuts ? energyEl > theElCut : energyEl > LowEnergyThreshold) {
           // Create a secondary electron and sample/compute directions.
 #ifdef ASYNC_MODE
-          Track &electron = secondaries.electrons.NextTrack(
+          Track &electron = particleManager.electrons.NextTrack(
               newRNG, energyEl, pos, eKin * dir - newEnergyGamma * newDirGamma, navState, currentTrack, globalTime);
 #else
           Track &electron = secondaries.electrons->NextTrack();
@@ -501,7 +489,7 @@ __global__ void TransportGammas(adept::TrackManager<Track> *gammas, Secondaries
 
           // Create a secondary electron and sample directions.
 #ifdef ASYNC_MODE
-          Track &electron = secondaries.electrons.NextTrack(
+          Track &electron = particleManager.electrons.NextTrack(
               newRNG, photoElecE, pos, vecgeom::Vector3D<double>{dirPhotoElec[0], dirPhotoElec[1], dirPhotoElec[2]},
               navState, currentTrack, globalTime);
 #else
diff --git a/include/AdePT/kernels/gammas_split.cuh b/include/AdePT/kernels/gammas_split.cuh
index cef659c5..23225d1c 100644
--- a/include/AdePT/kernels/gammas_split.cuh
+++ b/include/AdePT/kernels/gammas_split.cuh
@@ -21,20 +21,19 @@ using VolAuxData = adeptint::VolAuxData;
 
 namespace AsyncAdePT {
 
-__global__ void GammaHowFar(Track *gammas, Track *leaks, G4HepEmGammaTrack *hepEMTracks, const adept::MParray *active,
-                            Secondaries secondaries, adept::MParray *propagationQueue, adept::MParray *leakedQueue,
-                            Stats *InFlightStats, AllowFinishOffEventArray allowFinishOffEvent)
+__global__ void GammaHowFar(G4HepEmGammaTrack *hepEMTracks, ParticleManager particleManager,
+                            adept::MParray *propagationQueue, Stats *InFlightStats,
+                            AllowFinishOffEventArray allowFinishOffEvent)
 {
   constexpr unsigned short maxSteps        = 10'000;
   constexpr unsigned short kStepsStuckKill = 25;
-  int activeSize                           = active->size();
+  const int activeSize                     = particleManager.gammas.ActiveSize();
   for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < activeSize; i += blockDim.x * gridDim.x) {
-    const int slot      = (*active)[i];
-    auto &slotManager   = *secondaries.gammas.fSlotManager;
-    Track &currentTrack = gammas[slot];
+    const auto slot     = particleManager.gammas.ActiveAt(i);
+    Track &currentTrack = particleManager.gammas.TrackAt(slot);
 
     int lvolID                = currentTrack.navState.GetLogicalId();
-    VolAuxData const &auxData = AsyncAdePT::gVolAuxData[lvolID]; // FIXME unify VolAuxData
+    VolAuxData const &auxData = AsyncAdePT::gVolAuxData[lvolID];
 
     currentTrack.preStepEKin       = currentTrack.eKin;
     currentTrack.preStepGlobalTime = currentTrack.globalTime;
@@ -57,19 +56,8 @@ __global__ void GammaHowFar(Track *gammas, Track *leaks, G4HepEmGammaTrack *hepE
     auto survive = [&](LeakStatus leakReason = LeakStatus::NoLeak) {
       currentTrack.leakStatus = leakReason;
       if (leakReason != LeakStatus::NoLeak) {
-        // Get a slot in the leaks array
-        int leakSlot = secondaries.gammas.NextLeakSlot();
-        // Copy the track to the leaks array and store the index in the leak queue
-        leaks[leakSlot] = gammas[slot];
-        auto success    = leakedQueue->push_back(leakSlot);
-        if (!success) {
-          printf("ERROR: No space left in gammas leaks queue.\n\
-\tThe threshold for flushing the leak buffer may be too high\n\
-\tThe space allocated to the leak buffer may be too small\n");
-          asm("trap;");
-        }
-        // Free the slot in the tracks slot manager
-        slotManager.MarkSlotForFreeing(slot);
+        // Copy track at slot to the leaked tracks
+        particleManager.gammas.CopyTrackToLeaked(slot);
       }
     };
 
@@ -161,16 +149,14 @@ __global__ void GammaPropagation(Track *gammas, G4HepEmGammaTrack *hepEMTracks,
 }
 
 template <typename Scoring>
-__global__ void GammaSetupInteractions(Track *gammas, Track *leaks, G4HepEmGammaTrack *hepEMTracks,
-                                       const adept::MParray *active, Secondaries secondaries,
-                                       adept::MParray *nextActiveQueue, AllInteractionQueues interactionQueues,
-                                       adept::MParray *leakedQueue, Scoring *userScoring)
+__global__ void GammaSetupInteractions(G4HepEmGammaTrack *hepEMTracks, const adept::MParray *propagationQueue,
+                                       ParticleManager particleManager, AllInteractionQueues interactionQueues,
+                                       Scoring *userScoring)
 {
-  int activeSize = active->size();
+  int activeSize = propagationQueue->size();
   for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < activeSize; i += blockDim.x * gridDim.x) {
-    const int slot      = (*active)[i];
-    auto &slotManager   = *secondaries.gammas.fSlotManager;
-    Track &currentTrack = gammas[slot];
+    const int slot      = (*propagationQueue)[i];
+    Track &currentTrack = particleManager.gammas.TrackAt(slot);
 
     int lvolID = currentTrack.navState.GetLogicalId();
 
@@ -178,21 +164,10 @@ __global__ void GammaSetupInteractions(Track *gammas, Track *leaks, G4HepEmGamma
     auto survive = [&](LeakStatus leakReason = LeakStatus::NoLeak) {
       currentTrack.leakStatus = leakReason;
       if (leakReason != LeakStatus::NoLeak) {
-        // Get a slot in the leaks array
-        int leakSlot = secondaries.gammas.NextLeakSlot();
-        // Copy the track to the leaks array and store the index in the leak queue
-        leaks[leakSlot] = gammas[slot];
-        auto success    = leakedQueue->push_back(leakSlot);
-        if (!success) {
-          printf("ERROR: No space left in gammas leaks queue.\n\
-\tThe threshold for flushing the leak buffer may be too high\n\
-\tThe space allocated to the leak buffer may be too small\n");
-          asm("trap;");
-        }
-        // Free the slot in the tracks slot manager
-        slotManager.MarkSlotForFreeing(slot);
+        // Copy track at slot to the leaked tracks
+        particleManager.gammas.CopyTrackToLeaked(slot);
       } else {
-        nextActiveQueue->push_back(slot);
+        particleManager.gammas.EnqueueNext(slot);
       }
     };
 
@@ -227,17 +202,16 @@ __global__ void GammaSetupInteractions(Track *gammas, Track *leaks, G4HepEmGamma
 }
 
 template <typename Scoring>
-__global__ void GammaRelocation(Track *gammas, Track *leaks, G4HepEmGammaTrack *hepEMTracks, Secondaries secondaries,
-                                adept::MParray *nextActiveQueue, adept::MParray *relocatingQueue,
-                                adept::MParray *leakedQueue, Scoring *userScoring, const bool returnAllSteps,
+__global__ void GammaRelocation(G4HepEmGammaTrack *hepEMTracks, ParticleManager particleManager,
+                                adept::MParray *relocatingQueue, Scoring *userScoring, const bool returnAllSteps,
                                 const bool returnLastStep)
 {
   constexpr double kPushDistance = 1000 * vecgeom::kTolerance;
+  auto &slotManager              = *particleManager.gammas.fSlotManager;
   int activeSize                 = relocatingQueue->size();
   for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < activeSize; i += blockDim.x * gridDim.x) {
     const int slot      = (*relocatingQueue)[i];
-    auto &slotManager   = *secondaries.gammas.fSlotManager;
-    Track &currentTrack = gammas[slot];
+    Track &currentTrack = particleManager.gammas.TrackAt(slot);
 
     int lvolID = currentTrack.navState.GetLogicalId();
 
@@ -245,21 +219,10 @@ __global__ void GammaRelocation(Track *gammas, Track *leaks, G4HepEmGammaTrack *
     auto survive = [&](LeakStatus leakReason = LeakStatus::NoLeak) {
       currentTrack.leakStatus = leakReason;
       if (leakReason != LeakStatus::NoLeak) {
-        // Get a slot in the leaks array
-        int leakSlot = secondaries.gammas.NextLeakSlot();
-        // Copy the track to the leaks array and store the index in the leak queue
-        leaks[leakSlot] = gammas[slot];
-        auto success    = leakedQueue->push_back(leakSlot);
-        if (!success) {
-          printf("ERROR: No space left in gammas leaks queue.\n\
-\tThe threshold for flushing the leak buffer may be too high\n\
-\tThe space allocated to the leak buffer may be too small\n");
-          asm("trap;");
-        }
-        // Free the slot in the tracks slot manager
-        slotManager.MarkSlotForFreeing(slot);
+        // Copy track at slot to the leaked tracks
+        particleManager.gammas.CopyTrackToLeaked(slot);
       } else {
-        nextActiveQueue->push_back(slot);
+        particleManager.gammas.EnqueueNext(slot);
       }
     };
 
@@ -356,17 +319,15 @@ __global__ void GammaRelocation(Track *gammas, Track *leaks, G4HepEmGammaTrack *
 }
 
 template <typename Scoring>
-__global__ void GammaConversion(Track *gammas, G4HepEmGammaTrack *hepEMTracks, Secondaries secondaries,
-                                adept::MParray *nextActiveQueue, adept::MParray *interactingQueue, Scoring *userScoring,
-                                const bool returnAllSteps, const bool returnLastStep)
+__global__ void GammaConversion(G4HepEmGammaTrack *hepEMTracks, ParticleManager particleManager,
+                                adept::MParray *interactingQueue, Scoring *userScoring, const bool returnAllSteps,
+                                const bool returnLastStep)
 {
-  // int activeSize = active->size();
-  int activeSize = interactingQueue->size();
+  auto &slotManager = *particleManager.gammas.fSlotManager;
+  int activeSize    = interactingQueue->size();
   for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < activeSize; i += blockDim.x * gridDim.x) {
-    // const int slot      = (*active)[i];
     const int slot      = (*interactingQueue)[i];
-    auto &slotManager   = *secondaries.gammas.fSlotManager;
-    Track &currentTrack = gammas[slot];
+    Track &currentTrack = particleManager.gammas.TrackAt(slot);
 
     int lvolID                = currentTrack.navState.GetLogicalId();
     VolAuxData const &auxData = AsyncAdePT::gVolAuxData[lvolID]; // FIXME unify VolAuxData
@@ -375,7 +336,7 @@ __global__ void GammaConversion(Track *gammas, G4HepEmGammaTrack *hepEMTracks, S
     // Write local variables back into track and enqueue
     auto survive = [&]() {
       isLastStep = false; // particle survived
-      nextActiveQueue->push_back(slot);
+      particleManager.gammas.EnqueueNext(slot);
     };
 
     G4HepEmGammaTrack &gammaTrack = hepEMTracks[slot];
@@ -421,7 +382,7 @@ __global__ void GammaConversion(Track *gammas, G4HepEmGammaTrack *hepEMTracks, S
       // Deposit the energy here and kill the secondary
       edep = elKinEnergy;
     } else {
-      Track &electron = secondaries.electrons.NextTrack(
+      Track &electron = particleManager.electrons.NextTrack(
           newRNG, elKinEnergy, currentTrack.pos,
           vecgeom::Vector3D<double>{dirSecondaryEl[0], dirSecondaryEl[1], dirSecondaryEl[2]}, currentTrack.navState,
           currentTrack, currentTrack.globalTime);
@@ -455,7 +416,7 @@ __global__ void GammaConversion(Track *gammas, G4HepEmGammaTrack *hepEMTracks, S
       // Deposit: posKinEnergy + 2 * copcore::units::kElectronMassC2 and kill the secondary
       edep += posKinEnergy + 2 * copcore::units::kElectronMassC2;
     } else {
-      Track &positron = secondaries.positrons.NextTrack(
+      Track &positron = particleManager.positrons.NextTrack(
           currentTrack.rngState, posKinEnergy, currentTrack.pos,
           vecgeom::Vector3D<double>{dirSecondaryPos[0], dirSecondaryPos[1], dirSecondaryPos[2]}, currentTrack.navState,
           currentTrack, currentTrack.globalTime);
@@ -519,26 +480,24 @@ __global__ void GammaConversion(Track *gammas, G4HepEmGammaTrack *hepEMTracks, S
 }
 
 template <typename Scoring>
-__global__ void GammaCompton(Track *gammas, G4HepEmGammaTrack *hepEMTracks, Secondaries secondaries,
-                             adept::MParray *nextActiveQueue, adept::MParray *interactingQueue, Scoring *userScoring,
-                             const bool returnAllSteps, const bool returnLastStep)
+__global__ void GammaCompton(G4HepEmGammaTrack *hepEMTracks, ParticleManager particleManager,
+                             adept::MParray *interactingQueue, Scoring *userScoring, const bool returnAllSteps,
+                             const bool returnLastStep)
 {
-  // int activeSize = active->size();
-  int activeSize = interactingQueue->size();
+  auto &slotManager = *particleManager.gammas.fSlotManager;
+  int activeSize    = interactingQueue->size();
   for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < activeSize; i += blockDim.x * gridDim.x) {
-    // const int slot      = (*active)[i];
     const int slot      = (*interactingQueue)[i];
-    auto &slotManager   = *secondaries.gammas.fSlotManager;
-    Track &currentTrack = gammas[slot];
+    Track &currentTrack = particleManager.gammas.TrackAt(slot);
 
     int lvolID                = currentTrack.navState.GetLogicalId();
-    VolAuxData const &auxData = AsyncAdePT::gVolAuxData[lvolID]; // FIXME unify VolAuxData
+    VolAuxData const &auxData = AsyncAdePT::gVolAuxData[lvolID];
     bool isLastStep           = true;
 
     // Write local variables back into track and enqueue
     auto survive = [&]() {
       isLastStep = false; // particle survived
-      nextActiveQueue->push_back(slot);
+      particleManager.gammas.EnqueueNext(slot);
     };
 
     G4HepEmGammaTrack &gammaTrack = hepEMTracks[slot];
@@ -578,7 +537,7 @@ __global__ void GammaCompton(Track *gammas, G4HepEmGammaTrack *hepEMTracks, Seco
     // Check the cuts and deposit energy in this volume if needed
     if (ApplyCuts ? energyEl > theElCut : energyEl > LowEnergyThreshold) {
       // Create a secondary electron and sample/compute directions.
-      Track &electron = secondaries.electrons.NextTrack(
+      Track &electron = particleManager.electrons.NextTrack(
           newRNG, energyEl, currentTrack.pos, currentTrack.eKin * currentTrack.dir - newEnergyGamma * newDirGamma,
           currentTrack.navState, currentTrack, currentTrack.globalTime);
       electron.dir.Normalize();
@@ -656,17 +615,15 @@ __global__ void GammaCompton(Track *gammas, G4HepEmGammaTrack *hepEMTracks, Seco
 }
 
 template <typename Scoring>
-__global__ void GammaPhotoelectric(Track *gammas, G4HepEmGammaTrack *hepEMTracks, Secondaries secondaries,
-                                   adept::MParray *nextActiveQueue, adept::MParray *interactingQueue,
-                                   Scoring *userScoring, const bool returnAllSteps, const bool returnLastStep)
+__global__ void GammaPhotoelectric(G4HepEmGammaTrack *hepEMTracks, ParticleManager particleManager,
+                                   adept::MParray *interactingQueue, Scoring *userScoring, const bool returnAllSteps,
+                                   const bool returnLastStep)
 {
-  // int activeSize = active->size();
-  int activeSize = interactingQueue->size();
+  auto &slotManager = *particleManager.gammas.fSlotManager;
+  int activeSize    = interactingQueue->size();
   for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < activeSize; i += blockDim.x * gridDim.x) {
-    // const int slot      = (*active)[i];
     const int slot      = (*interactingQueue)[i];
-    auto &slotManager   = *secondaries.gammas.fSlotManager;
-    Track &currentTrack = gammas[slot];
+    Track &currentTrack = particleManager.gammas.TrackAt(slot);
 
     int lvolID                = currentTrack.navState.GetLogicalId();
     VolAuxData const &auxData = AsyncAdePT::gVolAuxData[lvolID]; // FIXME unify VolAuxData
@@ -675,7 +632,7 @@ __global__ void GammaPhotoelectric(Track *gammas, G4HepEmGammaTrack *hepEMTracks
     // Write local variables back into track and enqueue
     auto survive = [&]() {
       isLastStep = false; // particle survived
-      nextActiveQueue->push_back(slot);
+      particleManager.gammas.EnqueueNext(slot);
     };
 
     G4HepEmGammaTrack &gammaTrack = hepEMTracks[slot];
@@ -713,10 +670,10 @@ __global__ void GammaPhotoelectric(Track *gammas, G4HepEmGammaTrack *hepEMTracks
       G4HepEmGammaInteractionPhotoelectric::SamplePhotoElectronDirection(photoElecE, dirGamma, dirPhotoElec, &rnge);
 
       // Create a secondary electron and sample directions.
-      Track &electron =
-          secondaries.electrons.NextTrack(newRNG, photoElecE, currentTrack.pos,
-                                          vecgeom::Vector3D<double>{dirPhotoElec[0], dirPhotoElec[1], dirPhotoElec[2]},
-                                          currentTrack.navState, currentTrack, currentTrack.globalTime);
+      Track &electron = particleManager.electrons.NextTrack(
+          newRNG, photoElecE, currentTrack.pos,
+          vecgeom::Vector3D<double>{dirPhotoElec[0], dirPhotoElec[1], dirPhotoElec[2]}, currentTrack.navState,
+          currentTrack, currentTrack.globalTime);
 
       // if tracking or stepping action is called, return initial step
       if (returnLastStep) {