Issue 
A&A
Volume 666, October 2022



Article Number  A5  
Number of page(s)  18  
Section  Astrophysical processes  
DOI  https://doi.org/10.1051/00046361/202243634  
Published online  27 September 2022 
Particle acceleration and radiation reaction in a strongly magnetised rotating dipole
Université de Strasbourg, CNRS, Observatoire astronomique de Strasbourg, UMR 7550, 67000 Strasbourg, France
email: jerome.petri@astro.unistra.fr
Received:
25
March
2022
Accepted:
30
June
2022
Context. Neutron stars are surrounded by ultrarelativistic particles efficiently accelerated by ultrastrong electromagnetic fields. These particles copiously emit highenergy photons through curvature, synchrotron and inverse Compton radiation. To date, however, no numerical simulations have been able to handle such extreme regimes of very high Lorentz factors and magnetic field strengths close to or even above the quantum critical limit of 4.4 × 10^{9} T.
Aims. It is the purpose of this paper to study particle acceleration and radiation reaction damping in a rotating magnetic dipole with realistic field strengths of 10^{5} T–10^{10} T typical of millisecond and young pulsars and of magnetars.
Methods. To this end, we implemented an exact analytical particle pusher including radiation reaction in the reduced Landau–Lifshitz approximation where the electromagnetic field is assumed constant in time and uniform in space during one time step integration. The position update is performed using a velocity Verlet method. We extensively tested our algorithm against time independent background electromagnetic fields like the electric drift in cross electric and magnetic fields and the magnetic drift and mirror motion in a dipole. Finally, we apply it to realistic neutron star environments.
Results. We investigated particle acceleration and the impact of radiation reaction for electrons, protons, and iron nuclei inserted around millisecond pulsars, young pulsars, and magnetars, in comparison to situations without radiation reaction. We found that the maximum Lorentz factor depends on the particle species, but only weakly on the neutron star type. Electrons reach energies up to γ_{e} ≈ 10^{8} − 10^{9}, whereas protons reach energies up to γ_{p} ≈ 10^{5} − 10^{6} and iron up to γ ≈ 10^{4} − 10^{5}. While protons and iron are not affected by radiation reaction, electrons are drastically decelerated, reducing their maximum Lorentz factor by four orders of magnitude. We also found that the radiation reaction limit trajectories agree quite well with the reduced Landau–Lifshitz approximation in almost all cases.
Key words: magnetic fields / methods: analytical / stars: neutron / stars: rotation / pulsars: general
© J. Pétri 2022
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the SubscribetoOpen model. Subscribe to A&A to support open access publication.
1. Introduction
Neutron stars are known to harbour ultrastrong magnetic fields close to or even above the quantum critical limit of B_{c} ≈ 4.4 × 10^{9} T. The subclass of magnetars usually sustains field strengths well above this value of B_{c}. These stars are therefore able to accelerate leptons and hadrons to extremely relativistic regimes of very high Lorentz factors γ ≈ 10^{9}. In such an extreme environment, the radiation reaction is expected to drastically perturb their trajectory compared to the pure Lorentz force motion. Highenergy and very highenergy photons are produced and sometimes detected on Earth by Cerenkov telescopes.
To date, however, a quantitatively accurate study of these acceleration and radiation reaction mechanisms has failed due to the inability of current numerical algorithms to handle such strong fields. The problem is circumvented by artificially decreasing the magnetic field strength and other relevant physical parameters like the Lorentz factor, and meanwhile increasing the associated Larmor radius. Unfortunately, the high nonlinearity of the problem renders any extrapolation to realistic fields risky. The only satisfactory results must come from faithful simulations employing appropriate lengths and timescales observed around neutron stars.
The combination of strong fields and large Lorentz factors leads naturally to strong radiation reaction damping of the charged particle motion. These trajectories have been computed in the past for test particles, for instance by Finkbeiner et al. (1989) in the pulsar vacuum field. Finkbeiner et al. (1990) discussed the validity of the Lorentz–Dirac equation and the Landau–Lifshitz approximation used in such computations. Herold et al. (1985) integrated the equation of motion with the radiation reaction in the ultrarelativistic regime, and showed the difference between radiative damping and no damping for an aligned rotator. They also gave an estimate of the maximum Lorentz factor.
Exact analytical solutions of the Landau–Lifshitz equation have been found for the arbitrary plane waves reported by Piazza (2008) and Hadad et al. (2010). For constant and uniform electromagnetic fields, solutions have been known since the work of Heintzmann & Schrüfer (1973). These are special solutions found by removing the temporal and spatial derivatives from the Landau–Lifshitz approximation. This simplified version is sometimes called the reduced Landau–Lifshitz equation (LLR). We use this approximation to advance in time the position and velocity of charged particles.
Pushers based on exact analytical solutions have been implemented by several authors. For instance Laue & Thielheim (1986) evolved particles in an orthogonal magnetic dipole, whereas Ferrari & Trussoni (1974) investigated particle motion in a dipole field, neglecting the displacement current. Recently Pétri (2020) developed an algorithm to evolve particles in a strong electromagnetic field. Tomczak & Pétri (2020) applied it to a magnetic dipole associated with strongly magnetised rotating neutron stars. Gordon et al. (2017b,a) showed how to implement a fully covariant particle pusher and gave some hints that the radiation reaction should be included. Later Gordon & Hafizi (2021) developed a special unitary pusher for extreme fields achieving computation costs comparable to those for the Boris algorithm (Boris 1970).
In the ultrarelativistic regime, the radiation reaction almost exactly balances the electric field acceleration leading to a particle velocity only depending on the local electromagnetic field configuration. As shown by Mestel et al. (1985), the Lorentz factor can then be deduced from the trajectory curvature. Kelner et al. (2015) carefully studied the synchrocurvature radiation of ultrarelativistic particles evolving in a strongly curved electromagnetic field. The pitch angle plays a central role in controlling the synchrotron versus curvature regime.
Several different but not equivalent approaches have been designed to include the radiation reaction in a particle pusher for ultrastrong electromagnetic fields. Vranic et al. (2016) offers a comprehensive study of the most widely used techniques to implement the radiation reaction force in standard Lorentz force pushers. However, the numerical algorithms that explicitly solve the Landau–Lifshitz equation have some issues in satisfying conservation laws for long runs. Nevertheless, timesymmetric implicit methods seem to give better results (Elkina et al. 2014). Interestingly, exact analytical solutions of the reduced Landau–Lifshitz equation were found several decades ago by Heintzmann & Schrüfer (1973) for a constant electromagnetic field. These expressions are used by Li et al. (2021) for implementation in a particleincell (PIC) code following a projection onto an electric and a magnetic subspace (Boghosian 1987). Pétri (2021) also applied this exact solution to the acceleration of particles in a lowfrequency strong amplitude electromagnetic plane wave such as that launched by a strongly magnetised rotating neutron star.
In this paper we study particle acceleration in a realistic neutron star environment, using the exact scaling between the neutron star spin and the cyclotron frequency. In Sect. 2 we recall the equation of motion as derived by Landau–Lifshitz and its exact analytical solution, the appropriate normalisation, and the algorithm. Section 3 presents extensive tests of our algorithm in static fields showing its second order in time convergence. Section 4 describes an astrophysical application to neutron star electrodynamics and the upper limit of particle acceleration efficiency. Section 5 compares the radiation reaction limit regime to the exact motion. Finally, conclusions are drawn in Sect. 6.
2. Equation of motion
The selfforce produced by an accelerated charge is usually described by the LorentzAbrahamDirac (LAD) equation (Abraham 1902, 1904; Lorentz 1916; Dirac 1938). Unfortunately, this selfforce leads to runaway solutions because the associated equation of motion is of third order in time. Several remedies have been found to remove these unacceptable solutions (see e.g., Rohrlich 2007 for discussions). One approach often quoted in the literature is the Landau–Lifshitz formulation, a perturbative expansion of the LAD equation (Landau & Lifchitz 1989). In the remainder of this paper we adopt this point of view.
2.1. Landau–Lifshitz approximation
In order to eliminate the LAD flaw, Landau & Lifchitz (1989) derived an approximation valid in most configurations found in astrophysical applications. This new equation of motion is free of runaway instabilities and is largely employed in the plasma community. Their formulation leads to the equation of motion
where q and m are the particle charge and rest mass, u^{i} its fourvelocity, τ its proper time, F^{ik} the electromagnetic or Faraday tensor, c the speed of light, and τ_{m} the light crossing time across the particle classical radius r_{m} (within a factor unity)
It is advantageous to express it in terms of the electron classical radius r_{e} crossing time amounting to
The typical timescale for the radiation reaction is therefore
For instance, for protons this time is three orders of magnitude less than for leptons:
Interestingly, exact analytical solutions have been computed for Eq. (1) in some special configurations of electromagnetic fields, time dependent or time independent. We succinctly recall the useful results required for the present work.
2.2. Exact analytical solutions
An exact solution for LLR is based on the eigensystem expansion of the electromagnetic tensor . Earlier results were given by Heintzmann & Schrüfer (1973). Here we follow the notation of Li et al. (2021). Starting from the Lorentz force written as
where the electromagnetic tensor F has been replaced by G = q F/m to absorb the charge over mass ratio, we decompose the fourvelocity u into a magnetic and an electric part, denoted respectively as u_{B} and u_{E}, such that u = u_{E} + u_{B}. The real eigenvalues of are ±λ_{E}, whereas the imaginary eigenvalues are ±i λ_{B}, because λ_{E} and λ_{B} are real and positive numbers, with dimensions similar to pulsation thus in 1/s. Then each vector u_{E} and u_{B} remains in an eigensubspace satisfying
The vector components u_{E} and u_{B} are obtained by defining the projection operators onto the subspaces E and B by (Boghosian 1987)
where I is the identity matrix. These operators are well defined only if . If both electromagnetic invariants vanish, we retrieve a nulllike field, which requires a different treatment, as given for instance by Pétri (2021). In the non nulllike field we get
The equation of motion decouples into two parts given by
The exact analytical solutions with initial conditions and are
Adding the radiation reaction in the LLR limit leads to the exact expression
with . These expressions are similar to the original formulas found by Heintzmann & Schrüfer (1973). The radiation reaction effect becomes perceptible after a time τ ≈ 1/α. The component u_{E} is associated with the accelerating motion induced by the electric field, whereas the u_{B} component is related to the gyromotion in the magnetic field. When α vanishes, the radiation reaction effect disappears. The denominators in u_{E} and u_{B} reduce to unity and the solutions to the Lorentz force fourvelocity components are recovered.
2.3. Normalisation
The relevant physical parameters determining the particle trajectory is decided through normalisation procedures that imply the following useful quantities in order to write the equation of motion without dimensions. These primary fundamental variables are the speed of light c, a typical frequency ω involved in the problem, the particle electric charge q, and the particle rest mass m.
From these quantities we derive a typical timescale and length scale as well as electromagnetic field strengths such that the length scale L_{0} = c/ω, the timescale T_{0} = 1/ω, the magnetic field strength B_{0} = m ω/q, and the electric field strength E_{0} = c B_{0}. Normalised quantities will be overlaid with a tilde symbol.
The two important parameters defining the family of solutions are the field strength parameters a_{B} and a_{E} and the radiation reaction efficiency ω τ_{m} according to the following definitions:
Introducing the weighted and normalised electromagnetic field tensor by and a normalised time , the LandauLifshitz Eq. (1) is rewritten without dimensions as
The normalised and reduced Landau–Lifshitz equation reads
The particle fourvelocity depends only on the strength parameters a_{B} and a_{E} and on the radiation reaction strength parameter b. Therefore, it is unnecessary to compute trajectories for different particles possessing the same numbers a_{B}, a_{E}, b. The only differences are reflected in the physical timescale and space scales involved.
Generally speaking, we admit that the radiation reaction is negligible whenever the timescale of damping, given by 1/α, becomes larger than the characteristic timescale of our system, which is 1/ω. Expressed in quantities without dimension, we get . Therefore, the relevant parameter to quantify radiation reaction is not b, but the combination of b and the strength parameters a_{B} and a_{E}. Specific examples are given in the test Sect. 3.
2.4. Algorithm
For the remainder of this paper we use a Cartesian coordinate system (x, y, z) and the corresponding Cartesian orthonormal basis (e_{x}, e_{y}, e_{z}).
The velocity vector is integrated analytically following the previous discussion. Unfortunately, there is no simple analytical expression for the position vector, although some formulas can be found involving hypergeometric _{2}F_{1} functions with complex arguments (see Sect. 3 for an example in a constant magnetic field). The update in particle position is therefore performed by the velocityVerlet algorithm:
The superscript n refers to the proper time τ^{n} = n Δτ, and the same for the halfinteger superscript τ^{n + 1/2} = (n + 1/2) Δτ. We found this method more robust than the full analytical update in velocity and position. For particles trapped in a dipole magnetic field, undergoing bouncing motion with banana orbits typical of magnetic confinement devices for thermonuclear fusion reactors or in the Earth’s magnetosphere known as the Van Allen belt, the stability and convergence properties of the velocityVerlet algorithm is superior.
Before using our code to compute particle acceleration and radiation in the ultrastrong electromagnetic field of a dipole rotating in a vacuum, we test it against the exact analytical solutions in simple geometric configurations, but with very high Lorentz factors and/or very high fields. The results will also be compared to the radiation reaction limit regime, which is much less time consuming from a computational point of view, but also less accurate in some configurations (Sect. 5).
2.5. Radiation reaction limit
In ultrastrong electromagnetic fields, such as those present around neutron stars, radiation reaction plays an important role. In the asymptotic limit of ultrarelativistic motions, assuming that the radiation damping exactly balances the electric field acceleration, there is a simple analytical expression for the particle velocity depending only on the local values of the fields (Mestel et al. 1985). This velocity is decomposed into an electric drift motion, interpreted as the velocity required to switch to a frame where the electric and magnetic fields are aligned, and a motion along this common direction in this new frame. Denoting the velocity vector for positive charges as v_{+} and that for negative charges as v_{−}, we find
which corresponds to particles moving exactly at the speed of light. Here E_{0} and B_{0} are the strength of the electric and magnetic field in the frame where they are aligned. They are obtained from the electromagnetic invariants and ℐ_{2} = c E ⋅ B = c E_{0} B_{0}. Imposing E_{0} ≥ 0 we find
We compare the simulation results obtained from this simple prescription with the exact integration of the equation of motion according to LLR.
Applying this radiation reaction limit to neutron star magnetospheres, the velocity in Eq. (17) can be slightly simplified because of the presence of a plasma, the parallel electric field component (with respect to the magnetic field direction) being efficiently screened. In such a configuration, ℐ_{2}≪ℐ_{1} and ℐ_{1} < 0. The velocity then reduces to
The first term corresponds to the electric drift speed, whereas the second term is associated with the motion along the magnetic field lines, the particle gyromotion being absent in this picture. We note that the velocity component along the magnetic field reverses sign when crossing a point where E ⋅ B changes sign. These regions are able to trap particles depending on their charge and on the (E ⋅ B) configuration in the neighbourhood of this surface (Finkbeiner et al. 1989).
Several limiting cases are also useful to discuss. First, if the electric field vanishes, E_{0} = 0, the radiated power also vanishes, and the particle moves along the field lines with v_{±} = ±c B/B. Second, if the electric field is orthogonal to the magnetic field, E ⋅ B = 0 and E < c B, the particle motion is decomposed into an electric drift and a motion along B such that
This expression holds well within the lightcylinder of a forcefree magnetosphere.
3. Tests
We checked our algorithm against simple electromagnetic field configurations containing only an electric field, a magnetic field, or a cross electromagnetic field. Although the exact solutions are simple expressions, from a numerical point of view it is of paramount importance to ensure that the code is able to handle very high strength parameters and Lorentz factors such as those found in neutron star magnetospheres, that is about a_{B} ≈ 10^{20} and γ ≈ 10^{10}. Our main purpose in this section is to check that the results are not affected by roundoff errors.
3.1. Constant electric field
In a constant electric field, a charged particle is permanently accelerated in the direction of the electric field while it loses energy. Specialising the general solution (12) to a pure electric field aligned with the zaxis such that E = E e_{z}, we obtain
with the initial fourvelocity component along E, u_{⊥} the initial fourvelocity component perpendicular to E, and .
For a particle starting at rest, u_{∥} = u_{⊥} = 0 and γ_{0} = 1, the fourvelocity simplifies drastically to
This fourvelocity does not depend on the radiation reaction intensity. It accelerates as if it only experiences the Lorentz force. This peculiar situation is well known, and is discussed at length by Fulton & Rohrlich (1960) for a charge and its related classical radiation in a uniformly accelerating field.
The fourposition is given by introducing the two complex functions with the help of the hypergeometric functions _{2}F_{1} (Olver 2010) such that
Then the time and position are given by
with C_{0}, C_{2}, C_{3} complex constants of integration to satisfy the initial conditions.
Returning to Eq. (21), the typical electric acceleration timescale is τ_{acc} ∼ 1/ω_{E}. On the other hand, the radiation damping timescale is τ_{rad} ∼ 1/α. The ratio of the two timescales is therefore τ_{acc}/τ_{rad} ∼ τ_{m} ω_{E} = b. As expected, for small damping parameters b ≪ 1, the acceleration time is much shorter than the radiative damping and the particle is accelerated as if it would not radiate, until the time τ_{rad} ∼ τ_{acc}/b ≫ τ_{acc}. We note that this rough estimate needs to be corrected by taking into account the initial Lorentz factor, as discussed below.
To perform the simulations we use the characteristic frequency ω_{E} as normalisation, which leads to a normalised proper time . Therefore, the only relevant parameter apart from the initial conditions is b = τ_{m} ω_{E} and . Particles starting at rest or possessing an initial velocity directed along the electric field do not suffer from the radiative force. Consequently, as a typical example, particles start with an initial velocity perpendicular to the electric field, meaning u_{∥} = 0 and u_{⊥} ≠ 0. We chose different initial Lorentz factors in the set log γ_{0} = {0, 4, 8}. The damping factor is given by log b = {0, −5, −10, −15}. As output of the simulations, we plot the Lorentz factor increasing according to Eq. (21a) and shown in Fig. 1. For a particle starting at rest, whatever the damping parameter b, the solution is always equal to Eq. (22) with an acceleration arising around the time ω_{E} τ ∼ 1. This configuration is particular and is not impacted by the radiation reaction. More interesting are the particles starting with a substantial kick velocity and high Lorentz factors γ_{0} ≫ 1. The time derivative of the Lorentz factor is always negative for γ_{0} > 1 because
Fig. 1. Increase in the Lorentz factor due to radiation reaction for different initial Lorentz factors log γ_{0} = {0, 4, 8} and different damping factors log b = {0, −5, −10, −15}. The vertical lines show the time when the damping sets in before the electric acceleration phase starts. The coloured points show the simulation results and the black solid lines correspond to the exact analytical solutions. 
meaning that the particle first decelerates due to the radiative friction. At large times, when α τ ≫ 1 and ω_{E} τ ≫ 1, the Lorentz factor behaves as γ(τ)≈cosh(ω_{E} τ), losing its information about the initial state. It resembles the motion of a particle starting at rest, independently of γ_{0}. This occurs because the perpendicular motion is strongly damped, , and only the parallel velocity u_{∥} survives at large times with . In between, the normalised time remains small and the Lorentz factor can be approximated by
Therefore, before the acceleration phase starts, there is a deceleration step arising at time . These values agree with the curves in Fig. 1. If , the radiation reaction force has no time to set in and the motion tends to a purely accelerated regime given by Eq. (22). This is the case with γ_{0} = 10^{4} and b = 10^{−10} (orange dots) or γ_{0} = 10^{8} and b = 10^{−15}, which is just on the edge of this condition, showing a weak deceleration right before the electric boost (blue dots). The perpendicular momentum decrease is not necessarily significant before the acceleration; it is controlled by b and γ_{0} because at time ω_{E} τ ≈ 1 it braked to a momentum
Thus, the radiation reaction again impacts the motion if . Consequently, it is the combination that controls the damping efficiency, not b alone found from the simple arguments above.
Because the fourposition of the particle is computed numerically and not analytically according to Eq. (24), it is important to estimate the convergence rate of our scheme. To this end, Fig. 2 shows the error in the y and z positions and time t with decreasing proper time step Δτ for log γ_{0} = 4 and log b = −5 (in blue, orange, green, and red, respectively). The secondorder expectations are depicted by the green line. We conclude that the decrease in the relative error follows a second order in time scheme, as expected from the velocityVerlet algorithm shown in Sect. 2.
Fig. 2. Relative error of the position y, z and time t shown in the legend. The error decreases with second order in Δτ as given by the green line Δτ^{2} for log γ_{0} = 4 and log b = −5. The t and z errors overlap and are indistinguishable. 
3.2. Constant magnetic field
A charged particle orbiting in a constant magnetic field loses energy and decays until it rests. The rate of decay is controlled by the magnetic field strength only. The exact solution for the fourvelocity in a magnetic field directed along the zaxis with B = B e_{z} is given by
with the initial fourvelocity component along B, u_{⊥} the initial fourvelocity component perpendicular to B, and . Apart from the change in the gyrofrequency, the magnetic field strength impacts only the timescale for the decay via the exponential terms of arguments 2 α τ.
To perform simulations we use the characteristic frequency ω_{B} as normalisation with . Therefore, the only relevant parameters, apart from the initial conditions, are b = τ_{m} ω_{B} and . The length scale is therefore given in units of the nonrelativistic Larmor radius
Integrating the fourvelocity vector, an exact analytical expression for the fourposition is computed with the help of the hypergeometric functions _{2}F_{1}. Introducing the complex functions
the solution reads
where the C_{i} with i ∈ [0…2] are complex constants of integration to satisfy the initial conditions.
The particle trajectory follows a spiral, as shown in Fig. 3. The particle comes to rest after a typical time ω_{B} τ_{∞} ≫ 1/b.
Fig. 3. Particle orbit in a uniform and constant magnetic field and subject to radiation reaction. The initial Lorentz factor is log γ_{0} = 4. The inset shows the strong damped motion in green and even stronger damping in red where the spiralling is not seen. 
The corresponding Lorentz factor decreases according to Eq. (28a) and is shown in Fig. 4. The time when damping sets in is given approximately by . These times are shown as coloured vertical lines in Fig. 4. If the particle moves along the field line, it experiences no damping and keeps a uniform motion.
Fig. 4. Decrease in the Lorentz factor due to radiation reaction in a uniform and constant magnetic field associated with the orbits shown in Fig. 3. The coloured points show the simulation results and the black solid lines correspond to the exact analytical solutions. 
A comparison between the analytical trajectory (red solid line) and the numerical integration (blue dots) is shown in Fig. 5. A more quantitative agreement is proven in Fig. 6 where the relative error decreases with respect to the proper time step Δτ. Here the method is also second order in time, as expected.
Fig. 5. Comparison between the analytical solution, Eq. (31) (red solid line) and the numerical simulation (blue dots) for log γ_{0} = 4 and log b = −5. 
Fig. 6. Relative error of the position x and y as shown (see inset for legend). The error decreases with second order in Δτ, as given by the violet line Δτ^{2} for log γ_{0} = 4 and log b = { − 5, −10}. 
3.3. Cross electric and magnetic field
The cross electric and magnetic field configuration is a stringent test for an ultrarelativistic particle pusher. If the electric field strength is less than the magnetic field strength E < c B, then an appropriate Lorentz transform brings the problem to a frame where the electric field vanishes. We therefore return to the situation described in the last section with a constant and uniform magnetic field. For sufficiently long time the only remaining motion is the electric drift at speed v_{E} = E ∧ B/B^{2}. Therefore, the velocity is β_{E} = v_{E}/c = E/cB and the corresponding final Lorentz factor .
We performed simulations with E/cB = 0.999 and initial Lorentz factors log γ_{0} = {0, 4, 8}. The final Lorentz factor is γ_{∞} ≈ 22.3. Figure 7 shows the Lorentz factor (coloured dots) compared to the analytical expression (solid black lines). The agreement is excellent and demonstrates the high efficiency of our algorithm to capture ultrarelativistic motion with high precision.
Fig. 7. Decrease in the Lorentz factor due to radiation reaction in a cross electric and magnetic field. The coloured points show the simulation results and the black solid lines correspond to the exact analytical solutions. 
The quantitative agreement is checked by transforming the trajectory to the electric drift frame denoted by the coordinates (x′,y′). In this frame the orbital radius is decreasing, as shown in Fig. 8 for log γ_{0} = 4 and log b = −5, using different proper time steps such as log(ω_{B} Δτ) = {−1, −2}. The analytical solution found from the appropriate parameters in Eq. (31) is overlapped in red solid lines.
Fig. 8. Orbit in the electric drift frame (x′,y′) for log γ_{0} = 4 and log b = −5 for different proper time steps log(ω_{B} Δτ) = {−1, −2}. 
Figure 9 shows the relative error in the x′ and y′ position depending on the proper time step Δτ. The scheme converges to second order in proper time step.
Fig. 9. Relative error of the particle position associated with the electric drift motion shown in Fig. 8. 
3.4. Parallel electric and magnetic field
Another interesting configuration not reducible to any of the previous configurations are the parallel electric and magnetic fields. In this case the second electromagnetic invariant does not vanish ℐ_{2} ≠ 0. Consequently, there is no frame where either the electric or magnetic field vanishes. The electric and magnetic velocity components u_{E} and u_{B} decouple into an acceleration along the common direction and a gyration in the same direction. Assuming this direction to be e_{z}, the fourvelocity reads
We recognise the special cases of a pure electric field for u^{t}, u^{z} and a pure magnetic field for u^{x}, u^{y}, the only difference being the value of , including nonvanishing electric and magnetic contributions.
After a transition time, the gyromotion is significantly damped and the electric acceleration directs the velocity along its direction. The initial conditions are washed out and the particle moves as in a constant electric field with an almost constant acceleration leading to a hyperbolic motion, well known in special relativity kinematics. We estimate the duration of this transient stage. Either the particle is drastically accelerated before the gyration is damped or instead the orbit shrinks significantly before the electric field sensibly accelerates the particle. The situation depends on ordering of the eigenvalues λ_{E} and λ_{B}.
Figure 10 shows the analytic evolution of the Lorentz factor for a particle starting with only a perpendicular velocity component such that log γ_{0} = 8. The damping factor is set to log b = {0, −5, −10, −15} and the electric field strength E_{0} is varied relative to B_{0} such that log(E_{0}/c B_{0}) = {−2, 0, 2} (solid lines, dashed lines, and dotted lines, respectively). For a weak electric field E_{0} ≪ c B_{0} the particle trajectory follows a spiral motion similar to the previous case of a pure magnetic field until it almost rests. At later times the electric field starts to accelerate it quickly to ultrarelativistic speeds on a timescale 1/ω_{E} ≫ 1/ω_{B}. The particle performs many orbits before being deflected along the parallel direction (e_{z}). Increasing E_{0} decreases this timescale and the particle follows the common E and B direction before performing many gyrations. In the opposite limit of a strong electric field E_{0} ≫ c B_{0} electric acceleration quickly sets in. Figure 11 shows some results of numerical simulations with a weak electric field log(E_{0}/c B_{0}) = − 2, pertinent for almost forcefree neutron star magnetospheres, and initial Lorentz factors log γ_{0} = {0, 4, 8} and log b = {0, −5, −10, −15}.
Fig. 10. Analytic evolution of the Lorentz factor with initial condition log γ_{0} = 8, different damping factors log b = {0, −5, −10, −15}, and different electric field strengths E_{0} relative to B_{0} such that log(E_{0}/c B_{0}) = {−2, 0, 2} in solid lines, dashed lines, and dotted lines, respectively. 
Fig. 11. Evolution of the Lorentz factor for a parallel electromagnetic field configuration with initial Lorentz factors log γ_{0} = {0, 4, 8}, an electric field strength log(E_{0}/c B_{0}) = − 2, and log b = {0, −5, −10, −15}. The black solid lines correspond to the exact analytical solutions. 
Whenever there is a magnetic field aligned electric field component, at late times particles are always accelerated along the common direction. The timescale required is estimated by reckoning the proper time at which the parallel fourvelocity component becomes comparable to the perpendicular fourvelocity component for an initial velocity perpendicular to the magnetic field line. This represents the worst case, useful for comparison to the radiation reaction limit regime.
3.5. Almost cross field
In a plasmafilled magnetosphere, the electric field is efficiently screened meaning that the parallel component of the electric field is negligible with respect to its perpendicular component, E_{∥} ≪ E_{⊥}. As an application towards this configuration, we computed the motion of a particle in an almost cross electric and magnetic field with log(E_{∥} ≪ E_{⊥}) = {−1, −2, −3, −4} and different ratios E_{⊥}/c B = {0.1, 0.999}. The particle starts with an initial velocity in a plane perpendicular to B and Lorentz factor log γ_{0} = 4 in a field with log b = −5. Figure 12 shows the evolution to alignment of the velocity vector with the radiation reaction limit direction for a weak parallel electric field component (see the legend). We note that the angle θ should not be interpreted as the angle between the velocity vector and the magnetic field direction because the velocity in Eq. (17) is not along B. We observe that the time required for alignment is insensitive to the ratio E_{⊥}/c B, but depends strongly on the ratio E_{∥}/E_{⊥}. As expected, a weak parallel component tends to align the trajectory more slowly compared to a strong parallel component. If this alignment occurs on a length scale smaller than the magnetic field curvature radius, the radiation reaction limit could be used without significant loss of accuracy.
Fig. 12. Alignment of the velocity vector with the radiation reaction limit direction given by θ = 0° for different strengths of the parallel electric field component E_{∥} compared to the perpendicular component E_{⊥} for E_{⊥}/cB = 0.1 (dashed lines) and E_{⊥}/cB = 0.999 (solid lines). The damping parameter is log b = −5 and the initial Lorentz factor log γ_{0} = 4. 
3.6. Dipole magnetic field
As a step towards realistic configurations, we also investigate particle motion in a static magnetic dipole. Unfortunately, no simple exact analytical expressions are available for checking the algorithm; therefore, no quantitative accurate convergance test can be performed. First, we study the magnetic drift in the equatorial plane of a dipole field. Next, we look at trapped particles due to the mirror effect.
3.6.1. Magnetic drift
The magnetic drift motion in the equatorial plane of a dipole field is an interesting example to test our algorithm. The characteristic frequency is again ω_{B} and the particle initial Lorentz factor is γ_{0} = 10^{4}. The damping parameter is log b = {0, −5, −10, −15}.
Figure 13 shows the time evolution of the Lorentz factor that tends asymptotically to unity meaning the particle will rest. The particle returns to rest after a typical time controlled by b and shown as coloured vertical lines.
Fig. 13. Decrease in the Lorentz factor due to radiation reaction when drifting in a dipole magnetic field with initial Lorentz factor log γ_{0} = 4 and different damping constants. 
Figure 14 highlights the corresponding particle trajectory in the equatorial plane. For the weakest damping, the motion remains circular for the guiding centre. For the strongest damping, in green and red, the particle suffers from drastic radiative friction and tends to rest on a very short timescale compared to the drifting motion and orbital motion.
Fig. 14. Particle trajectory in the equatorial plane of a dipole magnetic field and associated with Fig. 13. The inset shows the strong damped motion in green and even stronger damping in red where the spiralling is not seen. 
3.6.2. Magnetic mirror
Due to the mirror effect, particles remain trapped in the dipole magnetic field of a star, as they do around the Earth in the Van Allen belt. However, when some dissipation occurs, for instance through the radiation reaction for high damping parameters, the particles quickly crash onto the surface of the magnetic object.
Figure 15 shows the evolution of the Lorentz factor for particles moving in the magnetic dipole field. For weak damping parameters log(τ_{m} ω_{B})≲ − 10, the radiation reaction remains negligible and the particle motion is almost adiabatic with the three characteristic periodic motions: gyration around the magnetic field, bouncing between the north and south magnetic pole, and precession in the azimuthal direction (blue and orange lines in Fig. 16). For log(τ_{m} ω_{B}) = − 5, the cyclotron motion is rapidly damped and the particle falls onto the star (green trajectory). For log(τ_{m} ω_{B}) = 0, the damping is even faster and the particle crashes onto the stellar surface, following a trajectory similar to the previous case (red solid line in Fig. 16).
Fig. 15. Decrease in the Lorentz factor of particles subject to the mirror effect. 
Fig. 16. Particle trajectories in the dipole magnetic field and associated with Fig. 15. For small damping parameters (orange and blue lines) the particle is trapped for a long time in the dipole, whereas for larger damping it quickly crashes onto the stellar surface. 
4. Application to neutron stars
After checking and testing our new algorithm, we are ready to apply it to realistic extreme cases of rotating magnetised neutron stars. The neutron star radius is fixed to R_{*} = 12 km. The accurate configuration of the electromagnetic field is taken from a rotating magnetic dipole in a vacuum and given by Deutsch (1955).
4.1. Relevant parameters without dimension
As a typical frequency we choose the stellar angular frequency ω = Ω_{*} and consider three populations of neutron stars: young pulsars with period P_{*} = 1 s and surface magnetic field strengths B_{*} = 10^{8} T; millisecond pulsars with period P_{*} = 5 ms and B_{*} = 10^{5} T; and magnetars with period P_{*} = 10 s and B_{*} = 10^{10} T. These quantities and their associated normalised strengths and damping parameters a_{B}, a_{E}, and b for electrons, protons, and iron are summarised in Table 1. The normalisation frequency is arbitrary, but from a microscopic point of view the most relevant frequencies are related to the electromagnetic tensor eigenfrequencies. Therefore, the low value of b should not be misinterpreted as a weak feedback of the radiation reaction. It is an artefact of the chosen typical frequency associated with the stellar rotation, which is many orders of magnitude smaller than the cyclotron frequency.
Typical period and surface magnetic field strength of millisecond pulsars, young pulsars, and magnetars.
We distinguish three kind of particles: a first group crashing onto the stellar surface, a second group of trapped particles, and a third group of escaping particles, all accelerated to high energies. Particles are considered trapped when they still have not crashed onto the surface or have not yet escaped the light cylinder. Particles are placed regularly within the light cylinder, starting at rest or with an initial velocity vector oriented along the magnetic field line, directed towards the star or towards infinity, with a Lorentz factor equal to γ_{0} = 10^{3} or starting at rest. The neutron star obliquity is set to χ = {0° ,30° ,60° ,90° ,120° ,150° ,180° }. We found that the final results are not very sensitive to the initial Lorentz factor because charges are immediately accelerated in the direction of the electric field and therefore lose memory about their initial state. Our simulation results are thus summarised for particles starting at rest only. We simulated a total number of 48 particles for each neutron star type and each obliquity, spread around three radii r_{0}: at the surface R_{*}; approximately halfway between the surface and the lightcylinder (a geometric average); and at the light cylinder, thus . For comparison, we performed simulations with and without the radiation reaction.
4.2. Orders of magnitude
Before presenting the accurate numerical simulations of particle trajectories and their radiation reaction in the vicinity of neutron stars, we recall the orders of magnitude of the maximum Lorentz factors expected when charges are accelerated in the electric potential produced by a rotating magnetised perfectly conducting star. The most optimistic view adopts the full potential drop between the pole and the equator as an estimate of the accelerating field, and thus
where r_{B} = c/ω_{B} is the nonrelativistic Larmor radius. If the accelerating potential is only available across the polar caps as expected from nearly forcefree magnetosphere models, the maximum energy corresponds to
which is a factor R_{*}/r_{L} smaller than for the former case. Table 2 summarises the maximum Lorentz factors for electrons, protons, and iron around millisecond pulsars, young pulsars, and magnetars. The values reported in this table for are at best upper limits for the vacuum case. Only an accurate numerical integration of the equation of motion gives robust results, as we now show.
Maximum Lorentz factor orders of magnitude from conservative arguments about neutron star magnetospheres.
4.3. Escaping particles
Particles reaching distances larger than 10 r_{L} are thought to be leaving the neutron star magnetosphere. The run halts when the particle reaches larger distances. Figure 17 shows the histogram of Lorentz factors for electrons, protons, and iron nuclei, irrespective of the magnetic field inclination angle χ. The left column corresponds to a motion with radiation reaction, and the right column to motion without radiation reaction. First, electrons are the most effectively accelerated particles reaching final Lorentz factors up to γ_{f} ∼ 10^{9} in the LLR approximation for millisecond pulsars. This is, however, two orders of magnitude less than without the radiation reaction where γ_{f} ∼ 10^{11}. Second, as expected, protons and iron nuclei acquire much less energy, only about γ_{f} ∼ 10^{6} for millisecond pulsars, whether LLR is used or not. For young pulsars, electrons also reach γ_{f} ∼ 10^{9} in the LLR regime instead of γ_{f} ∼ 10^{11} for the pure Lorentz force. Protons and iron nuclei are much less subject to the radiation reaction, showing no impact on the maximum Lorentz factor remaining at γ_{f} ∼ 10^{4} − 10^{4.5}. For magnetars, the radiation reaction remains negligible irrespective of the nature of each species. Electrons reach energies up to γ_{f} ∼ 10^{7.5}, whereas protons and iron nuclei γ_{f} ∼ 10^{3} − 10^{3.5}. Therefore, the radiation reaction does not significantly perturb the trajectories of particles with lower charge over mass ratios q/m. Contrary to electrons, protons, and iron do not suffer from radiation friction.
Fig. 17. Histogram of escaped particles for millisecond pulsars (top), for young pulsars (middle) and for magnetars (bottom). Electron Lorentz factors are shown in green, protons in red, and iron nuclei in blue. The left column includes the radiation reaction (RR); the right panel does not. 
4.4. Crashed particles
Closer to the star most species quickly crash onto the surface in a time much shorter than the neutron star spin period. Particles crashing onto the neutron star surface are easily recognised by the fact that their final position lies inside the star. Compared to escaping particles, the situation is now reversed; magnetars offer the highest energetic particles heating the surface and millisecond pulsars the lowest energetic particles (see left column of Fig. 18). This is accounted for by the lower surface magnetic field of millisecond pulsars, being three to five orders of magnitude lower than young pulsars or magnetars, respectively. Neglecting the radiation reaction, electrons are able to reach Lorentz factors up to γ_{f} ∼ 10^{11} for magnetars, but only γ_{f} ∼ 10^{8.5} for millisecond pulsars. The radiation reaction impact is strongest for magnetars. However, protons and iron are not perturbed by the radiation reaction except sensibly for magnetars. Nevertheless, we observe that with the radiation reaction protons remain the most energetic particles with final Lorentz factors of about γ_{f} ∼ 10^{7.5} − 10^{8.5} irrespective of the neutron star nature, millisecond, young, or magnetar. For electrons the situation is drastically different. They radiate copiously, decreasing the Lorentz factor by three orders of magnitude compared to the no radiation reaction case in the magnetar environment. The decrease is less pronounced for young or millisecond pulsars, but still perceptible.
4.5. Trapped particles
By default we assume that trapped particles are those not crashing onto the neutron star and not escaping to large distances outside the light cylinder within the simulation time span corresponding to several neutron star periods. Figure 19 summarises the distribution of Lorentz factors for electrons, protons, and iron in the LLR approximation and without radiation reaction. Protons and iron are still insensitive to the radiation reaction except for magnetars. Electrons are much more sensitive to the radiation reaction, decreasing their Lorentz factor by four orders of magnitude for millisecond pulsars, young pulsars, and magnetars. Millisecond pulsars produce trapped protons and iron with energies about γ_{f} ∼ 10^{7}, whereas young pulsars and magnetars one decade more up to γ_{f} ∼ 10^{8}, whether the radiation reaction is included or not. Electrons are trapped with similar Lorentz factors, although slightly less for millisecond pulsars.
4.6. Maximum Lorentz factor
For escaping particles in the wave zone, the gain in energy is limited by the spherical nature of the electromagnetic field, meaning decreasing in strength with distance as 1/r. For a null electromagnetic field Pétri (2021) showed that this severely limits the maximum Lorentz factor to values of
where a_{BL} is the strength parameter measured at the light cylinder. The radiation reaction also remains negligible in this wave zone. Table 3 summarises the relevant parameters at the light cylinder for the three kinds of neutrons stars. Generally speaking, we found no particles with Lorentz factors exceeding γ_{f} ≈ 10^{9.1} in the LLR regime. Because the electromagnetic vacuum used in our simulations corresponds to the one producing the strongest parallel electric field (with respect to the magnetic field), no particle should be created and moving with Lorentz factor higher than 10^{9.1} within the magnetosphere. Equation (35) is satisfied for nonnull electromagnetic waves like those launched by a rotating magnetic dipole. Instead, we found a simple linear relation between the strength parameter a_{BL} and the Lorentz factor such that
Maximum Lorentz factors γ_{max} for the three kinds of particles: electrons / protons / iron nuclei.
This increase in the acceleration efficiency is imputed to the presence of a still strong radial component of the electric field, which was absent in the study of Pétri (2021). The simulations performed in this section only followed a small number of particles, due to the stringent computation time required to accurately evolve the particle velocity and position. Describing the plasma feedback onto the electromagnetic field would require a much larger number of particles coupled to the evolution of the electromagnetic field via Maxwell equations, leading to a PIC code. To date, although PIC codes exist and have been adapted to simulate neutron star magnetospheres, none has yet been able to handle the parameter space explored in the present work. Therefore, let us contrast our results in the light of the existing kinetic descriptions of the magnetosphere.
4.7. Comparison to previous works
Several investigations of particle acceleration and radiation reaction around neutron stars have been attempted in the literature. However, due to severe numerical limitations, studies employing realistic field strengths for the neutron star are very rare. We mention, however, the pioneering work of Finkbeiner et al. (1989) and Finkbeiner et al. (1990), who employed a single test particle approach with the radiation reaction and found acceleration around the neutron star up to Lorentz factors of about γ_{f} ≲ 10^{9} for the Crab parameters. In a similar manner, at very large distances, in the wind zone, Michel & Li (1999) studied particle acceleration without the radiation reaction and found asymptotic values of γ_{f} ≲ 10^{9}. The flaw in these studies is that particles evolve in a prescribed external field without possible feedback due to the plasma around the neutron star. A fully kinetic description of the plasma and field started only recently using PIC schemes; earlier simulations did not take into account the radiation reaction. Unfortunately, the flaw in this approach is the use of unrealistically low field strengths. We mention some of these works here.
Cerutti et al. (2015) studied acceleration for an axisymmetric magnetosphere without the radiation reaction. They obtained a maximum energy for leptons γ_{f} ≲ 10^{3} related linearly to the magnetisation parameter in the plasma. Later, Cerutti et al. (2016) included the radiation reaction force and obtained maximum energies tat were one order of magnitude lower with γ_{f} ≲ 10^{2}. Dissipation in the striped wind due to magnetic reconnection led Cerutti et al. (2020) to the same conclusion. Other PIC simulations performed by Kalapotharakos et al. (2018) using similar algorithms with radiation reaction found similar results with γ_{f} ≲ 10^{3} for pairs. Nevertheless, these authors extrapolated to realistic energies by using rescaling techniques for field strengths, timescales, and space scales. How effectively and consistently this rescaling operates is not clear as the problem is highly nonlinear in a significant radiation reaction regime. General relativity does not significantly change these conclusions, as shown by Philippov & Spitkovsky (2018), who included framedragging effects and found γ_{f} ≲ 500 by extending their special relativistic results in Philippov et al. (2015).
When focusing on the near field of a dipole, Ferrari & Trussoni (1974) found an asymptotic Lorentz factor for electrons of about 10^{8} and slightly larger for protons (almost 10^{9}), but for faster rotations in a field of an oblique rotating dipole with strengths of 5 × 10^{6} T. When the radiation reaction remains irrelevant, their results agree with those of Kulsrud (1972), demonstrating a linear growth with the field strength parameter. Laue & Thielheim (1986) investigates the special case of an orthogonal rotator with radiation reaction, and for typical neutron star parameters they found a maximum energy for electrons of about 10^{9} and for protons of about 10^{6}.
Hadron acceleration has been much less discussed in this context, but it is equally important to understanding the origin of ultrahighenergy cosmic rays. To this end, Guépin et al. (2020) investigated proton and pair acceleration in an aligned neutron star magnetosphere with the radiation reaction. They drastically reduced the neutron star radius and the proton to electron mass ratio for computational purposes, and found the highest energies for pairs of about γ_{f} ≲ 700 and for protons of about γ_{f} ≲ 40. These stateofthe art results emphasise the difficulty in moving towards a realistic and selfconsistent description of neutron star magnetospheres. The main bottleneck is the particle pusher, which requires us to temporally resolve the gyromotion. This drawback is circumvented by employing an approximation called the radiation reaction limit, summarised in Sect. 2. It is therefore important to assess quantitatively the accuracy and efficiency of this alternative approach, which we do in the next section.
5. Comparison with the radiation reaction limit
The results obtained in the previous section rely on the numerical integration of the LLR equation accounting for realistic parameters introducing a huge gap between the gyrofrequency and the neutron star rotation period. The question arises then of how to improve our algorithm or to speed up the computation by several decades. To this end, in this section we compare the LLR results to the radiation reaction limit regime to assert the usefulness of the latter.
Integrating the exact LLR equations requires resolving the gyrofrequency, which is very stringent and impossible to use for a large sample of particles, which is required to perform kinetic simulations such as those done in the PIC or Vlasov codes. We therefore checked the accuracy of the much faster radiation reaction limit approximation where the particle velocity is expressed in terms of the local electromagnetic field (Eq. (17)). To this end, we computed trajectories for electrons and protons in the field of a millisecond pulsar for different magnetic moment inclination angles and different initial particle positions. Because by construction the speed in Eq. (17) is equal to the speed of light v_{±} = c, in the LLR approach particles are kicked with high initial Lorentz factors γ_{0} = 10^{3} and a velocity parallel to v_{±} in order to have comparable initial conditions for both sets of runs.
Figure 20 shows a sample of electron trajectories and demonstrates the reasonable results obtained by this asymptotic regime for a millisecond pulsar. However, the precision depends on the particle initial position. For motions starting at the surface (upper row in Fig. 20) some trajectories are well reproduced by the radiation reaction limit regime. The accuracy is less good for the brown and green paths, although the general trend is conserved. When starting at larger distances from the surface, for example at (middle row of Fig. 20), we observe better agreement between the two regimes. The best results are obtained for particles well away from the surface, starting at r = r_{L} (bottom row of Fig. 20). All trajectories computed in the radiation reaction limit regime overlap with the LLR integration.
Fig. 20. Sample of trajectories for electrons obtained with LLR in the field of a millisecond pulsar (solid lines) and in the radiation reaction limit (dotted symbols). The trajectories are projected along the xy plane in the left panel and on the xz plane in the right panel. Particles are launched from the stellar surface in the top row, at a distance in the middle row, and at a distance r = r_{L} in the bottom row. 
A comparison of trajectories for protons is shown in Fig. 21. Here the agreement is satisfactory within the light cylinder, close to the surface (top row) and at intermediate distances (middle row). For protons starting at the light cylinder radius r = r_{L} the results are more contrasted; some trajectories are well reproduced (in blue, yellow, and orange), and some are false (the brown and green motions) and expected to crash on the surface, but escaping in the radiation reaction limit regime. Iron shows trajectories that are very similar to protons because of the nearly identical mass to charge ratio q/m; the figures are therefore not shown in this almost identical case.
The radiation reaction limit regime is less accurate than the exact LLR integration scheme, but this is partially compensated for by the drastic decrease in computational time, lowered by several orders of magnitude in this approximation. Expression (17) could certainly be improved by carefully investigating these problematic cases, but we do not pursue this aim in this work. We demonstrated however that the velocity in Eq. (17) offers a valuable compromise between a timeconsuming full integration of the equation of motion in LLR and an artificial and unrealistic downscaling of the major physical parameters that make a neutron star a neutron star.
A convergence analysis of the radiation reaction limit integration scheme is shown in Fig. 22 for the relative error. We simulated a sample of 12 particles starting at different locations within the magnetosphere, and compared their last position to a reference solution. As no exact analytical solutions are known, we use as the reference numerical solution the one with the smallest time step. The integration scheme oscillates between the first and second order in time depending on the initial position of the particle, shown by the number r_{0}/r_{L} in the legend. For reference the Δt^{2} behaviour is shown as blue filled circles. In the previous simulations we fixed the time step to Δt ≈ 10^{−4}, and thus we expect a precision better than three digits in all cases. The discrepancy between the Landau–Lifshitz and radiation reaction regimes can therefore not be explained by a discretisation effect. We also checked that the initial condition of the velocity does not impact the trajectory in Landau–Lifshitz. The explanation must be searched for in the deficiency of the radiation reaction regime to satisfactorily account for all possible trajectories. This regime assumes a radiative friction force opposite to the threevelocity vector. However, the 3D version of the LLR equation also contains components along E ∧ B and along E and B when the linear term in velocity is retained. The main discrepancy arises from neglecting this linear term. In order to prove this argument, we designed a simplified version of the Landau–Lifshitz equation by only retaining in the new analytical solution the part of the radiation reaction force directed along the velocity (which is valid for ultrarelativistic speeds). The expressions for the fourvelocity then become
Fig. 22. Convergence of the radiation reaction limit regime showing the method of integration to be between the first and second order in time, depending on the initial position of the particle given by the number r_{0}/r_{L} in the legend. The Δt^{2} decrease in shown as blue filled circles. 
replacing Eq. (12). The results are shown in Fig. 23 for the exact Landau–Lifshitz equation in solid lines, the approximated Landau–Lifshitz equation in dashed lines, and the radiation reaction regime in dots. We observe some significant differences, notably in the middle right panel. We note however that the radiation reaction regime gives accurate results at low computational time expense for the majority of cases.
Fig. 23. Comparison between the radiation reaction limit (dotted symbols), the LLR (thick solid lines), and the approximated LLR motion (thin dashed lines). 
6. Conclusions
Strongly magnetised rotating neutron stars are powerful and efficient particle accelerators able to accelerate leptons and hadrons to Lorentz factors as high as 10^{9} for the former and slightly less for the latter. This upper limit remains largely independent of the nature of neutron star, millisecond pulsar, young pulsar, or magnetar. We achieved these results by implementing realistic parameters in our particle pusher based on the exact solution of the LLR approximation of the equation of motion. Through extensive numerical tests, we show that our scheme is second order in proper time.
The simulation results are accurate and robust, but at the expense of high computational cost because of the need to resolve the gyromotion, which is many decades smaller than the neutron star spin period. The radiation reaction limit regime offers a good compromise between accuracy and computational cost, but the simplistic expression used is unable to reproduce all trajectories satisfactorily. Nevertheless, it could be conceivable to improve this expression by taking into account a finite Lorentz factor and special electromagnetic field configuration when the radiation reaction is negligible due to a weak accelerating electric field. Nevertheless this extension is left for future work.
A straightforward implementation of the above pusher into a PIC code or Vlasov codes is prevented by the fact that LLR uses the proper time as integration parameter. However, its conversion into an inertial observer time is feasible, as shown by Pétri (2020). The next logical step would then be to shift from the test particle motion to a fully kinetic plasma simulation where the particle charge and current densities retroact to the electromagnetic field via Maxwell equations.
Acknowledgments
I am grateful to the referee for helpful comments and suggestions that helped to improved this work. This work has been supported by the CEFIPRA grant IFC/F5904B/2018 and ANR20CE310010. The authors would like to acknowledge the High Performance Computing Center of the University of Strasbourg for supporting this work by providing scientific support and access to computing resources. Part of the computing resources were funded by the Equipex Equip@Meso project (Programme Investissements d’Avenir) and the CPER Alsacalcul/Big Data.
References
 Abraham, M. 1902, Ann. Phys., 315, 105 [NASA ADS] [CrossRef] [Google Scholar]
 Abraham, M. 1904, Ann. Phys., 319, 236 [NASA ADS] [CrossRef] [Google Scholar]
 Boghosian, B. M. 1987, PhD Thesis, University of California [Google Scholar]
 Boris, J. 1970, Proceeding of Fourth Conference on Numerical Simulations of Plasmas [Google Scholar]
 Cerutti, B., Philippov, A., Parfrey, K., & Spitkovsky, A. 2015, MNRAS, 448, 606 [Google Scholar]
 Cerutti, B., Philippov, A. A., & Spitkovsky, A. 2016, MNRAS, 457, 2401 [Google Scholar]
 Cerutti, B., Philippov, A. A., & Dubus, G. 2020, A&A, 642, A204 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Deutsch, A. J. 1955, Ann. d’Astrophys., 18, 1 [Google Scholar]
 Dirac, P. A. M. 1938, Proc. R. Soc. Lond. Ser. A, 167, 148 [NASA ADS] [CrossRef] [Google Scholar]
 Elkina, N. V., Fedotov, A. M., Herzing, C., & Ruhl, H. 2014, Phys. Rev. E, 89, 053315 [NASA ADS] [CrossRef] [Google Scholar]
 Ferrari, A., & Trussoni, E. 1974, A&A, 36, 267 [NASA ADS] [Google Scholar]
 Finkbeiner, B., Herold, H., Ertl, T., & Ruder, H. 1989, A&A, 225, 479 [NASA ADS] [Google Scholar]
 Finkbeiner, B., Herold, H., & Ruder, H. 1990, A&A, 238, 462 [NASA ADS] [Google Scholar]
 Fulton, T., & Rohrlich, F. 1960, Ann. Phys., 9, 499 [NASA ADS] [CrossRef] [Google Scholar]
 Gordon, D. F., & Hafizi, B. 2021, Comput. Phys. Commun., 258, 107628 [NASA ADS] [CrossRef] [Google Scholar]
 Gordon, D. F., Hafizi, B., & Palastro, J. 2017a, AIP Conf. Proc., 1812, 050002 [CrossRef] [Google Scholar]
 Gordon, D. F., Palastro, J. P., & Hafizi, B. 2017b, Phys. Rev. A, 95, 033403 [NASA ADS] [CrossRef] [Google Scholar]
 Guépin, C., Cerutti, B., & Kotera, K. 2020, A&A, 635, A138 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hadad, Y., Labun, L., Rafelski, J., et al. 2010, Phys. Rev. D, 82, 096012 [NASA ADS] [CrossRef] [Google Scholar]
 Heintzmann, H., & Schrüfer, E. 1973, Phys. Lett. A, 43, 287 [NASA ADS] [CrossRef] [Google Scholar]
 Herold, H., Ertl, T., & Ruder, H. 1985, Mitteilungen der Astronomischen Gesellschaft Hamburg, 63, 174 [NASA ADS] [Google Scholar]
 Kalapotharakos, C., Brambilla, G., Timokhin, A., Harding, A. K., & Kazanas, D. 2018, ApJ, 857, 44 [NASA ADS] [CrossRef] [Google Scholar]
 Kelner, S. R., Prosekin, A. Y., & Aharonian, F. A. 2015, AJ, 149, 33 [Google Scholar]
 Kulsrud, R. M. 1972, ApJ, 174, L25 [NASA ADS] [CrossRef] [Google Scholar]
 Landau, L., & Lifchitz, E. 1989, Physique théorique: Tome 2, Théorie des champs (Moscou: Mir) [Google Scholar]
 Laue, H., & Thielheim, K. O. 1986, ApJS, 61, 465 [NASA ADS] [CrossRef] [Google Scholar]
 Li, F., Decyk, V. K., Miller, K. G., et al. 2021, J. Comput. Phys., 438, 110367 [NASA ADS] [CrossRef] [Google Scholar]
 Lorentz, H. A. H. A. 1916, The Theory of Electrons and its Applications to the Phenomena of Light and Radiant Heat (Leipzig: G.E. Stechert) [Google Scholar]
 Mestel, L., Robertson, J. A., Wang, Y. M., & Westfold, K. C. 1985, MNRAS, 217, 443 [NASA ADS] [CrossRef] [Google Scholar]
 Michel, F., & Li, H. 1999, Phys. Rep., 318, 227 [NASA ADS] [CrossRef] [Google Scholar]
 Olver, F. W. J. 2010, NIST Handbook of Mathematical Functions (Cambridge: Cambridge University Press) [Google Scholar]
 Philippov, A. A., & Spitkovsky, A. 2018, ApJ, 855, 94 [NASA ADS] [CrossRef] [Google Scholar]
 Philippov, A. A., Spitkovsky, A., & Cerutti, B. 2015, ApJ, 801, L19 [NASA ADS] [CrossRef] [Google Scholar]
 Piazza, A. D. 2008, Lett. Math. Phys., 83, 305 [NASA ADS] [CrossRef] [Google Scholar]
 Pétri, J. 2020, J. Plasma Phys., 86, 825860402 [CrossRef] [Google Scholar]
 Pétri, J. 2021, MNRAS, 503, 2123 [CrossRef] [Google Scholar]
 Rohrlich, F. 2007, Classical Charged Particles, 3rd edn. (Singapore: World Scientific Pub Co Inc) [CrossRef] [Google Scholar]
 Tomczak, I., & Pétri, J. 2020, J. Plasma Phys., 86, 825860401 [NASA ADS] [CrossRef] [Google Scholar]
 Vranic, M., Martins, J. L., Fonseca, R. A., & Silva, L. O. 2016, Comput. Phys. Commun., 204, 141 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
