diff --git a/src/struphy/kernel_arguments/pusher_args.cuh b/src/struphy/kernel_arguments/pusher_args.cuh
index 5c85f3d45..59815fced 100644
--- a/src/struphy/kernel_arguments/pusher_args.cuh
+++ b/src/struphy/kernel_arguments/pusher_args.cuh
@@ -12,7 +12,7 @@ struct MarkerArgs {
     int vdim;
     int weight_idx;
     int first_diagnostics_idx;
-    int first_init_idx;
+    int first_pusher_idx;
     int first_shift_idx;
     int residual_idx;
     int first_free_idx;
diff --git a/src/struphy/kernel_arguments/pusher_args_cuda.py b/src/struphy/kernel_arguments/pusher_args_cuda.py
index 0c9911d72..08fefff46 100644
--- a/src/struphy/kernel_arguments/pusher_args_cuda.py
+++ b/src/struphy/kernel_arguments/pusher_args_cuda.py
@@ -39,7 +39,7 @@ class CudaMarkerArguments(CudaStructArguments):
         ("vdim", "int"),
         ("weight_idx", "int"),
         ("first_diagnostics_idx", "int"),
-        ("first_init_idx", "int"),
+        ("first_pusher_idx", "int"),
         ("first_shift_idx", "int"),
         ("residual_idx", "int"),
         ("first_free_idx", "int"),
@@ -71,7 +71,7 @@ def __init__(
         self.weight_idx = weight_idx
         self.n_markers = markers.shape[0]
         self.first_diagnostics_idx = first_diagnostics_idx
-        self.first_init_idx = first_pusher_idx
+        self.first_pusher_idx = first_pusher_idx
         self.first_shift_idx = first_shift_idx
         self.residual_idx = residual_idx
         self.first_free_idx = first_free_idx
diff --git a/src/struphy/kernel_arguments/pusher_args_kernels.py b/src/struphy/kernel_arguments/pusher_args_kernels.py
index c5b3f52d8..2e66dc58e 100644
--- a/src/struphy/kernel_arguments/pusher_args_kernels.py
+++ b/src/struphy/kernel_arguments/pusher_args_kernels.py
@@ -71,7 +71,7 @@ def __init__(
 
         # useful indices
         self.first_diagnostics_idx = first_diagnostics_idx
-        self.first_init_idx = first_pusher_idx
+        self.first_pusher_idx = first_pusher_idx
         self.first_shift_idx = first_shift_idx  # starting idx for eta-shifts due to boundary conditions
         self.residual_idx = residual_idx  # residual in iterative solvers
         self.first_free_idx = first_free_idx  # index after which auxiliary saving is possible
diff --git a/src/struphy/pic/accumulation/kernels/cc_lin_mhd_5d_gradB/cc_lin_mhd_5d_gradB_kernels.py b/src/struphy/pic/accumulation/kernels/cc_lin_mhd_5d_gradB/cc_lin_mhd_5d_gradB_kernels.py
index 4428c130d..e50a5e7a4 100644
--- a/src/struphy/pic/accumulation/kernels/cc_lin_mhd_5d_gradB/cc_lin_mhd_5d_gradB_kernels.py
+++ b/src/struphy/pic/accumulation/kernels/cc_lin_mhd_5d_gradB/cc_lin_mhd_5d_gradB_kernels.py
@@ -99,7 +99,7 @@ def cc_lin_mhd_5d_gradB(
 
     markers = args_markers.markers
     n_markers = args_markers.n_markers
-    first_init_idx = args_markers.first_init_idx
+    first_pusher_idx = args_markers.first_pusher_idx
     mu_idx = args_markers.mu_idx
 
     # allocate for magnetic field evaluation
@@ -130,7 +130,7 @@ def cc_lin_mhd_5d_gradB(
             continue
 
         # if particle is refilled
-        if markers[ip, first_init_idx] == -1.0:
+        if markers[ip, first_pusher_idx] == -1.0:
             continue
 
         # marker positions
diff --git a/src/struphy/pic/accumulation/kernels/cc_lin_mhd_5d_gradB_dg/cc_lin_mhd_5d_gradB_dg_kernels.py b/src/struphy/pic/accumulation/kernels/cc_lin_mhd_5d_gradB_dg/cc_lin_mhd_5d_gradB_dg_kernels.py
index 7705768f5..015b4506a 100644
--- a/src/struphy/pic/accumulation/kernels/cc_lin_mhd_5d_gradB_dg/cc_lin_mhd_5d_gradB_dg_kernels.py
+++ b/src/struphy/pic/accumulation/kernels/cc_lin_mhd_5d_gradB_dg/cc_lin_mhd_5d_gradB_dg_kernels.py
@@ -71,7 +71,7 @@ def cc_lin_mhd_5d_gradB_dg(
 
     markers = args_markers.markers
     mu_idx = args_markers.mu_idx
-    first_init_idx = args_markers.first_init_idx
+    first_pusher_idx = args_markers.first_pusher_idx
 
     # allocate for magnetic field evaluation
     eta_diff = empty(3, dtype=float)
@@ -108,10 +108,10 @@ def cc_lin_mhd_5d_gradB_dg(
             continue
 
         # marker positions, mid point
-        eta_mid[:] = (markers[ip, 0:3] + markers[ip, first_init_idx : first_init_idx + 3]) / 2.0
+        eta_mid[:] = (markers[ip, 0:3] + markers[ip, first_pusher_idx : first_pusher_idx + 3]) / 2.0
         eta_mid[:] = mod(eta_mid[:], 1.0)
 
-        eta_diff[:] = markers[ip, 0:3] - markers[ip, first_init_idx : first_init_idx + 3]
+        eta_diff[:] = markers[ip, 0:3] - markers[ip, first_pusher_idx : first_pusher_idx + 3]
 
         # marker weight and velocity
         weight = markers[ip, 5]
diff --git a/src/struphy/pic/diagnostics/kernels/eval_gradB_ediff/eval_gradB_ediff_kernels.py b/src/struphy/pic/diagnostics/kernels/eval_gradB_ediff/eval_gradB_ediff_kernels.py
index c21807b7c..e33183e72 100644
--- a/src/struphy/pic/diagnostics/kernels/eval_gradB_ediff/eval_gradB_ediff_kernels.py
+++ b/src/struphy/pic/diagnostics/kernels/eval_gradB_ediff/eval_gradB_ediff_kernels.py
@@ -42,7 +42,7 @@ def eval_gradB_ediff(
     markers = args_markers.markers
     n_markers = args_markers.n_markers
     mu_idx = args_markers.mu_idx
-    first_init_idx = args_markers.first_init_idx
+    first_pusher_idx = args_markers.first_pusher_idx
     first_free_idx = args_markers.first_free_idx
 
     for ip in range(n_markers):
@@ -51,10 +51,10 @@ def eval_gradB_ediff(
             continue
 
         # marker positions, mid point
-        eta_mid[:] = (markers[ip, 0:3] + markers[ip, first_init_idx : first_init_idx + 3]) / 2.0
+        eta_mid[:] = (markers[ip, 0:3] + markers[ip, first_pusher_idx : first_pusher_idx + 3]) / 2.0
         eta_mid[:] = mod(eta_mid[:], 1.0)
 
-        eta_diff = markers[ip, 0:3] - markers[ip, first_init_idx : first_init_idx + 3]
+        eta_diff = markers[ip, 0:3] - markers[ip, first_pusher_idx : first_pusher_idx + 3]
 
         # marker weight and velocity
         weight = markers[ip, 5]
diff --git a/src/struphy/pic/pushing/kernels/bstar_2form/bstar_2form_kernels.py b/src/struphy/pic/pushing/kernels/bstar_2form/bstar_2form_kernels.py
index ea1ba6ad6..42bed58ef 100644
--- a/src/struphy/pic/pushing/kernels/bstar_2form/bstar_2form_kernels.py
+++ b/src/struphy/pic/pushing/kernels/bstar_2form/bstar_2form_kernels.py
@@ -50,7 +50,7 @@ def bstar_2form(
     markers = args_markers.markers
     n_markers = args_markers.n_markers
     mu_idx = args_markers.mu_idx
-    first_init_idx = args_markers.first_init_idx
+    first_pusher_idx = args_markers.first_pusher_idx
     first_shift_idx = args_markers.first_shift_idx
 
     for ip in range(n_markers):
@@ -59,13 +59,13 @@ def bstar_2form(
             continue
 
         eta_k[:] = markers[ip, 0:3] + markers[ip, first_shift_idx : first_shift_idx + 3]
-        eta_n[:] = markers[ip, first_init_idx : first_init_idx + 3]
+        eta_n[:] = markers[ip, first_pusher_idx : first_pusher_idx + 3]
 
         eta[:] = alpha[:3] * eta_k + (1.0 - alpha[:3]) * eta_n
         eta[:] = mod(eta, 1.0)
 
         v_k = markers[ip, 3]
-        v_n = markers[ip, first_init_idx + 3]
+        v_n = markers[ip, first_pusher_idx + 3]
         v = alpha[3] * v_k + (1.0 - alpha[3]) * v_n
 
         # spline evaluation
diff --git a/src/struphy/pic/pushing/kernels/bstar_parallel_3form/bstar_parallel_3form_kernels.py b/src/struphy/pic/pushing/kernels/bstar_parallel_3form/bstar_parallel_3form_kernels.py
index f72014dd2..b6356a1b0 100644
--- a/src/struphy/pic/pushing/kernels/bstar_parallel_3form/bstar_parallel_3form_kernels.py
+++ b/src/struphy/pic/pushing/kernels/bstar_parallel_3form/bstar_parallel_3form_kernels.py
@@ -45,7 +45,7 @@ def bstar_parallel_3form(
     markers = args_markers.markers
     n_markers = args_markers.n_markers
     mu_idx = args_markers.mu_idx
-    first_init_idx = args_markers.first_init_idx
+    first_pusher_idx = args_markers.first_pusher_idx
     first_shift_idx = args_markers.first_shift_idx
 
     for ip in range(n_markers):
@@ -54,13 +54,13 @@ def bstar_parallel_3form(
             continue
 
         eta_k[:] = markers[ip, 0:3] + markers[ip, first_shift_idx : first_shift_idx + 3]
-        eta_n[:] = markers[ip, first_init_idx : first_init_idx + 3]
+        eta_n[:] = markers[ip, first_pusher_idx : first_pusher_idx + 3]
 
         eta[:] = alpha[:3] * eta_k + (1.0 - alpha[:3]) * eta_n
         eta[:] = mod(eta, 1.0)
 
         v_k = markers[ip, 3]
-        v_n = markers[ip, first_init_idx + 3]
+        v_n = markers[ip, first_pusher_idx + 3]
         v = alpha[3] * v_k + (1.0 - alpha[3]) * v_n
 
         # evaluate Jacobian, result in dfm
diff --git a/src/struphy/pic/pushing/kernels/driftkinetic_hamiltonian/driftkinetic_hamiltonian_kernels.py b/src/struphy/pic/pushing/kernels/driftkinetic_hamiltonian/driftkinetic_hamiltonian_kernels.py
index 26acb7d20..d622133b3 100644
--- a/src/struphy/pic/pushing/kernels/driftkinetic_hamiltonian/driftkinetic_hamiltonian_kernels.py
+++ b/src/struphy/pic/pushing/kernels/driftkinetic_hamiltonian/driftkinetic_hamiltonian_kernels.py
@@ -44,7 +44,7 @@ def driftkinetic_hamiltonian(
     markers = args_markers.markers
     n_markers = args_markers.n_markers
     mu_idx = args_markers.mu_idx
-    first_init_idx = args_markers.first_init_idx
+    first_pusher_idx = args_markers.first_pusher_idx
     first_shift_idx = args_markers.first_shift_idx
 
     for ip in range(n_markers):
@@ -53,13 +53,13 @@ def driftkinetic_hamiltonian(
             continue
 
         eta_k[:] = markers[ip, 0:3] + markers[ip, first_shift_idx : first_shift_idx + 3]
-        eta_n[:] = markers[ip, first_init_idx : first_init_idx + 3]
+        eta_n[:] = markers[ip, first_pusher_idx : first_pusher_idx + 3]
 
         eta[:] = alpha[:3] * eta_k + (1.0 - alpha[:3]) * eta_n
         eta[:] = mod(eta, 1.0)
 
         v_k = markers[ip, 3]
-        v_n = markers[ip, first_init_idx + 3]
+        v_n = markers[ip, first_pusher_idx + 3]
         v = alpha[3] * v_k + (1.0 - alpha[3]) * v_n
 
         mu = markers[ip, mu_idx]
diff --git a/src/struphy/pic/pushing/kernels/grad_driftkinetic_hamiltonian/grad_driftkinetic_hamiltonian_kernels.py b/src/struphy/pic/pushing/kernels/grad_driftkinetic_hamiltonian/grad_driftkinetic_hamiltonian_kernels.py
index 1e89b6e36..417d384e7 100644
--- a/src/struphy/pic/pushing/kernels/grad_driftkinetic_hamiltonian/grad_driftkinetic_hamiltonian_kernels.py
+++ b/src/struphy/pic/pushing/kernels/grad_driftkinetic_hamiltonian/grad_driftkinetic_hamiltonian_kernels.py
@@ -58,7 +58,7 @@ def grad_driftkinetic_hamiltonian(
     markers = args_markers.markers
     n_markers = args_markers.n_markers
     mu_idx = args_markers.mu_idx
-    first_init_idx = args_markers.first_init_idx
+    first_pusher_idx = args_markers.first_pusher_idx
     first_shift_idx = args_markers.first_shift_idx
 
     for ip in range(n_markers):
@@ -67,7 +67,7 @@ def grad_driftkinetic_hamiltonian(
             continue
 
         eta_k[:] = markers[ip, 0:3] + markers[ip, first_shift_idx : first_shift_idx + 3]
-        eta_n[:] = markers[ip, first_init_idx : first_init_idx + 3]
+        eta_n[:] = markers[ip, first_pusher_idx : first_pusher_idx + 3]
 
         eta[:] = alpha[:3] * eta_k + (1.0 - alpha[:3]) * eta_n
         eta[:] = mod(eta, 1.0)
diff --git a/src/struphy/pic/pushing/kernels/push_deterministic_diffusion_stage/push_deterministic_diffusion_stage_kernels.py b/src/struphy/pic/pushing/kernels/push_deterministic_diffusion_stage/push_deterministic_diffusion_stage_kernels.py
index a12243267..3b891be49 100644
--- a/src/struphy/pic/pushing/kernels/push_deterministic_diffusion_stage/push_deterministic_diffusion_stage_kernels.py
+++ b/src/struphy/pic/pushing/kernels/push_deterministic_diffusion_stage/push_deterministic_diffusion_stage_kernels.py
@@ -49,7 +49,7 @@ def push_deterministic_diffusion_stage(
     # get marker arguments
     markers = args_markers.markers
     n_markers = args_markers.n_markers
-    first_init_idx = args_markers.first_init_idx
+    first_pusher_idx = args_markers.first_pusher_idx
     first_free_idx = args_markers.first_free_idx
 
     # get number of stages
@@ -121,7 +121,7 @@ def push_deterministic_diffusion_stage(
 
         # update positions for intermediate stages or last stage
         markers[ip, 0:3] = (
-            markers[ip, first_init_idx : first_init_idx + 3]
+            markers[ip, first_pusher_idx : first_pusher_idx + 3]
             + dt * a[stage] * k
             + last * markers[ip, first_free_idx : first_free_idx + 3]
         )
diff --git a/src/struphy/pic/pushing/kernels/push_eta_stage/push_eta_stage_cuda.cu b/src/struphy/pic/pushing/kernels/push_eta_stage/push_eta_stage_cuda.cu
index 978a5f1d5..0f5a1938c 100644
--- a/src/struphy/pic/pushing/kernels/push_eta_stage/push_eta_stage_cuda.cu
+++ b/src/struphy/pic/pushing/kernels/push_eta_stage/push_eta_stage_cuda.cu
@@ -14,7 +14,7 @@
  * @param c Butcher nodes c (n_stages entries; unused, as in the pyccel kernel).
  * @param n_stages Number of stages; a raw device pointer carries no length.
  *
- * Holes (first_init_idx == -1) and boundary particles (last column == -2) are skipped.
+ * Holes (first_pusher_idx == -1) and boundary particles (last column == -2) are skipped.
  * The kinetic boundary conditions are applied per marker after the update.
  */
 extern "C" __global__ void push_eta_stage(double dt, int stage, MarkerArgs args_markers, DomainArgs args_domain,
@@ -24,12 +24,12 @@ extern "C" __global__ void push_eta_stage(double dt, int stage, MarkerArgs args_
     int n_markers = args_markers.n_markers;
     if (ip >= n_markers) return;
 
-    int first_init_idx = args_markers.first_init_idx;
+    int first_pusher_idx = args_markers.first_pusher_idx;
     int first_free_idx = args_markers.first_free_idx;
     double last = stage == n_stages - 1 ? 1. : 0.;
 
     // check if marker is a hole or a boundary particle
-    if (args_markers.markers(ip, first_init_idx) == -1. || args_markers.markers(ip, args_markers.markers.shape[1] - 1) == -2.)
+    if (args_markers.markers(ip, first_pusher_idx) == -1. || args_markers.markers(ip, args_markers.markers.shape[1] - 1) == -2.)
         return;
 
     double e1 = args_markers.markers(ip, 0);
@@ -51,7 +51,7 @@ extern "C" __global__ void push_eta_stage(double dt, int stage, MarkerArgs args_
         // accumulation for last stage
         args_markers.markers(ip, first_free_idx + j) += dt * b[stage] * k[j];
         // update positions for intermediate stages or last stage
-        args_markers.markers(ip, j) = args_markers.markers(ip, first_init_idx + j) + dt * a[stage] * k[j] +
+        args_markers.markers(ip, j) = args_markers.markers(ip, first_pusher_idx + j) + dt * a[stage] * k[j] +
                                       last * args_markers.markers(ip, first_free_idx + j);
     }
 
diff --git a/src/struphy/pic/pushing/kernels/push_eta_stage/push_eta_stage_kernels.py b/src/struphy/pic/pushing/kernels/push_eta_stage/push_eta_stage_kernels.py
index b037d4746..e5f0245f1 100644
--- a/src/struphy/pic/pushing/kernels/push_eta_stage/push_eta_stage_kernels.py
+++ b/src/struphy/pic/pushing/kernels/push_eta_stage/push_eta_stage_kernels.py
@@ -43,7 +43,7 @@ def push_eta_stage(
     # get marker arguments
     markers = args_markers.markers
     n_markers = args_markers.n_markers
-    first_init_idx = args_markers.first_init_idx
+    first_pusher_idx = args_markers.first_pusher_idx
     first_free_idx = args_markers.first_free_idx
 
     if stage == n_stages - 1:
@@ -57,7 +57,7 @@ def push_eta_stage(
     # fmt: on
     for ip in range(n_markers):
         # check if marker is a hole or a boundary particle
-        if markers[ip, first_init_idx] == -1.0 or markers[ip, -1] == -2.0:
+        if markers[ip, first_pusher_idx] == -1.0 or markers[ip, -1] == -2.0:
             continue
 
         e1 = markers[ip, 0]
@@ -85,7 +85,7 @@ def push_eta_stage(
 
         # update positions for intermediate stages or last stage
         markers[ip, 0:3] = (
-            markers[ip, first_init_idx : first_init_idx + 3]
+            markers[ip, first_pusher_idx : first_pusher_idx + 3]
             + dt * a[stage] * k
             + last * markers[ip, first_free_idx : first_free_idx + 3]
         )
diff --git a/src/struphy/pic/pushing/kernels/push_gc_Bstar_discrete_gradient_1st_order/push_gc_Bstar_discrete_gradient_1st_order_kernels.py b/src/struphy/pic/pushing/kernels/push_gc_Bstar_discrete_gradient_1st_order/push_gc_Bstar_discrete_gradient_1st_order_kernels.py
index 32e0b3d64..4f0a4ac73 100644
--- a/src/struphy/pic/pushing/kernels/push_gc_Bstar_discrete_gradient_1st_order/push_gc_Bstar_discrete_gradient_1st_order_kernels.py
+++ b/src/struphy/pic/pushing/kernels/push_gc_Bstar_discrete_gradient_1st_order/push_gc_Bstar_discrete_gradient_1st_order_kernels.py
@@ -82,25 +82,25 @@ def push_gc_Bstar_discrete_gradient_1st_order(
     markers = args_markers.markers
     n_markers = args_markers.n_markers
     mu_idx = args_markers.mu_idx
-    first_init_idx = args_markers.first_init_idx
+    first_pusher_idx = args_markers.first_pusher_idx
     first_shift_idx = args_markers.first_shift_idx
     residual_idx = args_markers.residual_idx
     first_free_idx = args_markers.first_free_idx
 
     for ip in range(n_markers):
         # check if marker is converged or a hole
-        if markers[ip, first_init_idx] == -1.0:
+        if markers[ip, first_pusher_idx] == -1.0:
             continue
 
         eta_k[:] = markers[ip, 0:3] + markers[ip, first_shift_idx : first_shift_idx + 3]
-        eta_n[:] = markers[ip, first_init_idx : first_init_idx + 3]
+        eta_n[:] = markers[ip, first_pusher_idx : first_pusher_idx + 3]
 
         eta_mid[:] = (eta_k + eta_n) / 2.0
         eta_mid[:] = mod(eta_mid, 1.0)
         eta_diff[:] = eta_k - eta_n
 
         v_k = markers[ip, 3]
-        v_n = markers[ip, first_init_idx + 3]
+        v_n = markers[ip, first_pusher_idx + 3]
         v_mid = (v_k + v_n) / 2.0
         v_diff = v_k - v_n
 
diff --git a/src/struphy/pic/pushing/kernels/push_gc_Bstar_discrete_gradient_1st_order_newton/push_gc_Bstar_discrete_gradient_1st_order_newton_kernels.py b/src/struphy/pic/pushing/kernels/push_gc_Bstar_discrete_gradient_1st_order_newton/push_gc_Bstar_discrete_gradient_1st_order_newton_kernels.py
index fb88521fc..8f296d666 100644
--- a/src/struphy/pic/pushing/kernels/push_gc_Bstar_discrete_gradient_1st_order_newton/push_gc_Bstar_discrete_gradient_1st_order_newton_kernels.py
+++ b/src/struphy/pic/pushing/kernels/push_gc_Bstar_discrete_gradient_1st_order_newton/push_gc_Bstar_discrete_gradient_1st_order_newton_kernels.py
@@ -143,23 +143,23 @@ def push_gc_Bstar_discrete_gradient_1st_order_newton(
     markers = args_markers.markers
     n_markers = args_markers.n_markers
     mu_idx = args_markers.mu_idx
-    first_init_idx = args_markers.first_init_idx
+    first_pusher_idx = args_markers.first_pusher_idx
     first_shift_idx = args_markers.first_shift_idx
     residual_idx = args_markers.residual_idx
     first_free_idx = args_markers.first_free_idx
 
     for ip in range(n_markers):
         # check if marker is converged or a hole
-        if markers[ip, first_init_idx] == -1.0:
+        if markers[ip, first_pusher_idx] == -1.0:
             continue
 
         eta_k[:] = markers[ip, 0:3]
-        eta_n[:] = markers[ip, first_init_idx : first_init_idx + 3]
+        eta_n[:] = markers[ip, first_pusher_idx : first_pusher_idx + 3]
         eta_k_shifted[:] = eta_k + markers[ip, first_shift_idx : first_shift_idx + 3]
         eta_diff[:] = eta_k_shifted - eta_n
 
         v_k = markers[ip, 3]
-        v_n = markers[ip, first_init_idx + 3]
+        v_n = markers[ip, first_pusher_idx + 3]
         v_diff = v_k - v_n
 
         mu = markers[ip, mu_idx]
diff --git a/src/struphy/pic/pushing/kernels/push_gc_Bstar_discrete_gradient_2nd_order/push_gc_Bstar_discrete_gradient_2nd_order_kernels.py b/src/struphy/pic/pushing/kernels/push_gc_Bstar_discrete_gradient_2nd_order/push_gc_Bstar_discrete_gradient_2nd_order_kernels.py
index 5fe3aadac..3c6b8d82e 100644
--- a/src/struphy/pic/pushing/kernels/push_gc_Bstar_discrete_gradient_2nd_order/push_gc_Bstar_discrete_gradient_2nd_order_kernels.py
+++ b/src/struphy/pic/pushing/kernels/push_gc_Bstar_discrete_gradient_2nd_order/push_gc_Bstar_discrete_gradient_2nd_order_kernels.py
@@ -108,25 +108,25 @@ def push_gc_Bstar_discrete_gradient_2nd_order(
     markers = args_markers.markers
     n_markers = args_markers.n_markers
     mu_idx = args_markers.mu_idx
-    first_init_idx = args_markers.first_init_idx
+    first_pusher_idx = args_markers.first_pusher_idx
     first_shift_idx = args_markers.first_shift_idx
     residual_idx = args_markers.residual_idx
     first_free_idx = args_markers.first_free_idx
 
     for ip in range(n_markers):
         # check if marker is converged or a hole
-        if markers[ip, first_init_idx] == -1.0:
+        if markers[ip, first_pusher_idx] == -1.0:
             continue
 
         eta_k[:] = markers[ip, 0:3] + markers[ip, first_shift_idx : first_shift_idx + 3]
-        eta_n[:] = markers[ip, first_init_idx : first_init_idx + 3]
+        eta_n[:] = markers[ip, first_pusher_idx : first_pusher_idx + 3]
 
         eta_mid[:] = (eta_k + eta_n) / 2.0
         eta_mid[:] = mod(eta_mid, 1.0)
         eta_diff[:] = eta_k - eta_n
 
         v_k = markers[ip, 3]
-        v_n = markers[ip, first_init_idx + 3]
+        v_n = markers[ip, first_pusher_idx + 3]
         v_mid = (v_k + v_n) / 2.0
         v_diff = v_k - v_n
 
diff --git a/src/struphy/pic/pushing/kernels/push_gc_Bstar_explicit_multistage/push_gc_Bstar_explicit_multistage_kernels.py b/src/struphy/pic/pushing/kernels/push_gc_Bstar_explicit_multistage/push_gc_Bstar_explicit_multistage_kernels.py
index f47a562a9..03ac967b6 100644
--- a/src/struphy/pic/pushing/kernels/push_gc_Bstar_explicit_multistage/push_gc_Bstar_explicit_multistage_kernels.py
+++ b/src/struphy/pic/pushing/kernels/push_gc_Bstar_explicit_multistage/push_gc_Bstar_explicit_multistage_kernels.py
@@ -80,7 +80,7 @@ def push_gc_Bstar_explicit_multistage(
     markers = args_markers.markers
     n_markers = args_markers.n_markers
     mu_idx = args_markers.mu_idx
-    first_init_idx = args_markers.first_init_idx
+    first_pusher_idx = args_markers.first_pusher_idx
     first_free_idx = args_markers.first_free_idx
 
     # get number of stages
@@ -93,7 +93,7 @@ def push_gc_Bstar_explicit_multistage(
 
     for ip in range(n_markers):
         # check if marker is a hole
-        if markers[ip, first_init_idx] == -1.0:
+        if markers[ip, first_pusher_idx] == -1.0:
             continue
 
         # if stage == 0.:
@@ -207,11 +207,13 @@ def push_gc_Bstar_explicit_multistage(
 
         # update positions for intermediate stages or last stage
         markers[ip, 0:3] = (
-            markers[ip, first_init_idx : first_init_idx + 3]
+            markers[ip, first_pusher_idx : first_pusher_idx + 3]
             + dt * a[stage] * k
             + last * markers[ip, first_free_idx : first_free_idx + 3]
         )
-        markers[ip, 3] = markers[ip, first_init_idx + 3] + dt * a[stage] * k_v + last * markers[ip, first_free_idx + 3]
+        markers[ip, 3] = (
+            markers[ip, first_pusher_idx + 3] + dt * a[stage] * k_v + last * markers[ip, first_free_idx + 3]
+        )
 
         # apply kinetic boundary conditions
         pusher_utilities_kernels.apply_kinetic_bc_marker(ip, args_markers, args_domain, False)
diff --git a/src/struphy/pic/pushing/kernels/push_gc_bxEstar_dg_1st_order_newton/push_gc_bxEstar_dg_1st_order_newton_kernels.py b/src/struphy/pic/pushing/kernels/push_gc_bxEstar_dg_1st_order_newton/push_gc_bxEstar_dg_1st_order_newton_kernels.py
index f383d4329..3d9c1a3e2 100644
--- a/src/struphy/pic/pushing/kernels/push_gc_bxEstar_dg_1st_order_newton/push_gc_bxEstar_dg_1st_order_newton_kernels.py
+++ b/src/struphy/pic/pushing/kernels/push_gc_bxEstar_dg_1st_order_newton/push_gc_bxEstar_dg_1st_order_newton_kernels.py
@@ -129,18 +129,18 @@ def push_gc_bxEstar_dg_1st_order_newton(
     markers = args_markers.markers
     n_markers = args_markers.n_markers
     mu_idx = args_markers.mu_idx
-    first_init_idx = args_markers.first_init_idx
+    first_pusher_idx = args_markers.first_pusher_idx
     first_shift_idx = args_markers.first_shift_idx
     residual_idx = args_markers.residual_idx
     first_free_idx = args_markers.first_free_idx
 
     for ip in range(n_markers):
         # check if marker is converged or a hole
-        if markers[ip, first_init_idx] == -1.0:
+        if markers[ip, first_pusher_idx] == -1.0:
             continue
 
         eta_k[:] = markers[ip, 0:3]
-        eta_n[:] = markers[ip, first_init_idx : first_init_idx + 3]
+        eta_n[:] = markers[ip, first_pusher_idx : first_pusher_idx + 3]
         eta_k_shifted[:] = eta_k + markers[ip, first_shift_idx : first_shift_idx + 3]
         eta_diff[:] = eta_k_shifted - eta_n
 
diff --git a/src/struphy/pic/pushing/kernels/push_gc_bxEstar_discrete_gradient_1st_order/push_gc_bxEstar_discrete_gradient_1st_order_kernels.py b/src/struphy/pic/pushing/kernels/push_gc_bxEstar_discrete_gradient_1st_order/push_gc_bxEstar_discrete_gradient_1st_order_kernels.py
index ca3507264..7c26a86dd 100644
--- a/src/struphy/pic/pushing/kernels/push_gc_bxEstar_discrete_gradient_1st_order/push_gc_bxEstar_discrete_gradient_1st_order_kernels.py
+++ b/src/struphy/pic/pushing/kernels/push_gc_bxEstar_discrete_gradient_1st_order/push_gc_bxEstar_discrete_gradient_1st_order_kernels.py
@@ -75,18 +75,18 @@ def push_gc_bxEstar_discrete_gradient_1st_order(
     markers = args_markers.markers
     n_markers = args_markers.n_markers
     mu_idx = args_markers.mu_idx
-    first_init_idx = args_markers.first_init_idx
+    first_pusher_idx = args_markers.first_pusher_idx
     first_shift_idx = args_markers.first_shift_idx
     residual_idx = args_markers.residual_idx
     first_free_idx = args_markers.first_free_idx
 
     for ip in range(n_markers):
         # check if marker is converged or a hole
-        if markers[ip, first_init_idx] == -1.0:
+        if markers[ip, first_pusher_idx] == -1.0:
             continue
 
         eta_k[:] = markers[ip, 0:3] + markers[ip, first_shift_idx : first_shift_idx + 3]
-        eta_n[:] = markers[ip, first_init_idx : first_init_idx + 3]
+        eta_n[:] = markers[ip, first_pusher_idx : first_pusher_idx + 3]
 
         eta_mid[:] = (eta_k + eta_n) / 2.0
         eta_mid[:] = mod(eta_mid, 1.0)
diff --git a/src/struphy/pic/pushing/kernels/push_gc_bxEstar_discrete_gradient_2nd_order/push_gc_bxEstar_discrete_gradient_2nd_order_kernels.py b/src/struphy/pic/pushing/kernels/push_gc_bxEstar_discrete_gradient_2nd_order/push_gc_bxEstar_discrete_gradient_2nd_order_kernels.py
index 619512ae4..bd0457365 100644
--- a/src/struphy/pic/pushing/kernels/push_gc_bxEstar_discrete_gradient_2nd_order/push_gc_bxEstar_discrete_gradient_2nd_order_kernels.py
+++ b/src/struphy/pic/pushing/kernels/push_gc_bxEstar_discrete_gradient_2nd_order/push_gc_bxEstar_discrete_gradient_2nd_order_kernels.py
@@ -91,18 +91,18 @@ def push_gc_bxEstar_discrete_gradient_2nd_order(
     markers = args_markers.markers
     n_markers = args_markers.n_markers
     mu_idx = args_markers.mu_idx
-    first_init_idx = args_markers.first_init_idx
+    first_pusher_idx = args_markers.first_pusher_idx
     first_shift_idx = args_markers.first_shift_idx
     residual_idx = args_markers.residual_idx
     first_free_idx = args_markers.first_free_idx
 
     for ip in range(n_markers):
         # check if marker is converged or a hole
-        if markers[ip, first_init_idx] == -1.0:
+        if markers[ip, first_pusher_idx] == -1.0:
             continue
 
         eta_k[:] = markers[ip, 0:3] + markers[ip, first_shift_idx : first_shift_idx + 3]
-        eta_n[:] = markers[ip, first_init_idx : first_init_idx + 3]
+        eta_n[:] = markers[ip, first_pusher_idx : first_pusher_idx + 3]
 
         eta_mid[:] = (eta_k + eta_n) / 2.0
         eta_mid[:] = mod(eta_mid, 1.0)
diff --git a/src/struphy/pic/pushing/kernels/push_gc_bxEstar_explicit_multistage/push_gc_bxEstar_explicit_multistage_kernels.py b/src/struphy/pic/pushing/kernels/push_gc_bxEstar_explicit_multistage/push_gc_bxEstar_explicit_multistage_kernels.py
index c68033f5d..905bfb8e0 100644
--- a/src/struphy/pic/pushing/kernels/push_gc_bxEstar_explicit_multistage/push_gc_bxEstar_explicit_multistage_kernels.py
+++ b/src/struphy/pic/pushing/kernels/push_gc_bxEstar_explicit_multistage/push_gc_bxEstar_explicit_multistage_kernels.py
@@ -66,7 +66,7 @@ def push_gc_bxEstar_explicit_multistage(
     markers = args_markers.markers
     n_markers = args_markers.n_markers
     mu_idx = args_markers.mu_idx
-    first_init_idx = args_markers.first_init_idx
+    first_pusher_idx = args_markers.first_pusher_idx
     first_free_idx = args_markers.first_free_idx
 
     # get number of stages
@@ -79,7 +79,7 @@ def push_gc_bxEstar_explicit_multistage(
 
     for ip in range(n_markers):
         # check if marker is a hole
-        if markers[ip, first_init_idx] == -1.0:
+        if markers[ip, first_pusher_idx] == -1.0:
             continue
 
         eta1 = markers[ip, 0]
@@ -171,7 +171,7 @@ def push_gc_bxEstar_explicit_multistage(
 
         # update positions for intermediate stages or last stage
         markers[ip, 0:3] = (
-            markers[ip, first_init_idx : first_init_idx + 3]
+            markers[ip, first_pusher_idx : first_pusher_idx + 3]
             + dt * a[stage] * k
             + last * markers[ip, first_free_idx : first_free_idx + 3]
         )
diff --git a/src/struphy/pic/pushing/kernels/push_gc_cc_J2_dg_Hdiv/push_gc_cc_J2_dg_Hdiv_kernels.py b/src/struphy/pic/pushing/kernels/push_gc_cc_J2_dg_Hdiv/push_gc_cc_J2_dg_Hdiv_kernels.py
index 624ab42b5..c10e36e0a 100644
--- a/src/struphy/pic/pushing/kernels/push_gc_cc_J2_dg_Hdiv/push_gc_cc_J2_dg_Hdiv_kernels.py
+++ b/src/struphy/pic/pushing/kernels/push_gc_cc_J2_dg_Hdiv/push_gc_cc_J2_dg_Hdiv_kernels.py
@@ -80,7 +80,7 @@ def push_gc_cc_J2_dg_Hdiv(
     markers = args_markers.markers
     n_markers = args_markers.n_markers
     mu_idx = args_markers.mu_idx
-    first_init_idx = args_markers.first_init_idx
+    first_pusher_idx = args_markers.first_pusher_idx
     first_free_idx = args_markers.first_free_idx
 
     for ip in range(n_markers):
@@ -90,7 +90,7 @@ def push_gc_cc_J2_dg_Hdiv(
 
         # marker positions, mid point
         eta_old[:] = markers[ip, 0:3]
-        eta_mid[:] = (markers[ip, 0:3] + markers[ip, first_init_idx : first_init_idx + 3]) / 2.0
+        eta_mid[:] = (markers[ip, 0:3] + markers[ip, first_pusher_idx : first_pusher_idx + 3]) / 2.0
         eta_mid[:] = mod(eta_mid[:], 1.0)
 
         v = markers[ip, 3]
@@ -204,6 +204,6 @@ def push_gc_cc_J2_dg_Hdiv(
         e /= abs_b_star_para
         e /= det_df
 
-        markers[ip, 0:3] = markers[ip, first_init_idx : first_init_idx + 3] - dt * e
+        markers[ip, 0:3] = markers[ip, first_pusher_idx : first_pusher_idx + 3] - dt * e
         markers[ip, 0:3] *= alpha
         markers[ip, 0:3] += eta_old * (1.0 - alpha)
diff --git a/src/struphy/pic/pushing/kernels/push_gc_cc_J2_dg_init_Hdiv/push_gc_cc_J2_dg_init_Hdiv_kernels.py b/src/struphy/pic/pushing/kernels/push_gc_cc_J2_dg_init_Hdiv/push_gc_cc_J2_dg_init_Hdiv_kernels.py
index 8fcaafd1c..388c069cc 100644
--- a/src/struphy/pic/pushing/kernels/push_gc_cc_J2_dg_init_Hdiv/push_gc_cc_J2_dg_init_Hdiv_kernels.py
+++ b/src/struphy/pic/pushing/kernels/push_gc_cc_J2_dg_init_Hdiv/push_gc_cc_J2_dg_init_Hdiv_kernels.py
@@ -67,12 +67,12 @@ def push_gc_cc_J2_dg_init_Hdiv(
     markers = args_markers.markers
     n_markers = args_markers.n_markers
     mu_idx = args_markers.mu_idx
-    first_init_idx = args_markers.first_init_idx
+    first_pusher_idx = args_markers.first_pusher_idx
     first_free_idx = args_markers.first_free_idx
 
     for ip in range(n_markers):
         # check if marker is a hole
-        if markers[ip, first_init_idx] == -1.0:
+        if markers[ip, first_pusher_idx] == -1.0:
             continue
 
         eta1 = markers[ip, 0]
diff --git a/src/struphy/pic/pushing/kernels/push_gc_cc_J2_stage_H1vec/push_gc_cc_J2_stage_H1vec_kernels.py b/src/struphy/pic/pushing/kernels/push_gc_cc_J2_stage_H1vec/push_gc_cc_J2_stage_H1vec_kernels.py
index a3cb0e4e5..d40c24fc7 100644
--- a/src/struphy/pic/pushing/kernels/push_gc_cc_J2_stage_H1vec/push_gc_cc_J2_stage_H1vec_kernels.py
+++ b/src/struphy/pic/pushing/kernels/push_gc_cc_J2_stage_H1vec/push_gc_cc_J2_stage_H1vec_kernels.py
@@ -85,7 +85,7 @@ def push_gc_cc_J2_stage_H1vec(
     # get marker arguments
     markers = args_markers.markers
     n_markers = args_markers.n_markers
-    first_init_idx = args_markers.first_init_idx
+    first_pusher_idx = args_markers.first_pusher_idx
     first_free_idx = args_markers.first_free_idx
 
     # get number of stages
@@ -98,7 +98,7 @@ def push_gc_cc_J2_stage_H1vec(
 
     for ip in range(n_markers):
         # check if marker is a hole
-        if markers[ip, first_init_idx] == -1.0:
+        if markers[ip, first_pusher_idx] == -1.0:
             continue
 
         eta1 = markers[ip, 0]
@@ -203,7 +203,7 @@ def push_gc_cc_J2_stage_H1vec(
 
         # update positions for intermediate stages or last stage
         markers[ip, 0:3] = (
-            markers[ip, first_init_idx : first_init_idx + 3]
+            markers[ip, first_pusher_idx : first_pusher_idx + 3]
             - dt * a[stage] * e
             + last * markers[ip, first_free_idx : first_free_idx + 3]
         )
diff --git a/src/struphy/pic/pushing/kernels/push_gc_cc_J2_stage_Hdiv/push_gc_cc_J2_stage_Hdiv_kernels.py b/src/struphy/pic/pushing/kernels/push_gc_cc_J2_stage_Hdiv/push_gc_cc_J2_stage_Hdiv_kernels.py
index ab5e5f89d..c2a82a5f5 100644
--- a/src/struphy/pic/pushing/kernels/push_gc_cc_J2_stage_Hdiv/push_gc_cc_J2_stage_Hdiv_kernels.py
+++ b/src/struphy/pic/pushing/kernels/push_gc_cc_J2_stage_Hdiv/push_gc_cc_J2_stage_Hdiv_kernels.py
@@ -81,7 +81,7 @@ def push_gc_cc_J2_stage_Hdiv(
     markers = args_markers.markers
     n_markers = args_markers.n_markers
     mu_idx = args_markers.mu_idx
-    first_init_idx = args_markers.first_init_idx
+    first_pusher_idx = args_markers.first_pusher_idx
     first_free_idx = args_markers.first_free_idx
 
     # get number of stages
@@ -98,7 +98,7 @@ def push_gc_cc_J2_stage_Hdiv(
             continue
 
         # if particle is refilled
-        if markers[ip, first_init_idx] == -1.0:
+        if markers[ip, first_pusher_idx] == -1.0:
             continue
 
         eta1 = markers[ip, 0]
@@ -204,7 +204,7 @@ def push_gc_cc_J2_stage_Hdiv(
 
         # update positions for intermediate stages or last stage
         markers[ip, 0:3] = (
-            markers[ip, first_init_idx : first_init_idx + 3]
+            markers[ip, first_pusher_idx : first_pusher_idx + 3]
             - dt * a[stage] * e
             + last * markers[ip, first_free_idx : first_free_idx + 3]
         )
diff --git a/src/struphy/pic/pushing/kernels/push_pc_eta_stage_H1vec/push_pc_eta_stage_H1vec_kernels.py b/src/struphy/pic/pushing/kernels/push_pc_eta_stage_H1vec/push_pc_eta_stage_H1vec_kernels.py
index 5946bdf4b..e9fe1c0c5 100644
--- a/src/struphy/pic/pushing/kernels/push_pc_eta_stage_H1vec/push_pc_eta_stage_H1vec_kernels.py
+++ b/src/struphy/pic/pushing/kernels/push_pc_eta_stage_H1vec/push_pc_eta_stage_H1vec_kernels.py
@@ -66,7 +66,7 @@ def push_pc_eta_stage_H1vec(
     # get marker arguments
     markers = args_markers.markers
     n_markers = args_markers.n_markers
-    first_init_idx = args_markers.first_init_idx
+    first_pusher_idx = args_markers.first_pusher_idx
     first_free_idx = args_markers.first_free_idx
 
     # get number of stages
@@ -131,7 +131,7 @@ def push_pc_eta_stage_H1vec(
 
         # update markers for the next stage
         markers[ip, 0:3] = (
-            markers[ip, first_init_idx : first_init_idx + 3]
+            markers[ip, first_pusher_idx : first_pusher_idx + 3]
             + dt * k * a[stage]
             + last * markers[ip, first_free_idx : first_free_idx + 3]
         )
diff --git a/src/struphy/pic/pushing/kernels/push_pc_eta_stage_Hcurl/push_pc_eta_stage_Hcurl_kernels.py b/src/struphy/pic/pushing/kernels/push_pc_eta_stage_Hcurl/push_pc_eta_stage_Hcurl_kernels.py
index 33004306d..55e02b6ff 100644
--- a/src/struphy/pic/pushing/kernels/push_pc_eta_stage_Hcurl/push_pc_eta_stage_Hcurl_kernels.py
+++ b/src/struphy/pic/pushing/kernels/push_pc_eta_stage_Hcurl/push_pc_eta_stage_Hcurl_kernels.py
@@ -59,7 +59,7 @@ def push_pc_eta_stage_Hcurl(
     # get marker arguments
     markers = args_markers.markers
     n_markers = args_markers.n_markers
-    first_init_idx = args_markers.first_init_idx
+    first_pusher_idx = args_markers.first_pusher_idx
     first_free_idx = args_markers.first_free_idx
 
     # get number of stages
@@ -127,7 +127,7 @@ def push_pc_eta_stage_Hcurl(
 
         # update markers for the next stage
         markers[ip, 0:3] = (
-            markers[ip, first_init_idx : first_init_idx + 3]
+            markers[ip, first_pusher_idx : first_pusher_idx + 3]
             + dt * k * a[stage]
             + last * markers[ip, first_free_idx : first_free_idx + 3]
         )
diff --git a/src/struphy/pic/pushing/kernels/push_pc_eta_stage_Hdiv/push_pc_eta_stage_Hdiv_kernels.py b/src/struphy/pic/pushing/kernels/push_pc_eta_stage_Hdiv/push_pc_eta_stage_Hdiv_kernels.py
index 6ced6a2f9..0cae6016c 100644
--- a/src/struphy/pic/pushing/kernels/push_pc_eta_stage_Hdiv/push_pc_eta_stage_Hdiv_kernels.py
+++ b/src/struphy/pic/pushing/kernels/push_pc_eta_stage_Hdiv/push_pc_eta_stage_Hdiv_kernels.py
@@ -59,7 +59,7 @@ def push_pc_eta_stage_Hdiv(
     # get marker arguments
     markers = args_markers.markers
     n_markers = args_markers.n_markers
-    first_init_idx = args_markers.first_init_idx
+    first_pusher_idx = args_markers.first_pusher_idx
     first_free_idx = args_markers.first_free_idx
 
     # get number of stages
@@ -128,7 +128,7 @@ def push_pc_eta_stage_Hdiv(
 
         # update markers for the next stage
         markers[ip, 0:3] = (
-            markers[ip, first_init_idx : first_init_idx + 3]
+            markers[ip, first_pusher_idx : first_pusher_idx + 3]
             + dt * k * a[stage]
             + last * markers[ip, first_free_idx : first_free_idx + 3]
         )
diff --git a/src/struphy/pic/pushing/kernels/push_vxb_analytic/push_vxb_analytic_cuda.cu b/src/struphy/pic/pushing/kernels/push_vxb_analytic/push_vxb_analytic_cuda.cu
index bc96bfd77..d14f81ec6 100644
--- a/src/struphy/pic/pushing/kernels/push_vxb_analytic/push_vxb_analytic_cuda.cu
+++ b/src/struphy/pic/pushing/kernels/push_vxb_analytic/push_vxb_analytic_cuda.cu
@@ -28,10 +28,10 @@ extern "C" __global__ void push_vxb_analytic(double dt, int stage, MarkerArgs ar
     int n_markers = args_markers.n_markers;
     if (ip >= n_markers) return;
 
-    int first_init_idx = args_markers.first_init_idx;
+    int first_pusher_idx = args_markers.first_pusher_idx;
 
     // check if marker is a hole
-    if (args_markers.markers(ip, first_init_idx) == -1. || args_markers.markers(ip, args_markers.markers.shape[1] - 1) == -2.)
+    if (args_markers.markers(ip, first_pusher_idx) == -1. || args_markers.markers(ip, args_markers.markers.shape[1] - 1) == -2.)
         return;
 
     double e1 = args_markers.markers(ip, 0);
diff --git a/src/struphy/pic/pushing/kernels/push_vxb_analytic/push_vxb_analytic_kernels.py b/src/struphy/pic/pushing/kernels/push_vxb_analytic/push_vxb_analytic_kernels.py
index 3341e4794..defa85cc3 100644
--- a/src/struphy/pic/pushing/kernels/push_vxb_analytic/push_vxb_analytic_kernels.py
+++ b/src/struphy/pic/pushing/kernels/push_vxb_analytic/push_vxb_analytic_kernels.py
@@ -54,7 +54,7 @@ def push_vxb_analytic(
     # get marker arguments
     markers = args_markers.markers
     n_markers = args_markers.n_markers
-    first_init_idx = args_markers.first_init_idx
+    first_pusher_idx = args_markers.first_pusher_idx
 
     # fmt: off
     #$ omp parallel private (ip, e1, e2, e3, v, dfm, det_df, span1, span2, span3, b_form, b_cart, b_abs, b_norm, vpar, vxb_norm, vperp, b_normxvperp)
@@ -62,7 +62,7 @@ def push_vxb_analytic(
     # fmt: on
     for ip in range(n_markers):
         # check if marker is a hole
-        if markers[ip, first_init_idx] == -1.0 or markers[ip, -1] == -2.0:
+        if markers[ip, first_pusher_idx] == -1.0 or markers[ip, -1] == -2.0:
             continue
 
         e1 = markers[ip, 0]
diff --git a/src/struphy/pic/pushing/kernels/push_vxb_implicit/push_vxb_implicit_cuda.cu b/src/struphy/pic/pushing/kernels/push_vxb_implicit/push_vxb_implicit_cuda.cu
index 1190f3e25..8fbb4a7e2 100644
--- a/src/struphy/pic/pushing/kernels/push_vxb_implicit/push_vxb_implicit_cuda.cu
+++ b/src/struphy/pic/pushing/kernels/push_vxb_implicit/push_vxb_implicit_cuda.cu
@@ -27,10 +27,10 @@ extern "C" __global__ void push_vxb_implicit(double dt, int stage, MarkerArgs ar
     int n_markers = args_markers.n_markers;
     if (ip >= n_markers) return;
 
-    int first_init_idx = args_markers.first_init_idx;
+    int first_pusher_idx = args_markers.first_pusher_idx;
 
     // check if marker is a hole
-    if (args_markers.markers(ip, first_init_idx) == -1.) return;
+    if (args_markers.markers(ip, first_pusher_idx) == -1.) return;
 
     double e1 = args_markers.markers(ip, 0);
     double e2 = args_markers.markers(ip, 1);
diff --git a/src/struphy/pic/pushing/kernels/push_vxb_implicit/push_vxb_implicit_kernels.py b/src/struphy/pic/pushing/kernels/push_vxb_implicit/push_vxb_implicit_kernels.py
index 505959518..2b1c5facf 100644
--- a/src/struphy/pic/pushing/kernels/push_vxb_implicit/push_vxb_implicit_kernels.py
+++ b/src/struphy/pic/pushing/kernels/push_vxb_implicit/push_vxb_implicit_kernels.py
@@ -65,7 +65,7 @@ def push_vxb_implicit(
     # get marker arguments
     markers = args_markers.markers
     n_markers = args_markers.n_markers
-    first_init_idx = args_markers.first_init_idx
+    first_pusher_idx = args_markers.first_pusher_idx
 
     # fmt: off
     #$ omp parallel firstprivate(b_prod) private (ip, e1, e2, e3, v, dfm, det_df, span1, span2, span3, b_form, b_cart, rhs, lhs, lhs_inv, vec, res)
@@ -73,7 +73,7 @@ def push_vxb_implicit(
     # fmt: on
     for ip in range(n_markers):
         # check if marker is a hole
-        if markers[ip, first_init_idx] == -1.0:
+        if markers[ip, first_pusher_idx] == -1.0:
             continue
 
         e1 = markers[ip, 0]
diff --git a/src/struphy/pic/pushing/kernels/unit_b_1form/unit_b_1form_kernels.py b/src/struphy/pic/pushing/kernels/unit_b_1form/unit_b_1form_kernels.py
index 07460ad4e..77cd4aa0a 100644
--- a/src/struphy/pic/pushing/kernels/unit_b_1form/unit_b_1form_kernels.py
+++ b/src/struphy/pic/pushing/kernels/unit_b_1form/unit_b_1form_kernels.py
@@ -39,7 +39,7 @@ def unit_b_1form(
     markers = args_markers.markers
     n_markers = args_markers.n_markers
     mu_idx = args_markers.mu_idx
-    first_init_idx = args_markers.first_init_idx
+    first_pusher_idx = args_markers.first_pusher_idx
     first_shift_idx = args_markers.first_shift_idx
 
     for ip in range(n_markers):
@@ -48,7 +48,7 @@ def unit_b_1form(
             continue
 
         eta_k[:] = markers[ip, 0:3] + markers[ip, first_shift_idx : first_shift_idx + 3]
-        eta_n[:] = markers[ip, first_init_idx : first_init_idx + 3]
+        eta_n[:] = markers[ip, first_pusher_idx : first_pusher_idx + 3]
 
         eta[:] = alpha[:3] * eta_k + (1.0 - alpha[:3]) * eta_n
         eta[:] = mod(eta, 1.0)
diff --git a/src/struphy/pic/pushing/prepare_push_cuda.cu b/src/struphy/pic/pushing/prepare_push_cuda.cu
new file mode 100644
index 000000000..71c02556d
--- /dev/null
+++ b/src/struphy/pic/pushing/prepare_push_cuda.cu
@@ -0,0 +1,28 @@
+/**
+ * Prepare the marker buffer for a push, as the three column-slice assignments at the start of Pusher._push:
+ *
+ *     markers[:, first_pusher_idx:first_shift_idx] = markers[:, :n_init]
+ *     markers[:, first_shift_idx:residual_idx] = 0.0
+ *     markers[:, residual_idx:-2] = 0.0
+ *
+ * One thread per entry of the contiguous column range [first_pusher_idx, n_cols - 2) of every row, so that
+ * neighbouring threads write neighbouring addresses of the row-major buffer and the whole range is written
+ * in a single pass (instead of one strided pass over the buffer per slice).
+ *
+ * @param markers Marker buffer (n_rows x n_cols, row-major), every row including holes.
+ * @param n_rows Number of rows of the buffer.
+ * @param n_cols Number of columns of the buffer.
+ * @param first_pusher_idx First column of the saved initial phase space coordinates.
+ * @param n_init Number of saved coordinates, 3 + vdim (first_shift_idx = first_pusher_idx + n_init).
+ */
+extern "C" __global__ void prepare_push(double* markers, long long n_rows, int n_cols, int first_pusher_idx,
+                                        int n_init) {
+    long long width = n_cols - 2 - first_pusher_idx;
+    long long t = (long long)blockDim.x * blockIdx.x + threadIdx.x;
+    if (t >= n_rows * width) return;
+
+    long long ip = t / width;
+    int col = first_pusher_idx + (int)(t - ip * width);
+    double* row = markers + ip * n_cols;
+    row[col] = col < first_pusher_idx + n_init ? row[col - first_pusher_idx] : 0.0;
+}
diff --git a/src/struphy/pic/pushing/pusher.py b/src/struphy/pic/pushing/pusher.py
index 7f3d64961..762dd3bd9 100644
--- a/src/struphy/pic/pushing/pusher.py
+++ b/src/struphy/pic/pushing/pusher.py
@@ -1,9 +1,11 @@
 "Accelerated particle pushing."
 
 import logging
+from functools import cache
+from pathlib import Path
 
 import cunumpy as xp
-from cunumpy.kernels import Kernel, PyccelKernel
+from cunumpy.kernels import CudaKernel, Kernel, PyccelKernel
 from line_profiler import profile
 from maybempi import MPI
 from scope_profiler import ProfileManager
@@ -16,6 +18,13 @@
 logger = logging.getLogger("struphy")
 
 
+@cache
+def _prepare_push_cuda() -> CudaKernel:
+    """The CUDA kernel that fuses the column-slice assignments at the start of :meth:`Pusher._push`
+    (compiled on first call)."""
+    return CudaKernel.from_file(Path(__file__).with_name("prepare_push_cuda.cu"), "prepare_push")
+
+
 class Pusher:
     r"""
     Class for solving particle ODEs
@@ -208,17 +217,29 @@ def _push(self, dt: float):
         logger.debug(f"{residual_idx =}")
         logger.debug(f"{self.particles.n_cols =}")
 
-        init_slice = slice(first_pusher_idx, first_shift_idx)
-        shift_slice = slice(first_shift_idx, residual_idx)
+        if self._cuda:
+            # the three assignments below in one pass over the buffer (each slice is a strided pass on its own)
+            n_rows, n_cols = markers.shape
+            _prepare_push_cuda()(
+                markers,
+                n_rows,
+                n_cols,
+                first_pusher_idx,
+                3 + vdim,
+                n_threads=n_rows * (n_cols - 2 - first_pusher_idx),
+            )
+        else:
+            init_slice = slice(first_pusher_idx, first_shift_idx)
+            shift_slice = slice(first_shift_idx, residual_idx)
 
-        # save initial phase space coordinates
-        markers[:, init_slice] = markers[:, : 3 + vdim]
+            # save initial phase space coordinates
+            markers[:, init_slice] = markers[:, : 3 + vdim]
 
-        # set boundary shifts to zero
-        markers[:, shift_slice] = 0.0
+            # set boundary shifts to zero
+            markers[:, shift_slice] = 0.0
 
-        # clear buffer columns starting from residual index, dont clear ID (last column) and loc_box
-        markers[:, residual_idx:-2] = 0.0
+            # clear buffer columns starting from residual index, dont clear ID (last column) and loc_box
+            markers[:, residual_idx:-2] = 0.0
 
         rank = self.particles.mpi_rank
         logger.debug(f"rank {rank}: starting {self.kernel} ...")
diff --git a/src/struphy/pic/pushing/pusher_utilities_kernels.cuh b/src/struphy/pic/pushing/pusher_utilities_kernels.cuh
index a40dc406a..48b06ef19 100644
--- a/src/struphy/pic/pushing/pusher_utilities_kernels.cuh
+++ b/src/struphy/pic/pushing/pusher_utilities_kernels.cuh
@@ -38,7 +38,7 @@ __device__ inline void reflect_velocity(long long ip, Array2D<double> markers,
  * @param newton Accumulate periodic shifts if true; overwrite them otherwise.
  *
  * Call after updating the position. Boundary type 0 wraps periodically,
- * type 1 mirrors position and logical velocity and sets first_init_idx to -1,
+ * type 1 mirrors position and logical velocity and sets first_pusher_idx to -1,
  * type 2 clears all columns except the ID to -1, and type 3 is handled on
  * the host. The marker view carries shape and strides; j is a CUDA-only
  * column index replacing Pyccel array slices.
@@ -46,7 +46,7 @@ __device__ inline void reflect_velocity(long long ip, Array2D<double> markers,
 __device__ inline void apply_kinetic_bc_marker(int ip, const MarkerArgs& args_markers,
                                              const DomainArgs& args_domain, bool newton) {
     const long long* bc_type = args_markers.bc_type;
-    int first_init_idx = args_markers.first_init_idx;
+    int first_pusher_idx = args_markers.first_pusher_idx;
     int first_shift_idx = args_markers.first_shift_idx;
 
     // Remove markers before applying periodic or reflecting boundaries.
@@ -92,7 +92,7 @@ __device__ inline void apply_kinetic_bc_marker(int ip, const MarkerArgs& args_ma
         }
     }
     if (n_reflected == 0) return;
-    args_markers.markers(ip, first_init_idx) = -1.;
+    args_markers.markers(ip, first_pusher_idx) = -1.;
     for (int axis = 0; axis < 3; ++axis) {
         if (reflected[axis]) reflect_velocity(ip, args_markers.markers, args_domain, axis);
     }
diff --git a/src/struphy/pic/pushing/pusher_utilities_kernels.py b/src/struphy/pic/pushing/pusher_utilities_kernels.py
index 78b9c7d70..b21c1bfb0 100644
--- a/src/struphy/pic/pushing/pusher_utilities_kernels.py
+++ b/src/struphy/pic/pushing/pusher_utilities_kernels.py
@@ -27,7 +27,7 @@ def apply_kinetic_bc_marker(
     * 0 (periodic): wrap :math:`\eta_i` into :math:`[0, 1)` and store the shift :math:`\pm 1`
       at ``first_shift_idx + axis`` (set for ``newton=False``, accumulated for ``newton=True``).
     * 1 (reflect): mirror :math:`\eta_i` at the boundary, flip the logical velocity component
-      and mark the marker as done for multi-stage/iterative pushers (``first_init_idx = -1``).
+      and mark the marker as done for multi-stage/iterative pushers (``first_pusher_idx = -1``).
     * 2 (remove): turn the marker into a hole (all columns except the ID set to -1).
     * 3 (refill): skipped here, handled by :meth:`~struphy.pic.base.Particles.apply_kinetic_bc`.
 
@@ -48,7 +48,7 @@ def apply_kinetic_bc_marker(
 
     markers = args_markers.markers
     bc_type = args_markers.bc_type
-    first_init_idx = args_markers.first_init_idx
+    first_pusher_idx = args_markers.first_pusher_idx
     first_shift_idx = args_markers.first_shift_idx
 
     # remove
@@ -94,7 +94,7 @@ def apply_kinetic_bc_marker(
     if n_reflected == 0:
         return
 
-    markers[ip, first_init_idx] = -1.0
+    markers[ip, first_pusher_idx] = -1.0
 
     # reflect velocities (Jacobian at the reflected position)
     dfm = empty((3, 3), dtype=float)
diff --git a/src/struphy/propagators/current_coupling_5d_gradb.py b/src/struphy/propagators/current_coupling_5d_gradb.py
index 9fac1ec42..a596074df 100644
--- a/src/struphy/propagators/current_coupling_5d_gradb.py
+++ b/src/struphy/propagators/current_coupling_5d_gradb.py
@@ -459,14 +459,14 @@ def __call__(self, dt):
         holes = particles.holes
         args_markers = particles.args_markers
         markers = args_markers.markers
-        first_init_idx = args_markers.first_init_idx
+        first_pusher_idx = args_markers.first_pusher_idx
         first_free_idx = args_markers.first_free_idx
 
         # clear buffer
-        markers[:, first_init_idx:-2] = 0.0
+        markers[:, first_pusher_idx:-2] = 0.0
 
         # save old marker positions
-        markers[:, first_init_idx : first_init_idx + 3] = markers[:, :3]
+        markers[:, first_pusher_idx : first_pusher_idx + 3] = markers[:, :3]
 
         # sum up total magnetic field b_full1 = b_eq + b_tilde (in-place)
         b_full = self._b2.copy(out=self._b_full)
@@ -530,7 +530,7 @@ def __call__(self, dt):
             diffs = self.update_feec_variables(u=u_new)
 
             # clear the buffer
-            markers[:, first_init_idx:-2] = 0.0
+            markers[:, first_pusher_idx:-2] = 0.0
 
             # update_weights
             if self.variables.energetic_ions.species.weights_params.control_variate:
@@ -661,7 +661,7 @@ def __call__(self, dt):
                 sum_u_diff_loc = xp.sum((u_diff.toarray() ** 2))
 
                 sum_H_diff_loc = xp.sum(
-                    (markers[~holes, :3] - markers[~holes, first_init_idx : first_init_idx + 3]) ** 2,
+                    (markers[~holes, :3] - markers[~holes, first_pusher_idx : first_pusher_idx + 3]) ** 2,
                 )
 
                 buffer_array = xp.array([sum_u_diff_loc])
@@ -853,7 +853,7 @@ def __call__(self, dt):
             diffs = self.update_feec_variables(u=u_new)
 
             # clear the buffer
-            markers[:, first_init_idx:-2] = 0.0
+            markers[:, first_pusher_idx:-2] = 0.0
 
             # update_weights
             if self.variables.energetic_ions.species.weights_params.control_variate: