{'fraction_non_alphanumeric': 0.040702314445331206, 'fraction_numerical': 0.014139649746269332, 'mean_word_length': 4.4816343375980185, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A coupled guiding center-Boris particle pusher for magnetized plasmas in compact-object magnetospheres

We present a novel numerical scheme for simulating the motion of relativistic charged particles in magnetospheres of compact objects, typically filled with highly magnetized collisionless plasmas. The new algorithm is based on a dynamic switch between the full system of equations of motion and a guiding center approximation. The switch between the two formulations is based on the magnetization of the plasma particles, such that the dynamics are accurately captured by the guiding center motion even when the gyro-frequency is under-resolved by the time step. For particles with a large gyro-radius, due to acceleration in, e.g., reconnecting current sheets, the algorithm adaptively switches to solve the full equations of motion instead. The new scheme is directly compatible with standard Particle-in-Cell codes, and is readily applicable in curved spacetimes via a dedicated covariant formulation. We test the performance of the coupled algorithm by evolving charged particles in electromagnetic configurations of reconnecting current sheets in magnetized plasma, obtained from special- and general-relativistic Particle-in-Cell simulations. The new coupled pusher is capable of producing highly accurate particle trajectories even when the time step is many orders of magnitude larger than the gyro-period, substantially reducing the restrictions of the temporal resolution.


Introduction
Magnetospheres of compact objects, such as black holes (BHs) and neutron stars (NSs), are typically filled with highly magnetized collisionless plasma. The collective plasma dynamics is best studied with a kinetic description governed by the Vlasov-Maxwell equations. The Particle-in-Cell (PiC) algorithm is the most widely employed method to solve the kinetic equations, owing its success to its simplicity, reliability, and remarkable performance on parallel architectures. Special-relativistic (SR) PiC simulations are used to study the dynamics of magnetized relativistic plasma, and have illuminated the physics of fundamental kinetic processes such as particle acceleration in shocks (Spitkovsky 2008;Sironi & Spitkovsky 2009a,b), magnetic reconnection (Guo et al. 2014;Sironi & Spitkovsky 2014;Werner & Uzdensky 2017;Werner et al. 2018;Hakobyan et al. 2019), and turbulence (Zhdankin et al. 2018;Comisso & Sironi 2018Zhdankin et al. 2019). More recently, PiC simulations have been applied to model the plasma magnetospheres of compact objects, resulting in first-principles models of coherent radio emission from non-stationary reconnection ) and pair discharges , γ-ray emission from reconnecting current sheets (Cerutti et al. 2016;Philippov & Spitkovsky 2018;Kalapotharakos et al. 2018), and X-ray emission in magnetar magnetospheres (Chen & Beloborodov 2017). General-relativistic (GR) PiC simulations have investigated the dynamics of plasmas in curved spacetime, showing that dragging of inertial frames is an essential ingredient for pair creation at the pulsar's polar cap (Philippov et al. 2015), and unveiling the dynamics of pair discharges and jet launching in BH magnetospheres (Levinson & Cerutti 2018;Chen & Yuan 2019;Parfrey et al. 2019;Crinquand et al. 2020).
Collisionless plasma magnetospheres of compact objects typically display a very large separation between the global spatial and temporal scales (e.g., the size of the object and its rotational frequency) and the typical dynamic scales of the particle motion (e.g., the gyro-radius ρ C and the gyro-frequency Ω C ). This imposes constraints on PiC methods, reducing the parameter range that can be realistically probed with simulations. Furthermore, in compact-object magnetospheres where gravity may play an important role, it is imperative to include the effects of curved spacetime, further complicating the solution algorithms. For highly magnetized magnetospheres the scale separation can become extreme, with gyro-radii ∼ 10 10 times smaller than the system size. This difference can be further enhanced by fast (on time scales of order 10 −14 s near surfaces of NSs) synchrotron losses, pushing the particles to the zeroth Landau level. To mitigate the problem, a rescaling of the system's physical parameters (e.g., to lower magnetic field strengths) is typically employed, artificially decreasing the scale separation.
A compromise is provided by applying a guiding center approximation (GCA; see e.g., Northrop 1961Northrop , 1963 reducing the particle motion to that of its guiding center. In the GCA equations, the particle dynamics is averaged over its gyro-motion, and therefore no constraints related to particle gyration are imposed on the time step. This approach has been successfully applied to study test particle motion in the solar corona (Rosdahl & Galsgaard 2009;Gordovskyy et al. 2010;Gordovskyy & Browning 2011a,b;Gordovskyy et al. 2014;Zhou et al. 2015;Pinto et al. 2016;Zhou et al. 2016;Ripperda et al. 2017a,b;Threlfall et al. 2017), in the Earth's magnetotail (Leroy et al. 2019), and has been used in PiC simulations of pair discharges . However, the GCA is based on the assumption that the particle gyro-radius is much smaller than the typical length scale of electromagnetic fluctuations. This implies that particles whose gyro-radius increases significantly, e.g., due to acceleration or a local decrease of the magnetic field strength, may not be accurately described by the GCA equations. This scenario is expected in compact-object magnetospheres, where particles mainly travel along magnetic field lines with negligible gyro-radius, but may get accelerated and gain significant transverse momentum in reconnecting current sheets. The GCA alone is therefore inadequate to accurately model the acceleration of particles in such situations.
In this paper, we present a reliable particle pusher that both accurately captures particle dynamics when its gyro-frequency is significantly under-resolved by the time step, and when its gyro-radius becomes large due to acceleration in, e.g., reconnecting current sheets. We develop a novel algorithm that couples a leapfrog GCA solver with a standard Boris (Boris 1970) scheme, and switches between the two approaches according to criteria based on the magnetization of the plasma particles. The coupled method is also applicable in curved spacetimes, via a dedicated covariant formulation. The paper is organized as follows: in Section 2 we review the relativistic equations of motion for charged particles in electromagnetic fields, and the GCA equations for the motion of the guiding center. We discuss the numerical methods employed for the solution of the equations of motion, presenting two leapfrog-like numerical schemes for the solution of the SR-GCA equations, and a method for the solution of the GR-GCA system. In Section 3 we present the new coupled Boris-GCA pusher, providing the necessary steps for the construction and implementation of the algorithm. We also present a GR version of the scheme which is applicable to the motion of charged particles in strong gravitational fields. In Section 4, we test the new coupled pusher by simulating the motion of charged particles in an isolated reconnecting current sheet in flat spacetime and a current sheet formed in a BH magnetosphere. Finally, in Section 5 we summarize our results and present our main conclusions.

Equations of motion and numerical methods for charged particles in electromagnetic fields
In this Section we briefly review the theoretical framework describing the motion of charged particles in electromagnetic fields. We consider both the case of motion in flat spacetime (i.e. in the absence of gravity) and in curved spacetimes. In each case we first present the relevant equations of motion, and then describe the numerical algorithms we employ to solve them.

Full system of equations in flat spacetime
The evolution equations for the particle position, x, and 4-velocity, u, read where Γ = 1 + u 2 /c 2 is the Lorentz factor, c is the speed of light, and q and m are the particle charge and mass. Equations (1)-(2) are typically approximated with a second-order time discretization, where x and u are staggered in time. The solution of equation (2) is usually carried out with a structure-preserving algorithm such as the Boris approach (Boris 1970), which possesses desirable qualities (e.g., the conservation of phase-space volume). Several other strategies have been presented over the years, to address specific numerical requirements (see  and references therein for a comprehensive comparison of particle pushers). All approaches are however constrained on the time step, which must be typically kept to a value comparable to the gyro-period t C = 2πmcΓ/(|q|B). This condition may become extremely demanding for strong magnetic fields, as is often the case for compact-object magnetospheres. Hereafter, we will refer to equations (1)-(2) as the "SREOM" system, for clarity.

Guiding Center Approximation in flat spacetime
In the guiding center formalism, the guiding center position R is evolved in place of x (see e.g., Northrop 1961Northrop , 1963. The particle 4-velocity u is split into the parallel and perpendicular (to the magnetic field) components u and u ⊥ . The perpendicular 4-velocity can be further divided into a drift component u ⊥d and a gyration component u ⊥g , describing the guiding center drift across magnetic field lines and the particle gyro-motion (at the cyclotron frequency Ω C = |q|B/(mcΓ)) in the plane perpendicular to B, respectively. The projection of u along the parallel and perpendicular directions is given by where b = B/B. While the parallel component u is directly evolved in time, the gyration component u ⊥g is instead determined by the evolution of the magnetic moment µ = mu 2 ⊥g /(2Bκ), where κ = 1/ 1 − v 2 E /c 2 is the Lorentz factor associated with a frame of reference moving with velocity with w E = cE × B/(E 2 + B 2 ). This velocity is chosen such that in the frame moving with v E , E is parallel to B1.
In the nonrelativistic limit E B we retrieve v E = cE × B/B 2 , the usual E × B velocity. The magnetic moment µ is an adiabatic invariant of the motion (to first order) in the GCA equations. The validity of the GCA relies (among others) on the assumption that the particle gyro-radius ρ C = mcu ⊥g /(|q|B) is much smaller than the length scales characterizing spatial variations of the electromagnetic fields. Therefore, numerical methods based on the GCA require careful checks on the validity of this assumption, which may be violated when particles experience strong acceleration (e.g., in reconnection regions, around X-points, where E > B) in the direction perpendicular to B, resulting in non-negligible gyro-motion.
The GCA equations are derived by means of an expansion of the SREOM system in powers of ρ C /L (where L is the characteristic length scale of variation of electromagnetic fields), resulting in drift terms (perpendicular to magnetic field lines) proportional to spatial and temporal variations of E and B (Vandervoort 1960;Northrop 1961Northrop , 1963. In the expansion, the v E drift represents the only zeroth-order term; if E B and temporal derivatives of the electromagnetic fields are neglected (under the assumption that particle dynamics takes place on much faster time scales than the electromagnetic field evolution), two additional first-order terms arise, namely the curvature and ∇B drifts. The curvature drift is given by and it can be shown that the relative magnitude of where ρ C = mcΓ/(|q|B) is the typical gyro-radius of particles traveling with Lorentz factor Γ. It can be seen that, as the particle magnetization (parametrized by L/ρ C ) increases, the curvature drift becomes progressively less 1When E · B = 0, this expression reduces to a familiar drift velocity with , which corresponds to vanishing electric (if E < B) or magnetic (if E > B) field in the frame moving with v E . important than the E × B drift, which instead does not scale with ρ C . Similarly, in the evolution equation for u the expansion introduces zeroth-order terms (given by parallel electric fields) and first-order curvature acceleration terms. Again, the relative magnitude of these terms scales such that, for large magnetizations, the effects of curvature are negligible.
The ∇B ("mirror") drift velocity v ∇B is proportional to µ, and it can be shown that the relative magnitude of the mirror and curvature drift velocities scales as v ∇B /v c ∼ v 2 ⊥g /v 2 , where v ⊥g = u ⊥g /Γ and v = u /Γ. First-order mirror force terms, which enter in the evolution equation for u , scale similarly with respect to curvature terms as v 2 ⊥g /(v E v ). In practice, for moderately relativistic particles with v ∼ c, mirror effects can be neglected as long as v ⊥g c (which implies small µ). We verify that in all our experiments this condition is respected at all times; this scenario is representative of physical situations where gyration has negligible importance in the overall particle dynamics (e.g., in the case of compact-object magnetospheres). For these reasons, we choose to discard mirror drift and acceleration terms in the following equations. Note that, in principle, the condition for which mirror effects are negligible should be rigorously checked in numerical methods based on the GCA equations, to avoid large inaccuracies in the results.
By these arguments, in this work we choose to retain v E and v c as the only relevant drift terms (as well as curvature acceleration terms), bearing in mind that the latter may not play a significant role for scenarios in which the plasma magnetization is sufficiently high. Hence, we solve the equations of motion for this "reduced" GCA, where E = E · b. The Lorentz factor here is calculated by averaging energy fluctuations over a gyro-period, We will hereafter refer to equations (8)-(10) as the "SRGCA" system for clarity. Note that neglecting mirror force terms in equation (9) may in principle introduce errors in the particle energy: in the case where E = 0 and curvature terms are neglected, the particle energy as expressed by equation (11) should remain constant, but since there is no mechanism of conversion of perpendicular velocity into parallel velocity, solving equation (10) for u ⊥g will introduce a deviation from exact energy conservation. This effect is in practice negligible as long as mirror forces (determined by µ) are negligible, as discussed above. An alternative approach consists of avoiding the solution of equation (10) and keeping u ⊥g as a constant. In this case, energy is conserved exactly, at the expense of numerically altering µ. In all our experiments this deviation from exact energy conservation does not significantly influence the results, since mirror forces remain negligible at all times. We therefore rely on the former approach, including the solution of equation (10). For compatibility with standard PiC codes, we wish to solve the system above with a leapfrog-like, second-order scheme. Here we present such a scheme consisting of a semi-implicit discretization of the equations of motion, where the position R is centered on integer time steps, while u is centered on half steps. We analyze two subcases, namely the case where curvature effects are negligible and the case that includes curvature. In compact-object magnetospheres, µ is typically quickly brought to zero via synchrotron radiation. Therefore, in the algorithms presented below µ = 0 can be assumed for specific applications.

Zero curvature case
In cases where it is possible to neglect all gradient terms in the SRGCA system (e.g., for highly magnetized plasmas where v c and curvature acceleration terms become negligible), the equations of motion take the simple form These equations can be solved with a leapfrog-like approach which relies on the second-order discrete equations where we use a short-hand notation for quantities that depend on R only, e.g., b n = b(R n ), v n E = v E (R n ), etc. The Lorentz factor is a function of both R and u , and µ is a conserved quantity across a time step. Equations (14) and (15) can be solved, at each time step, in a semi-implicit fashion. While equation (15) can be solved explicitly, equation (14) is nonlinearly implicit in the new position R n+1 . At each time step, starting from R n and u n−1/2 , we therefore apply the following solution procedure: 1. Solve the 4-velocity equation (15) explicitly, 2. Solve equation (14) for R n+1 iteratively, e.g., with a fixed-point method. An appropriate initial guess is given by the previous value R n . At the generic k-th iteration, the nonlinear residual to be evaluated reads In our experiments, 3-4 iterations are sufficient to reach convergence within an absolute iteration error of 10 −7 .
This approach has been successfully applied in multi-dimensional simulations of pair discharges in pulsar polar caps , in a scenario where v E acts as the predominant drift mechanism in the motion of particles.

Nonzero curvature case
If the magnetic field curvature cannot be neglected, it is necessary to include the curvature terms in the SRGCA system. Applying a second-order discretization yields where the curvature drift is a function of both R and u , While the system (19)- (20) is now manifestly more complex than in the zero curvature case, the solution procedure remains analogous and can be carried out according to the following steps: 1. Solve the 4-velocity equation (20) explicitly, via The Lorentz factor Γ(R n , u n+1/2 ) needed in equation (22) can be calculated as the only positive (larger than or equal to 1) solution of a second-order polynomial k 1 Γ 2 + k 2 Γ + k 3 = 0 with coefficients 2. Solve equation (19) for R n+1 iteratively. At the generic k-th iteration, the nonlinear residual to be evaluated reads Again, 3-4 iterations usually suffice to reach convergence within an absolute iteration error of 10 −7 .

Full system of equations in curved spacetimes
For compact-object environments where strong gravitational fields are present, spacetime curvature effects may become important, influencing the motion of particles. A general-relativistic formulation of the equations of motion for a charged particle in the presence of gravitational and electromagnetic fields reads where the 3+1 split of the ADM formalism (Arnowitt et al. 1959) has been employed to define the evolution of the particle position x i and 4-velocity u i with respect to coordinate time t (Roman indices run from 1 to 3). The right-hand side terms correspond to the particle motion and geodesic terms, and to the covariant Lorentz force, Here, α, β i , and γ i j represent the metric functions (lapse, shift, and spatial 3-metric) of a Riemannian manifold g µν (Greek indices run from 0 to 3), with associated Levi-Civita pseudotensor e i jk = √ γ i jk . In this section we employ geometrized units such that c = G = 1. The Lorentz factor in this formalism is defined as Γ = 1 + γ i j u i u j , and the electric and magnetic fields are constructed by projecting the electromagnetic field tensor F µν , For clarity, we will refer to equations (28)-(29) as the "GREOM" system from now on. The highly nonlinear character of the GREOM system does not allow for the application of explicit, leapfroglike discretization schemes. However, the motion along geodesics combined with the Lorentz force still demands for the accurate preservation of phase-space trajectories. For this reason, an appropriate choice for the solution of the system above is the implicit midpoint rule (IMR), which is a symplectic method (see e.g., Hairer et al. 2006). Alternative numerical schemes have been proposed, e.g., based on a Hamiltonian formulation that preserves energy exactly during the simulation Bacchini et al. 2019b). However, the Lorentz force and the geodesic motion may impose dynamics on vastly different time scales (e.g., for particles whose gyro-motion is much faster than the motion along geodesics), resulting in the GREOM system becoming stiff, which complicates the solution procedure.
An alternative approach consists of Strang splitting the equations above, by discretizing the geodesic motion with a symplectic scheme (e.g., IMR), while applying the electromagnetic push of the Lorentz force in a local inertial reference frame (in which spacetime curvature can be ignored for instantaneous pointwise interactions). The projection to and from the local frame is obtained via a set of orthonormal tetrads. In this frame, the Lorentz force in the 4-velocity update simply takes the standard special-relativistic form (equation (2)), and can be solved explicitly, e.g., with a Boris algorithm (see Parfrey et al. 2019;Crinquand et al. 2020; the method will be presented in detail in a future publication (Parfrey et al. 2020, in prep.). The overall scheme retains second-order accuracy. The advance of the particle position x = (x 1 , x 2 , x 3 ) and 4-velocity u = (u 1 , u 2 , u 3 ) from time level n to n + 1 (note that the variables are non-staggered) can be summarized in the following substeps: 1. First half-push: the 4-velocity u n is first projected to a local inertial frame, evaluating the metric functions at x n . Then, an updated 4-velocity u is obtained by applying a Boris push on u n over half a time step. The electromagnetic fields needed for the update are calculated by projecting to the same local frame. After the update, u is projected back to the global coordinate basis.
2. Position update and geodesic motion: the particle position and 4-velocity are evolved in the global coordinates via the 3+1 split geodesic equations, updating u with the contribution from geometric sources. The equations of motion are discretized with an IMR scheme, wherex = (x n+1 + x n )/2 andū = (u + u )/2. The system above is implicit in x n+1 and u , and therefore requires a nonlinear iteration, which usually converges quickly provided that the metric functions do not introduce any stiffness. In our experiments, 3-4 fixed-point iterations with x n and u as initial guesses are generally sufficient to reach convergence within an absolute iteration error of 10 −7 .
3. Second half-push: the 4-velocity u is projected to a local inertial frame, and a second Boris push is carried out over half a time step, bringing u to u n+1 (with metric and fields calculated at x n+1 ).
We note that the half pushes (steps 1. and 3. above) are applied in a local inertial frame only for the sake of simplicity; in fact, a covariant Boris-like push in the global coordinate basis (without the need to project to a local frame) can still be applied explicitly, although it involves more evaluations of the metric functions. This option could be considered for experiments in spacetimes where calculating a suitable tetrad basis is non-straightforward, or too expensive.

Guiding Center Approximation in curved spacetime
When spacetime curvature is considered, a formal treatment of the gyrokinetic approximation leading to the GCA equations becomes much more complicated than the flat-spacetime case presented in Section 2.2 (see e.g., Beklemishev & Tessarotto 1999, 2004Cremaschini et al. 2008). Here, we present a substantially simplified version of the covariant GCA, constructed such that the resulting equations can be solved with an approach similar to that discussed for the GREOM system in Section 2.3. In this way we ensure that the two systems of equations can be employed together in a coupled numerical scheme. We consider the following equations of motion, where R i is the guiding center position, and In these equations, u i has been split into the usual parallel, drift, and gyration components, is the unit vector in the direction parallel to the magnetic field. Here, we again assume the D × B drift as the predominant cross-field motion, and neglect all other drift terms. We therefore take the drift velocity where w E,i = e i jk D j B k /(D l D l + B l B l ), for a drifting reference frame with associated Lorentz factor κ = 1/ 1 − γ i j v E,i v E, j . With these definitions, equation (36) describes the motion of the guiding center along the parallel and drift directions, and equation (37) describes the change in magnitude of the parallel 4-velocity component. The time-variation of the parallel 4-velocity is given by the projection of both the geodesic (apparent) force term G i and the Lorentz force L i . In flat spacetime, the geodesic term vanishes and the SR limit of the SRGCA system is retrieved (see Section 2.2). Finally, equation (38) expresses the conservation of the magnetic moment µ = mγ i j u ⊥g,i u ⊥g, j /(2Bκ). This directly provides u ⊥g ≡ γ i j u ⊥g,i u ⊥g, j , analogously to the special-relativistic case. The Lorentz factor can be calculated from the GCA variables as We will refer to equations (36)-(38) as the "GRGCA" system from here on. The choice of only including the D × B drift may not be applicable for scenarios where, for example, accurate calculation of particle energy losses to curvature radiation is important. However, if radiation losses are negligible, curvature drifts as well as other drift mechanisms will be of second-order importance compared to D × B drifts. Our tests show that our simplified approach already provides satisfactory results for particle motion in the upstream of current sheets in a BH magnetosphere (see Section 4). The numerical solution of the GRGCA system can be carried out with a Strang split approach similar to the scheme presented in the previous Section. Once again, we apply the contribution of the Lorentz force to the 4-velocity update in a local inertial reference frame, via a set of orthonormal tetrads. The position update and the geodesic push are computed iteratively with an IMR scheme. The resulting algorithm for the update of R i and u , applied between time levels n and n + 1, is composed of steps below. For simplicity, in the following we will employ the usual notation for 3-vectors R = (R 1 , R 2 , R 3 ), u = (u 1 , u 2 , u 3 ), etc.
1. First half-push: the parallel 4-velocity u n is projected to a local inertial frame, evaluating the metric functions at R n . In this frame, we update u n to u over half a time step, evaluating the acceleration determined by the parallel electric field D n = D n · b n . Here, we again rely on a short-hand notation for quantities that depend on R only, e.g., D n = D(R n ), etc. The updated parallel 4-velocity u = u b n is then projected back to the global coordinate basis.
2. Position update and geodesic motion: the particle position and 4-velocity are evolved in the global coordinates via the 3+1 split geodesic equations. The equations of motion are discretized with an IMR scheme, where the mid-point iteration variables areR = (R n+1 + R n )/2 andū = (u + u )/2. The system above is implicit in R n+1 and u , and can be solved iteratively. At each iteration step, in order to evaluate the G term, we reconstruct the full velocity vector asū =ū b(R) + v E (R)Γ(R,ū ) +ū ⊥g . The gyration part of the 4-velocity is updated in magnitude,ū ⊥g = 2µκ(R)B(R)/m, but kept fixed in direction, i.e.ū ⊥g =ū ⊥g b n ⊥g where b n ⊥g = u n ⊥g /u n ⊥g . TheṼ term and the Lorentz factor are given by equations (39) and (41) respectively, as functions of R and u . Similarly to the iterative scheme employed for the GREOM system, 3-4 fixed-point iterations with R n and u as initial guesses suffice to reach convergence within an absolute iteration error of 10 −7 , in all our experiments.
3. Second half-push: the updated parallel 4-velocity u is again projected to a local inertial frame, evaluating the metric functions at R n+1 . The updated parallel 4-velocity is calculated in this frame by adding the contribution of the parallel electric field D n+1 = D n+1 · b n+1 . Finally, u n+1 = u n+1 b n+1 is projected back to the global coordinate basis.
We note again that, similarly to the numerical solution of the GREOM system (see previous Section 2.3), the half-push steps 1. and 3. could in principle be carried out in the global coordinates, without the need to project to a local inertial frame. Here, we will employ the same projection procedure in the spirit of keeping the two schemes compatible for a simultaneous use in a coupled algorithm.

Coupling the GCA and the full system of equations
In this work we propose a new coupling of particle pushers that mitigates the issues related to small-scale gyration, by averaging over the gyro-motion and therefore loosening the constraints on the time time step imposed by high magnetizations. For both special-and general-relativistic simulations, the strategy consists of a dynamic switch between the SREOM/GREOM and the SRGCA/GRGCA systems of equations, based on local parameters measured at the particle position. At each time step we determine which set of equations is solved according to the following criteria: • If the particle gyro-radius is found to be smaller than the local cell size ∆ (or a chosen fraction of it), then the SRGCA or GRGCA system is employed for the current time-advancing step. An additional check is adopted to ensure that the value of E is smaller than (a fraction of) B locally, such that strong particle acceleration can be detected prior to a violation of the GCA assumptions.
• If the particle gyro-radius is larger than the local cell size (or if E > B), then the SREOM or GREOM system is adopted.
In numerical simulations where the fields quantities are spatially resolved on the computational grid, the minimum length over which these fields (such as B) can vary corresponds to the grid spacing ∼ ∆ . Therefore, by switching to the SREOM/GREOM system when ρ C ∆ , we also automatically filter out the case where (B/∇B) < ρ C , i.e. the particle gyro-radius exceeds the length scale of variation of B, which constitutes another violation of the GCA assumptions.
The coupled strategy has the advantage that particle gyro-motion on subgrid scales is automatically averaged over in the SRGCA/GRGCA system, since electromagnetic fluctuations at such scales are not captured. In this way, no constraints related to accuracy loss due to fast gyration are imposed on the time step. At the same time, the motion of particles with non-negligible (larger than the cell spacing) gyro-radius is accurately captured by solving the SREOM/GREOM system. The switch between the two sets of equations, applied dynamically, allows for capturing particles transitioning to different physical regimes: particles with negligible gyro-motion can be accurately tracked as they gain energy, until the gyro-radius grows above the cell spacing; while energetic particles with large gyro-radii can be captured as they experience energy losses, transitioning to a regime of negligible gyro-motion.
We note that an effective cell size ∆ to be compared with ρ C can be reasonably defined in relation to the local cell volume, e.g., ∆ = (∆V) 1/n for n-dimensional simulations. In this way, it is possible to account for the general case of non-uniform grid spacing. This criterion remains applicable in general-relativistic simulations, where the cell volume is formally defined via the proper integral evaluated over the cell's spatial extent. In the more common case of special-relativistic simulations with uniform grid spacing, one can simply take ∆ = ∆x.

Special-relativistic coupled scheme
For simulations in flat spacetimes, we apply a Boris push to solve the SREOM system, and the leapfrog pushers described in Sections 2.2.1 and 2.2.2 to solve the SRGCA system. The position and velocity update between consecutive time levels is therefore carried out as follows: 1. Check for the validity of the GCA assumptions at the current particle location; if the SRGCA system can be employed, switch to the GCA variables R, u , µ. 3. After the update, store the new position and 4-velocity. If a GCA push was applied, retrieve x and u from the updated GCA variables.
A key aspect of this coupled approach lies in making the SRGCA system and the full SREOM system compatible with each other, that is, a back-and-forth switch between the two formulations must not introduce inconsistencies. Since the GCA set of equations represents an averaging of the full equations of motion, there is an inevitable loss of information when switching to the former. Specifically, the three components of the particle velocity in threedimensional space are reduced to the components parallel and perpendicular to the magnetic field, and information on the gyro-motion is lost. In our experiments, when switching from x and u to the GCA variables at the beginning of a time step, we first assume R ≈ x (which is a valid approximation in the limit of small gyro-radius), and then decompose the velocity as u n−1/2 = u n−1/2 b n−1/2 + v n−1/2 E + v c (R n−1/2 , u n−1/2 ) Γ(R n−1/2 , u n−1/2 ) + u n−1/2 ⊥g , where u = u · b, and b, v E , and v c are calculated with the electromagnetic fields evaluated at R n−1/2 = (R n + R n−1 )/2. This directly provides u and µ for the integration of the SRGCA system. Our assumption therefore consists of reducing the drift velocity to v E and v c only (which will be proven satisfactory in the following Sections). Similarly, at the end of a GCA integration step, we switch back to the full set of variables by taking the updated particle position as x ≈ R. Further assumptions have to be made to reconstruct the full particle velocity from the GCA variables. The integration of the SRGCA system yields only the magnitude of the new parallel and gyration 4-velocity, u and u ⊥g ; the approximation introduced by the GCA becomes manifest as a lack of information on the direction of the updated 4-velocity. Specifically, while the new direction parallel to the magnetic field is given by R, the direction of the updated u ⊥g is unknown. We therefore choose to store b n−1/2 ⊥g = u n−1/2 ⊥g /u n−1/2 ⊥g at the previous time step, and reconstruct the full 4-velocity (assuming R n+1/2 = (R n+1 + R n )/2) as hence assuming that the direction of u ⊥g (and therefore the gyration phase) is unchanged. In principle, it is possible to separately evolve the gyration phase by adding an extra equation to the SRGCA system. Here, our choice of keeping the direction of u ⊥g is justified by the fact that, if the GCA is applicable at the current integration step, the direction of u ⊥g does not play a significant role in the particle dynamics.

General-relativistic coupled scheme
The coupling of the GREOM and GRGCA systems (see Sections 2.3 and 2.4) can be carried out in a fashion analogous to the special-relativistic case. Notable differences lie in the fact that position and 4-velocity are now non-staggered, and that a Strang split is employed to separate the Lorentz force contribution from that of the geodesic terms. One integration step is carried out as follows: 1. Check for the validity of the GCA assumptions at the current particle location; if the GRGCA system can be employed, switch to the GCA variables R i ≈ x i , u , µ.
2. If solving the GREOM system, update the particle 4-velocity with a Boris push in a local inertial frame over half a time step; if solving the GRGCA system, update the parallel 4-velocity in a local frame over half a time step.
3. Update the particle position and calculate the geodesic contribution to the 4-velocity (either the full u i for the GREOM or the parallel and drift components for the GRGCA) with an IMR scheme, iterating on the updated variables until convergence.
4. Apply a second half-update to the full particle 4-velocity (if solving the GREOM system) or to the parallel 4-velocity (if solving the GRGCA system) in a local frame.
5. After the update, store the new position and 4-velocity. If a GCA push was applied, retrieve x i ≈ R i and u i from the updated GCA variables (keeping the direction of the gyration 4-velocity fixed across a time step).

Numerical tests
In this Section we test the new coupled particle pusher, assessing the accuracy and computational cost of the algorithm. We choose to run our tests in two distinct setups, with and without the effect of strong gravity. We consider static (time-invariant) electromagnetic field configurations given by PiC simulations of an isolated current sheet (in flat spacetime; see Hakobyan et al. 2019) and of a current sheet in a BH magnetosphere (in Kerr spacetime; see Parfrey et al. 2019). We then push test particles (which do not provide feedback to the electromagnetic fields) in these configurations, comparing the particle trajectories produced by our coupled approach with the results given by solving the SREOM/GREOM system with the Boris pusher. Electromagnetic fields are evaluated at the particle position x via interpolation from N surrounding grid points, e.g., where W is a linear shape function. Gradients of grid quantities are obtained at the particle position by differentiating W with respect to x, e.g., ∇B(x) = N B N ∇W(x − x N ). In order to probe the dynamics of rapidly-gyrating particles, we consider a range of values for the ratio ρ C /L. This allows us to probe increasing scale separations between spatial field variations and particle gyro-motion.

Flat spacetime: Particle motion in reconnecting current sheets
Our electromagnetic field configuration is a fully developed reconnecting current sheet. This state represents a snaphsot from a PiC simulation of a magnetized pair plasma in flat spacetime (Hakobyan et al. 2019) with TRISTAN-MP (Buneman 1993;Spitkovsky 2005). The initial configuration consists of a Harris-like profile for the magnetic field in a two-dimensional box of size x × y = 5000d e × 5600d e (with grid spacing ∆x = 0.2d e ), where d e = c/ω p is the electron skin depth and ω p is the plasma frequency. The plasma magnetization in the ambient region, upstream of the current sheet, is σ = 100. Along the polarity inversion region at x = 2500d e , magnetic field lines break and reconnect, producing chains of plasmoids which travel in the y-direction. Particles travel toward the reconnection region with speed ∼ v E , gaining energy as they approach the X-line region. There, electric fields may grow significantly to strengths E > B, providing a source of strong particle acceleration.
In this setup, we initialize 5 positrons in the ambient region at x = 2540d e , ∼ 40d e away from the current sheet. The initial velocity is set to v = v E calculated with the value of the electromagnetic fields at the initial particle positions. The integration is carried out until t = 100ω −1 p . The original PiC run was parametrized such that particles have an initial gyro-radius ρ C,0 /d e = 10 (and gyro-frequency Ω C,0 /ω p ∼ 10) away from the current sheet. In our experiments, we find that for ρ C,0 /d e 10 −1 the Boris pusher can accurately capture the particle motion at all times, without the need for a switch to the GCA. Therefore, we explore a parameter range below this limit, ρ C,0 /d e = 10 −2 , 10 −3 , 10 −4 (corresponding to Ω C,0 /ω p = 10 2 , 10 3 , 10 4 ), in order to probe increasingly larger scale separations between spatial field variations and particle gyro-motion.
In order to obtain a reference solution, we run a first set of simulations with the standard Boris pusher, dynamically adapting the time step at each iteration as ∆t gyr = 2πΩ −1 C /60, such that the particle gyration is always resolved with at least 60 steps (hence away from the current sheet, ∆t gyr,0 0.1Ω −1 C,0 ). We then perform the same set of simulations with the Boris pusher, this time limiting the time step such that ∆t = 0.45ω −1 p (hence ∆t/∆t gyr,0 Ω C,0 /ω p ). In this way, we always enforce a time step that is comparable to that employed for the PiC simulation from which the field configuration is taken (obtained from a CFL-type condition). A final set of runs with the same criterion for the time step is performed with the coupled pusher (considering both curvature and E × B drifts). We note that, with this choice, the ratio ∆t gyr,0 /∆t decreases linearly as q/m increases, making the Figure 1: Left: trajectories of the ρ C,0 /d e = 10 −2 particles produced with the Boris pusher resolving the gyration (∆t = ∆t gyr , black lines and squares), and Boris pusher (cyan lines) and coupled pusher (red lines) with with ∆t/∆t gyr,0 10 2 . Green diamonds indicate, for the coupled pusher, the first switch between GCA and Boris schemes. Right: evolution in time of the Lorentz factor Γ and theρ C /∆x and E/B ratios (measured at the particle position), which determine the value of the switch S, for particles P3 and P4. For the coupled pusher, we show the value of S indicating the switch between the GCA (S = 0) and Boris (S = 1) algorithms. The results of the coupled scheme completely overlap with the reference Boris trajectories. The Boris scheme, instead, produces slightly diverging trajectories unless the gyration is accurately resolved at all times.
reference Boris runs progressively more expensive, at least as long as the particle is in the ambient region away from the current sheet.
To handle the switch between GCA and Boris pushers in the coupled approach, we employ the following criterion: at each time step, we measure the local value ofρ C = mcΓ/(|q|B) , i.e. the upper limit for a particle travelling with u ⊥g ∼ Γ. We then compare this value with a fraction of the grid spacing, f ρ (note that we take ∆ = ∆x). Additionally, we measure the local value of E/B as a proxy for the magnitude of v E , and we compare it to a constant factor f E . We therefore obtain a switch step function and apply a GCA integration step if S = 0 or a Boris push if S = 1. In the following we find that f ρ = 0.4 and f E = 1 produce the best results for the coupled pusher.
Our results are presented in Figures 1-2, where we show, for ρ C,0 /d e = 10 −2 and ρ C,0 /d e = 10 −4 and for all methods employed, the trajectory of all 5 particles (labeled at their initial positions as P1-P5) in space (left-hand panels). Green diamonds indicate, for each trajectory produced with the coupled pusher, the first switch between GCA and Boris schemes. The evolution in time of Γ,ρ C /∆x, and the E/B ratio for selected particles is shown in the right-hand panels. For the coupled pusher, we show the value of the switch function S in time in the bottom-right panels. The switch between the GCA and Boris algorithms is indicated by the values ofρ C /∆x and E/B falling above and below the f ρ and f E thresholds, respectively (dashed black lines in theρ C /∆x and E/B plots).
All particles travel from their initial position toward the current sheet, encountering spatial oscillations of the electromagnetic fields (caused by fast magnetosonic waves excited by plasmoid mergers; see Philippov et al. 2019) along the way. As a consequence, mild oscillations in the Lorentz factor arise, before the particles reach the reconnection region. There, particles can cross regions where strong electric fields E > B cause acceleration to Γ 1 (e.g., particle P4 in Figure 1 and P3 in Figure 2. In other cases, particles can fall inside plasmoids (visible as closed, circular loops along the current sheet), where they remain trapped around the central "O-point" (a magnetic null), experiencing mild acceleration (e.g., particle P3 in Figure 1). Figure 2: Left: trajectories of the ρ C,0 /d e = 10 −4 particles produced with the Boris pusher resolving the gyration (∆t = ∆t gyr , black lines and squares), and Boris pusher (cyan lines) and coupled pusher (red lines) with ∆t/∆t gyr,0 10 4 . Green diamonds indicate, for the coupled pusher, the first switch between GCA and Boris schemes. Right: evolution in time of the Lorentz factor Γ and theρ C /∆x and E/B ratios (measured at the particle position), which determine the value of the switch S, for particles P3 and P4. For the coupled pusher, we show the value of S indicating the switch between the GCA (S = 0) and Boris (S = 1) algorithms. The results of the coupled scheme completely overlap with the reference Boris trajectories. The Boris scheme, instead, produces substantially diverging trajectories unless the gyration is accurately resolved at all times.
The Boris pusher with ∆t/∆t gyr,0 Ω C,0 /ω p (cyan lines), which does not accurately capture gyration at all times, shows increasingly evident discrepancies with respect to the reference Boris results with ∆t = ∆t gyr (black lines/squares) as ρ C,0 /d e decreases. For ρ C,0 /d e = 10 −4 (i.e. Ω C,0 ∆t 10 4 ), the inaccuracy of the Boris pusher grows so large that it produces completely diverging particle trajectories. The results of the coupled pusher with ∆t/∆t gyr,0 Ω C,0 /ω p (red lines) are, instead, very consistent with the reference Boris solution in all cases. In our setup, all particles start with purely drifting motion (zero gyro-radius) and therefore their trajectories are initially captured by the GCA pusher. As they travel toward the current sheet, oscillations in theρ C /∆x and E/B ratios caused by fast magnetosonic waves can determine multiple switches between the GCA and Boris pushers. When particles reach the current sheet and cross E > B regions, they can experience a significant increase in Γ and therefore inρ C , which invalidates the GCA assumptions. At that point, the Boris scheme is capable of accurately resolving the particle motion, and the GCA pusher is no longer employed (e.g., particle P3 in Figure 2). Particles that approach o-points, instead, may not experience strong acceleration, such that the GCA pusher can be employed for longer times (e.g., particle P3 in Figure 1).
In a second ensemble of runs, we compare the reference Boris results with those of the coupled pusher with and without curvature terms (drift and acceleration) as illustrated in Sections 2.2.1 and 2.2.2. We exemplify our results for ρ C,0 /d e = 10 −2 and ρ C,0 /d e = 10 −4 in Figures 3-4, where we observe that, for the lowest ρ C,0 /d e ratio, neglecting curvature (orange lines) has a significant effect on several particle trajectories, which can diverge from the reference solution (black lines and squares). For such trajectories, quantifying the curvature acceleration terms in equation (9) reveals that these are comparable in magnitude to the parallel acceleration qE /m. Therefore, if curvature drifts are not taken into account, errors can arise in the particle trajectory from the solution of the GCA equations. These errors can then propagate and affect the results also when the Boris pusher is called after a switch. Instead, including curvature terms (magenta lines) reproduces the reference Boris results very accurately. As ρ C,0 /d e decreases, the importance of curvature effects in producing the correct trajectory decreases, and for ρ C,0 /d e = 10 −4 results obtained via the coupled pusher with and without curvature exhibit very small differences. This is expected, as the magnitude of the curvature terms (relative to zeroth-order terms) scales linearly with ρ C /L, and is therefore 100 times smaller for ρ C,0 /d e = 10 −4 . Therefore, the persisting minor differences in the trajectories where curvature drifts are not accounted for are also expected to be reduced as the magnetization is Figure 3: Left: trajectories of the ρ C,0 /d e = 10 −2 particles produced with the Boris pusher resolving the gyration (∆t = ∆t gyr , black lines and squares) and the coupled pushers with and without curvature drifts (magenta and orange lines respectively). Right: evolution in time of the Lorentz factor Γ and theρ C /∆x and E/B ratios (measured at the particle position), which determine the value of the switch S, for particles P2 and P5. For the coupled pushers, we show the value of S indicating the switch between the GCA (S = 0) and Boris (S = 1) algorithms. For this value of ρ C,0 /d e , several trajectories diverge from the reference solution when neglecting curvature terms in the coupled pusher.
further increased. We conclude that, for strongly magnetized plasmas in the upstream of current sheets, retaining v E as the only drift mechanism (as in the simple leapfrog GCA scheme of Section 2.2.2) suffices in producing accurately resolved particle trajectories with a coupled approach.
Our results show that, for large magnetizations (parametrized here via the ratio ρ C,0 /d e ) imposing small gyroradii, the coupled pushers retains very high accuracy even with a time step that may be several orders of magnitude larger than the gyro-period. This condition occurs primarily in regions away from current sheets, where particles essentially travel in the direction of E × B with negligible gyro-motion. In these regions, the coupled scheme is as accurate as the Boris pusher with ∆t = ∆t gyr (which fully resolves the gyration), but it allows for much larger (up to ∼ 10 4 times in our ρ C,0 /d e = 10 −4 run) time steps. Conversely, the results of the Boris pusher with such large time steps exhibit substantial inaccuracies.

Curved spacetime: Particle motion in a black-hole magnetosphere
Our electromagnetic configuration results from an initial axisymmetric Wald solution of Maxwell's equations in Kerr spacetime (Wald 1974). In this setup, a BH of mass M = 1 and spin a = 0.999 (in geometrized units where M = r g = 1) is embedded in an external magnetic field. Magnetic field lines thread the event horizon, stretching to infinity in the direction parallel to the BH spin axis. Due to frame-dragging effects, the solution includes a nonzero electric field. The magnetic field strength at infinity is set such that away from the central object, particles have an initial gyro-radius ρ C,0 /r g ∼ 10 −3 (corresponding to Ω C,0 /(c/r g ) ∼ 10 3 ). Starting from this initial state, a general-relativistic PiC simulation (in spherical Kerr-Schild coordinates) was run, where pair plasma is continuously injected (Parfrey et al. 2019). The simulation employed a computational grid covering the domain 0.985r g ≤ r ≤ 8r g , 0 ≤ θ ≤ π with 1280 2 cells. The evolution of the system leads to a quasi-steady state, where the plasma is globally force-free, and a thin current sheet of length ∼ r g forms in the equatorial plane. Analogously to the flat-spacetime case in the previous Section, particles entering the current sheet can be accelerated to relativistic energies by means of strong electric fields. The static (time-independent) electromagnetic fields taken from this configuration are employed to carry out our experiments.
We test our general-relativistic coupled scheme (which only includes E × B drifts in all cases) by initializing 5 electrons in the vicinity of the current sheet (at a distance ∼ r g from the equatorial plane). The particle velocity is Figure 4: Left: trajectories of the ρ C,0 /d e = 10 −4 particles produced with the Boris pusher resolving the gyration (∆t = ∆t gyr , black lines and squares) and the coupled pushers with and without curvature drifts (magenta and orange lines respectively). Right: evolution in time of the Lorentz factor Γ and theρ C /∆x and E/B ratios (measured at the particle position), which determine the value of the switch S, for particles P1 and P3. For the coupled pushers, we show the value of S indicating the switch between the GCA (S = 0) and Boris (S = 1) algorithms. For this value of ρ C,0 /d e , trajectories produced considering only E × B drifts show only minor differences with respect to those obtained when including curvature terms, indicating that E × B drifts are predominant. set equal to v E at the initial particle positions. We carry out the numerical integration until t ∼ 20r g /c, although most particles fall inside the event horizon well before this time. For consistency with the parameters employed in the previous Section, we explore the range ρ C,0 /r g = 10 −2 , 10 −3 , 10 −4 (Ω C,0 /(c/r g ) = 10 2 , 10 3 , 10 4 ). The choice of ∆t for each method is analogous to the flat-spacetime case, and we again obtain a reference solution by running the Boris scheme with a time step ∆t gyr = 2πΩ −1 C /60. The switch between the Strang split Boris and GCA schemes illustrated in Section 3 is handled similarly to the previous case, with a switch function S assuming values of 0 or 1 corresponding to a GCA or Boris push respectively. In this case, the grid spacing ∆ is computed from the volume integral (43), and we take f ρ = 1 and f E = 1.
Our results are presented in Figures 5-6 for ρ C,0 /r g = 10 −2 and ρ C,0 /r g = 10 −4 . We again show the trajectories for all particles with all three methods in the left-hand panels, and the time evolution of the Lorentz factor Γ,ρ/∆ , and D/B ratios (with dashed black lines indicating the f ρ and f E thresholds), and value of S (for the coupled pusher) for selected particles in the right-hand panels. In all cases, the particles initially drift along the D × B direction with zero gyro-radius. Particles that do not immediately fall into the event horizon end up in the current sheet, where they experience strong acceleration (up to Γ ∼ 10 3 for the smallest ρ C,0 /r g ). Differently from the flat-spacetime case, here fast magnetosonic waves from plasmoid mergers do not produce significant oscillations in E/B away from the current sheet, due to a smaller scale separation in the original global GRPiC simulation.
In agreement with the flat-spacetime test of Section 4.1 (see Figure 1), for ρ C,0 /r g = 10 −2 , we observe slight discrepancies between the trajectories produced with the coupled pusher (red lines with green diamonds marking the first switch between Boris and GCA) and the reference Boris results (black lines and squares). Consistently with the flat-spacetime case, we observe a progressive improvement in the results such that, for ρ C,0 /r g = 10 −4 , all particle trajectories are accurately captured by the coupled pusher. This indicates a transition to a sufficiently magnetized regime where D × B drifts dominate the particle motion away from the current sheet. In this regime, the inclusion of only v E in the GCA equations suffices in producing accurate trajectories. Conversely, as ρ C,0 /r g decreases the Boris scheme with ∆t/∆t gyr,0 Ω C,0 /(c/r g ) produces progressively larger inaccuracies in the results, which are completely diverging in the ρ C,0 /r g = 10 −4 case. We note that, for the coupled pusher, the thresholds f ρ = f E = 1 work well to activate the switch and to obtain accurate trajectories. For large charge-to-mass ratios, these thresholds act such that the GCA equations are employed until shortly before the particles reach the current sheet, where strong acceleration sharply raises theρ C /∆ ratio. At that point, the Boris pusher is capable of fully Figure 5: Left: trajectories of the ρ C,0 /r g = 10 −2 particles produced with the Boris pusher setting ∆t = ∆t gyr (black lines and squares), and the Boris pusher (cyan lines) and coupled pusher (red lines) with ∆t/∆t gyr,0 10 2 . Green diamonds indicate, for the coupled pusher, the first switch between the GCA and Boris schemes. Right: evolution of the Lorentz factor Γ, ratios for the switch functionρ C /∆ and D/B, and value of S (for the coupled pusher; S = 0 for GCA, S = 1 for Boris) in time for particles P1 and P5. Consistently with the special-relativistic case, for this value of ρ C,0 /r g the coupled pusher introduces slight inaccuracies in the particle trajectory, due to the fact that only D × B drifts are included in the GCA equations. Figure 6: Left: trajectories of the ρ C,0 /r g = 10 −4 particles produced with the Boris pusher setting ∆t = ∆t gyr (black lines and squares), and the Boris pusher (cyan lines) and coupled pusher (red lines) with ∆t/∆t gyr,0 10 4 . Green diamonds indicate, for the coupled pusher, the first switch between the GCA and Boris schemes. Right: evolution of the Lorentz factor Γ, ratios for the switch functionρ C /∆ and D/B, and value of S (for the coupled pusher; S = 0 for GCA, S = 1 for Boris) in time for particles P1 and P4. For this case, the coupled pusher reproduces the reference results very well, while the Boris algorithm produces large inaccuracies unless the gyro-motion is a accurately resolved with small time steps ∆t gyr = 2πΩ C /60 at all times. Figure 7: Trajectories of 4 particles (ρ C,0 /r g = 10 −1 , 10 −2 , 10 −3 , 10 −4 ) in the ambient region away from the central BH. The Boris pusher with ∆t = ∆t gyr (black lines and squares) provides reference results where, for large magnetizations (i.e. ρ C,0 /r g < 10 −1 ), particles travel along straight magnetic field lines with negligible gyro-motion. When ∆t/∆t gyr,0 Ω C,0 /(c/r g ), the Boris pusher (cyan lines) remains accurate only as long as ρ C,0 /r g ≥ 10 −1 ; for larger magnetizations, the small-scale gyration is not accurately resolved at all times, and the scheme introduces large inaccuracies in the trajectories, which substantially diverge from the reference solutions. In the coupled pusher (red lines), instead, small-scale gyration always triggers a call to the GCA solver only, and the trajectories perfectly overlap with the reference results. For the case ρ C,0 /r g = 10 −1 , where the gyro-radius remains well above the grid spacing, we overplot the guiding center trajectory (green line) for reference. resolving the particle motion, and numerical errors caused by overstepping the gyro-period are avoided.
In a final set of runs, we consider particle motion in the ambient region (outside of the region where the BH rotation twists magnetic field lines, r sin θ > 2r g ), where highly magnetized particles travel along (essentially) straight field lines with almost constant velocity and negligible gyro-motion. For this test, we initialize 4 electrons such that ρ C,0 /r g = 10 −1 , 10 −2 , 10 −3 , 10 −4 (i.e. Ω C,0 /(c/r g ) = 10 1 , 10 2 , 10 3 , 10 4 ). All particles are initially placed at r = 5.5r g , θ = π/2 with zero velocity, and we carry out the integration until t = 20r g /c. The thresholds f ρ = f E = 1 are employed for the coupled pusher.
The results are presented in Figure 7, where we show part of the reference Boris trajectories (black lines and squares) obtained by accurately resolving the gyro-motion, and the results of the Boris and coupled schemes with ∆t/∆t gyr,0 Ω C,0 /(c/r g ) (cyan and red lines respectively). For ρ C,0 /r g = 10 −1 , the gyro-radius is well above the grid spacing, such that the gyro-motion is clearly visible. In this case, the Boris scheme with ∆t/∆t gyr,0 Ω C,0 /(c/r g ) can accurately resolve the gyration, hence reproducing the reference solution very well. The coupled scheme performs equally well, never switching to the GCA equations, and overlaps the Boris solution perfectly. As a reference, for this case we overplot the guiding center trajectory produced by integrating the GRGCA equations (green line). As the magnetization increases such that ρ C,0 /r g ≤ 10 −2 , the gyro-radius progressively shrinks below the grid spacing, therefore requiring extremely small time steps to accurately resolve the gyration. In the Boris scheme with ∆t/∆t gyr,0 Ω C,0 /(c/r g ) the gyro-period becomes significantly overstepped, introducing progressively larger inaccuracies which result in completely diverging trajectories. In the coupled scheme, instead, the subgrid gyro-motion is detected via the f ρ threshold, such that the Boris pusher is never employed, and the trajectory is calculated via the GCA equations at all times. As a consequence, the results of the coupled pusher reproduce the reference trajectories extremely well also for large magnetizations, with no signs of the inaccuracies that characterize the under-resolved Boris results. The case ρ C,0 /r g = 10 −4 does not represent a limiting value for the applicability of the coupled scheme: in our experiments, we further pushed the scale separation to values ρ C,0 /r g < 10 −7 , with completely analogous results, indicating that the coupled scheme is capable of retaining high accuracy even for gyro-motion on extremely small scales.
Our tests confirm that, for strongly magnetized plasmas and even for scenarios where strong gravity plays an important role, the coupled pusher can produce reliable physical results, capturing the motion of particles with gyro-radius much smaller the grid spacing. Conversely, the standard Boris pusher produces large inaccuracies that can only be mitigated by taking much smaller time steps (we observe large inaccuracies for Ω C,0 ∆t 10), increasing the computational requirements dramatically (up to a factor ∼ r g /ρ C,0 larger than the cost of the coupled pusher).

Conclusions
Plasma in collisionless magnetospheres of NSs and BHs is usually highly magnetized, such that plasma particles gyrate with gyro-frequencies that are orders of magnitude larger than the frequency of rotation of the compact object. In this regime, particle gyration becomes negligible, and the particle motion practically reduces to that of the guiding center. However, in dissipation regions, e.g., magnetic reconnection layers where the local magnetic field strength goes to zero, particles can accelerate and gain significant transverse momentum resulting in non-negligible gyro-radii. These conditions put stringent constraints on the time step of numerical algorithms for plasma simulations of compact-object magnetospheres, and the system's physical parameters are in practice typically rescaled to artificially decrease the scale separation.
To mitigate this problem we designed a novel coupled particle pusher that can accurately resolve particle dynamics both in highly magnetized plasmas and in acceleration regions where the local magnetic field strength decreases. The algorithm dynamically switches between solving the full equations of motion with a Boris scheme, and the guiding center equations of motion with a newly developed semi-implicit leap-frog scheme, depending on the magnetization of the plasma particles. In magnetized regions we solve only for the guiding center motion of the orbit, averaging over the particle's gyration. In this way the restriction on the time step due to very small gyro-frequency is alleviated, while the particle dynamics is accurately resolved both for large and small gyro-radii. This is appropriate for example for typical field strengths near pulsars and magnetars where synchrotron losses very quickly push particles to the zeroth Landau level.
We tested the new algorithm for typical plasma parameters in compact-object magnetospheres. We considered two cases, namely i) an isolated reconnecting current sheet in magnetized plasma, and ii) a reconnection layer that formed in the magnetized magnetosphere of a BH. We compared the motion of charged particles with our new coupled pusher to a reference solution of the full equations of motion where the gyro-frequency is fully resolved, for a range of particle magnetizations. In this way, we probed increasingly larger scale separations between the particle gyro-motion and the length scales of variation of the electromagnetic fields (which are resolved on the computational grid). We find that the new coupled scheme can reproduce the reference results very well for the large scale separations we considered here, as applicable in magnetospheres of compact objects. We also show that for particles moving in the magnetized upstream (away from current sheets) all drifts associated with field curvature can be safely neglected for the strong magnetic fields near compact objects. Our particle pusher is presented such that it can be readily implemented in existing PiC algorithms, both in flat spacetime and in a dedicated covariant formulation for curved spacetime. The scheme can also be used to trace test particles in magnetohydrodynamic simulations, for example to study particle acceleration due to resistive relativistic reconnection in compact-object magnetospheres, accretion disks, coronae, and jets Ripperda et al. 2019aRipperda et al. ,b, 2020. In future applications we will focus on global simulations of particle acceleration in current sheets that form in highly magnetized compact-object magnetospheres. In these environments, when particles accelerate to high energies, they are typically subject to radiative (e.g. synchrotron) losses. In order to probe the fast-cooling regime of energetic particles in the current sheet (Cerutti et al. 2016;Philippov & Spitkovsky 2018), previous PiC simulations of pulsar magnetospheres required to take very small time steps such that the Larmor frequency is resolved everywhere, including regions near the stellar surface where gyration can become extremely fast2. Our new coupled pusher removes this limitation, as only high-energy particles, pushed with a Boris algorithm, are subject to synchrotron cooling.
As a way forward to truly multi-scale simulations of relativistic magnetospheres of compact objects, a particularly attractive approach for these magnetized plasmas is to couple a PiC method with a force-free description, where magnetic fields are considered so strong that the inertia of the plasma can be ignored; highly magnetized ideal (i.e., where E · B = 0) parts of the computational domain can be treated with a simple force-free algorithm, while particle acceleration (and the associated feedback on electromagnetic fields) in dissipation regions, like current 2This is necessary because accurate accounting of the synchrotron losses requires resolving the Larmor frequency in time. Under-resolving Larmor gyrations in time, while keeping the same synchrotron cooling strength, results in numerical instabilities.
sheets, is accurately captured by a PiC method. While our coupled GCA-Boris approach alleviates constraints on the time step based on the gyro-radius and gyro-frequency, the skin depth and associated plasma frequency still require to be resolved in explicit PiC simulations. Ultimately, an implicit PiC approach can alleviate constraints to resolve the plasma frequency when dynamics take place on scales much larger than the skin depth (Bacchini et al. 2019a). Such an approach can be combined with an adaptive or stretched mesh, that refines particularly in dissipation regions, while remaining coarse in the magnetized upstream. We will work on developing these new numerical methods in the future.