Typical period and surface magnetic field strength of millisecond pulsars, young pulsars, and magnetars.
Maximum Lorentz factor orders of magnitude from conservative arguments about neutron star magnetospheres.
Maximum Lorentz factors γ_{max} for the three kinds of particles: electrons / protons / iron nuclei.
All Figures
Fig. 1. Increase in the Lorentz factor due to radiation reaction for different initial Lorentz factors log γ_{0} = {0, 4, 8} and different damping factors log b = {0, −5, −10, −15}. The vertical lines show the time when the damping sets in before the electric acceleration phase starts. The coloured points show the simulation results and the black solid lines correspond to the exact analytical solutions. 

In the text 
Fig. 2. Relative error of the position y, z and time t shown in the legend. The error decreases with second order in Δτ as given by the green line Δτ^{2} for log γ_{0} = 4 and log b = −5. The t and z errors overlap and are indistinguishable. 

In the text 
Fig. 3. Particle orbit in a uniform and constant magnetic field and subject to radiation reaction. The initial Lorentz factor is log γ_{0} = 4. The inset shows the strong damped motion in green and even stronger damping in red where the spiralling is not seen. 

In the text 
Fig. 4. Decrease in the Lorentz factor due to radiation reaction in a uniform and constant magnetic field associated with the orbits shown in Fig. 3. The coloured points show the simulation results and the black solid lines correspond to the exact analytical solutions. 

In the text 
Fig. 5. Comparison between the analytical solution, Eq. (31) (red solid line) and the numerical simulation (blue dots) for log γ_{0} = 4 and log b = −5. 

In the text 
Fig. 6. Relative error of the position x and y as shown (see inset for legend). The error decreases with second order in Δτ, as given by the violet line Δτ^{2} for log γ_{0} = 4 and log b = { − 5, −10}. 

In the text 
Fig. 7. Decrease in the Lorentz factor due to radiation reaction in a cross electric and magnetic field. The coloured points show the simulation results and the black solid lines correspond to the exact analytical solutions. 

In the text 
Fig. 8. Orbit in the electric drift frame (x′,y′) for log γ_{0} = 4 and log b = −5 for different proper time steps log(ω_{B} Δτ) = {−1, −2}. 

In the text 
Fig. 9. Relative error of the particle position associated with the electric drift motion shown in Fig. 8. 

In the text 
Fig. 10. Analytic evolution of the Lorentz factor with initial condition log γ_{0} = 8, different damping factors log b = {0, −5, −10, −15}, and different electric field strengths E_{0} relative to B_{0} such that log(E_{0}/c B_{0}) = {−2, 0, 2} in solid lines, dashed lines, and dotted lines, respectively. 

In the text 
Fig. 11. Evolution of the Lorentz factor for a parallel electromagnetic field configuration with initial Lorentz factors log γ_{0} = {0, 4, 8}, an electric field strength log(E_{0}/c B_{0}) = − 2, and log b = {0, −5, −10, −15}. The black solid lines correspond to the exact analytical solutions. 

In the text 
Fig. 12. Alignment of the velocity vector with the radiation reaction limit direction given by θ = 0° for different strengths of the parallel electric field component E_{∥} compared to the perpendicular component E_{⊥} for E_{⊥}/cB = 0.1 (dashed lines) and E_{⊥}/cB = 0.999 (solid lines). The damping parameter is log b = −5 and the initial Lorentz factor log γ_{0} = 4. 

In the text 
Fig. 13. Decrease in the Lorentz factor due to radiation reaction when drifting in a dipole magnetic field with initial Lorentz factor log γ_{0} = 4 and different damping constants. 

In the text 
Fig. 14. Particle trajectory in the equatorial plane of a dipole magnetic field and associated with Fig. 13. The inset shows the strong damped motion in green and even stronger damping in red where the spiralling is not seen. 

In the text 
Fig. 15. Decrease in the Lorentz factor of particles subject to the mirror effect. 

In the text 
Fig. 16. Particle trajectories in the dipole magnetic field and associated with Fig. 15. For small damping parameters (orange and blue lines) the particle is trapped for a long time in the dipole, whereas for larger damping it quickly crashes onto the stellar surface. 

In the text 
Fig. 17. Histogram of escaped particles for millisecond pulsars (top), for young pulsars (middle) and for magnetars (bottom). Electron Lorentz factors are shown in green, protons in red, and iron nuclei in blue. The left column includes the radiation reaction (RR); the right panel does not. 

In the text 
Fig. 18. Same as Fig. 17, but for crashed particles. 

In the text 
Fig. 19. Same as Fig. 17, but for trapped particles. 

In the text 
Fig. 20. Sample of trajectories for electrons obtained with LLR in the field of a millisecond pulsar (solid lines) and in the radiation reaction limit (dotted symbols). The trajectories are projected along the xy plane in the left panel and on the xz plane in the right panel. Particles are launched from the stellar surface in the top row, at a distance in the middle row, and at a distance r = r_{L} in the bottom row. 

In the text 
Fig. 21. Same as Fig. 20, but for protons. 

In the text 
Fig. 22. Convergence of the radiation reaction limit regime showing the method of integration to be between the first and second order in time, depending on the initial position of the particle given by the number r_{0}/r_{L} in the legend. The Δt^{2} decrease in shown as blue filled circles. 

In the text 
Fig. 23. Comparison between the radiation reaction limit (dotted symbols), the LLR (thick solid lines), and the approximated LLR motion (thin dashed lines). 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.