# Particle acceleration and radiation reaction in a strongly magnetized rotating dipole

Pétri Jérôme

► To cite this version:

Pétri Jérôme. Particle acceleration and radiation reaction in a strongly magnetized rotating dipole. Astronomy & Astrophysics - A&A, 2022, 666, pp.A5. 10.1051/0004-6361/202243634. hal-03722272

HAL Id: hal-03722272

https://hal.science/hal-03722272v1

Submitted on 12 Apr 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Particle acceleration and radiation reaction in a strongly magnetised rotating dipole

J. Pétri

Université de Strasbourg, CNRS, Observatoire astronomique de Strasbourg, UMR 7550, 67000 Strasbourg, France
e-mail: jerome.petri@astro.unistra.fr

Received 25 March 2022 / Accepted 30 June 2022

## ABSTRACT

**Context.** Neutron stars are surrounded by ultra-relativistic particles efficiently accelerated by ultra-strong electromagnetic fields. These particles copiously emit high-energy 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 \times 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 $\gamma_e \approx 10^8-10^9$, whereas protons reach energies up to $\gamma_p \approx 10^5-10^6$ and iron up to $\gamma \approx 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

## 1. Introduction

Neutron stars are known to harbour ultra-strong magnetic fields close to or even above the quantum critical limit of $B_c \approx 4.4 \times 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 $\gamma \approx 10^9$. In such an extreme environment, the radiation reaction is expected to drastically perturb their trajectory compared to the pure Lorentz force motion. High-energy and very high-energy 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 non-linearity 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 ultra-relativistic 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 ultra-relativistic 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 synchro-curvature radiation of ultra-relativistic 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 ultra-strong 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, time-symmetric 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 particle-in-cell (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 low-frequency 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 self-force produced by an accelerated charge is usually described by the Lorentz-Abraham-Dirac (LAD) equation (Abraham 1902, 1904; Lorentz 1916; Dirac 1938). Unfortunately, this self-force 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

$$ \frac{du^i}{d\tau} = \frac{q}{m} F^{ik} u_k + \frac{q \tau_m}{m} g^i, \quad (1a) $$

$$ g^i = \partial_\ell F^{ik} u_k u^\ell + \frac{q}{m} \left( F^{ik} F_{k\ell} u^\ell + (F^{\ell m} u_m) (F_{\ell k} u^k) \frac{u^i}{c^2} \right), \quad (1b) $$

where $q$ and $m$ are the particle charge and rest mass, $u^i$ its four-velocity, $\tau$ its proper time, $F^{ik}$ the electromagnetic or Faraday tensor, $c$ the speed of light, and $\tau_m$ the light crossing time across the particle classical radius $r_m$ (within a factor unity)

$$ \tau_m = \frac{q^2}{6 \pi \epsilon_0 m c^3}. \quad (2) $$

It is advantageous to express it in terms of the electron classical radius $r_e$ crossing time amounting to

$$ \tau_e = \frac{2}{3} \frac{r_e}{c} = 6.26 \times 10^{-24} \text{ s}. \quad (3) $$

The typical timescale for the radiation reaction is therefore

$$ \tau_m = \frac{2}{3} \frac{r_m}{c} = \left( \frac{q^2/e^2}{m/m_e} \right) \tau_e. \quad (4) $$

For instance, for protons this time is three orders of magnitude less than for leptons:

$$ \tau_p = \frac{m_e}{m_p} \tau_e = 3.41 \times 10^{-27} \text{ s}. \quad (5) $$

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 $F^{i_k}$. 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

$$ \frac{du}{d\tau} = G u, \quad (6) $$

where the electromagnetic tensor $F$ has been replaced by $G = qF/m$ to absorb the charge over mass ratio, we decompose the four-velocity $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 $G_{i_k}$ are $\pm \lambda_E$, whereas the imaginary eigenvalues are $\pm i\lambda_B$, because $\lambda_E$ and $\lambda_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

$$ G u_E = \pm \lambda_E u_E, \quad (7a) $$

$$ G u_B = \pm i \lambda_B u_B. \quad (7b) $$
The vector components $u_E$ and $u_B$ are obtained by defining the projection operators onto the subspaces $E$ and $B$ by (Boghosian 1987)

$$P = \frac{\lambda_B^2 I + G^2}{\lambda_E^2 + \lambda_B^2}, \qquad (8a)$$

$$Q = \frac{\lambda_E^2 I - G^2}{\lambda_E^2 + \lambda_B^2}, \qquad (8b)$$

where $I$ is the identity matrix. These operators are well defined only if $\lambda_E^2 + \lambda_B^2 \neq 0$. If both electromagnetic invariants vanish, we retrieve a null-like field, which requires a different treatment, as given for instance by Pétri (2021). In the non null-like field we get

$$u_E = P u, \qquad (9a)$$

$$u_B = Q u. \qquad (9b)$$

The equation of motion decouples into two parts given by

$$\frac{d^2 u_E}{d\tau^2} = +\lambda_E^2 u_E, \qquad (10a)$$

$$\frac{d^2 u_B}{d\tau^2} = -\lambda_B^2 u_B. \qquad (10b)$$

The exact analytical solutions with initial conditions $u_E^0 = P u^0$ and $u_B^0 = Q u^0$ are

$$u_E(\tau) = u_E^0 \cosh(\lambda_E \tau) + G u_E^0 \frac{\sinh(\lambda_E \tau)}{\lambda_E}, \qquad (11a)$$

$$u_B(\tau) = u_B^0 \cos(\lambda_B \tau) + G u_B^0 \frac{\sin(\lambda_B \tau)}{\lambda_B}. \qquad (11b)$$

Adding the radiation reaction in the LLR limit leads to the exact expression

$$\frac{u_E(\tau)}{c} = \frac{u_E^0 \cosh(\lambda_E \tau) + G u_E^0 \sinh(\lambda_E \tau) / \lambda_E}{\sqrt{|u_E^0|^2 + |u_B^0|^2 e^{-2\alpha \tau}}}, \qquad (12a)$$

$$\frac{u_B(\tau)}{c} = \frac{u_B^0 \cos(\lambda_B \tau) + G u_B^0 \sin(\lambda_B \tau) / \lambda_B}{\sqrt{|u_B^0|^2 + |u_E^0|^2 e^{2\alpha \tau}}}, \qquad (12b)$$

with $\alpha = \tau_m (\lambda_E^2 + \lambda_B^2)$. These expressions are similar to the original formulas found by Heintzmann & Schrüfer (1973). The radiation reaction effect becomes perceptible after a time $\tau \approx 1/\alpha$. The component $u_E$ is associated with the accelerating motion induced by the electric field, whereas the $u_B$ component is related to the gyro-motion in the magnetic field. When $\alpha$ vanishes, the radiation reaction effect disappears. The denominators in $u_E$ and $u_B$ reduce to unity and the solutions to the Lorentz force four-velocity 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 $\omega$ 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/\omega$, the timescale $T_0 = 1/\omega$, the magnetic field strength $B_0 = m\omega/q$, and the electric field strength $E_0 = cB_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 $\omega \tau_m$ according to the following definitions:

$$a_B = \frac{B}{B_0} = \frac{\omega_B}{\omega}, \qquad (13a)$$

$$a_E = \frac{E}{E_0} = \frac{\omega_E}{\omega}, \qquad (13b)$$

$$b = \omega \tau_m. \qquad (13c)$$

Introducing the weighted and normalised electromagnetic field tensor by $\tilde{F}^{ik} = q F^{ik} / m\omega$ and a normalised time $\tilde{\tau} = \omega \tau$, the Landau-Lifshitz Eq. (1) is rewritten without dimensions as

$$\frac{d\tilde{u}^i}{d\tilde{\tau}} = \tilde{F}^{ik} \tilde{u}_k + b \tilde{g}^i, \qquad (14a)$$

$$\tilde{g}^i = \tilde{\partial}_\ell \tilde{F}^{ik} \tilde{u}_k \tilde{u}^\ell + (\tilde{F}^{ik} \tilde{F}_{k\ell} \tilde{u}^\ell + (\tilde{F}^{\ell m} \tilde{u}_m) (\tilde{F}_{\ell k} \tilde{u}^k) \tilde{u}^i). \qquad (14b)$$

The normalised and reduced Landau-Lifshitz equation reads

$$\frac{d\tilde{u}^i}{d\tilde{\tau}} = \tilde{F}^{ik} \tilde{u}_k + b \left( \tilde{F}^{ik} \tilde{F}_{k\ell} \tilde{u}^\ell + (\tilde{F}^{\ell m} \tilde{u}_m) (\tilde{F}_{\ell k} \tilde{u}^k) \tilde{u}^i \right). \qquad (15)$$

The particle four-velocity 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/\alpha$, becomes larger than the characteristic timescale of our system, which is $1/\omega$. Expressed in quantities without dimension, we get $\tau_m \omega (a_E^2 + a_B^2) = b (a_E^2 + a_B^2) \ll 1$. 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 $_2F_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 velocity-Verlet algorithm:

$$\mathbf{u}^{n+1/2} = \mathbf{L}(\Delta\tau/2, \mathbf{u}^n, \mathbf{E}(\mathbf{x}^n), \mathbf{B}(\mathbf{x}^n)), \qquad (16a)$$

$$\mathbf{x}^{n+1} = \mathbf{x}^n + \mathbf{u}^{n+1/2} \Delta\tau, \qquad (16b)$$

$$\mathbf{u}^{n+1} = \mathbf{L}(\Delta\tau/2, \mathbf{u}^{n+1/2}, \mathbf{E}(\mathbf{x}^{n+1}), \mathbf{B}(\mathbf{x}^{n+1})). \qquad (16c)$$

The superscript $n$ refers to the proper time $\tau^n = n \Delta\tau$, and the same for the half-integer superscript $\tau^{n+1/2} = (n + 1/2) \Delta\tau$. 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 velocity-Verlet algorithm is superior.

Before using our code to compute particle acceleration and radiation in the ultra-strong 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 ultra-strong electromagnetic fields, such as those present around neutron stars, radiation reaction plays an important role. In the asymptotic limit of ultra-relativistic 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

$$v_{\pm} = \frac{\mathbf{E} \wedge \mathbf{B} \pm (E_0 \mathbf{E}/c + c B_0 \mathbf{B})}{E_0^2/c^2 + B^2}, \quad (17)$$

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 $I_1 = \mathbf{E}^2 - c^2 \mathbf{B}^2 = E_0^2 - c^2 B_0^2$ and $I_2 = c \mathbf{E} \cdot \mathbf{B} = c E_0 B_0$. Imposing $E_0 \ge 0$ we find

$$E_0^2 = \frac{1}{2} (I_1 + \sqrt{I_1^2 + 4 I_2}), \quad (18a)$$

$$c B_0 = \operatorname{sign}(I_2) \sqrt{E_0^2 - I_1}. \quad (18b)$$

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, $|I_2| \ll |I_1|$ and $I_1 < 0$. The velocity then reduces to

$$v_{\pm} \approx \frac{\mathbf{E} \wedge \mathbf{B}}{B^2} \pm \operatorname{sign}(\mathbf{E} \cdot \mathbf{B}) \frac{\sqrt{c^2 B^2 - E^2}}{B^2} \mathbf{B}. \quad (19)$$

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 gyro-motion being absent in this picture. We note that the velocity component along the magnetic field reverses sign when crossing a point where $\mathbf{E} \cdot \mathbf{B}$ changes sign. These regions are able to trap particles depending on their charge and on the $(\mathbf{E} \cdot \mathbf{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_{\pm} = \pm c \mathbf{B}/B$.

Second, if the electric field is orthogonal to the magnetic field, $\mathbf{E} \cdot \mathbf{B} = 0$ and $E < c \mathbf{B}$, the particle motion is decomposed into an electric drift and a motion along $\mathbf{B}$ such that

$$v_{\pm} = \frac{\mathbf{E} \wedge \mathbf{B}}{B^2} \pm \frac{\sqrt{c^2 B^2 - E^2}}{B^2} \mathbf{B}. \quad (20)$$

This expression holds well within the light-cylinder of a force-free 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 \approx 10^{20}$ and $\gamma \approx 10^{10}$. Our main purpose in this section is to check that the results are not affected by round-off 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 z-axis such that $\mathbf{E} = E \mathbf{e}_z$, we obtain

$$\frac{u^t}{c} = \gamma(\tau) = \frac{\gamma_0 c \cosh(\omega_E \tau) + u_z^0 \sinh(\omega_E \tau)}{\sqrt{\gamma_0^2 c^2 - u_\parallel^2 - u_\perp^2 e^{-2\alpha\tau}}}, \quad (21a)$$

$$\frac{u^x}{c} = \frac{u_x^0}{\sqrt{(\gamma_0^2 c^2 - u_\parallel^2) e^{2\alpha\tau} - u_\perp^2}}, \quad (21b)$$

$$\frac{u^y}{c} = \frac{u_y^0}{\sqrt{(\gamma_0^2 c^2 - u_\parallel^2) e^{2\alpha\tau} - u_\perp^2}}, \quad (21c)$$

$$\frac{u^z}{c} = \frac{u_z^0 \cosh(\omega_E \tau) + \gamma_0 c \sinh(\omega_E \tau)}{\sqrt{\gamma_0^2 c^2 - u_\parallel^2 - u_\perp^2 e^{-2\alpha\tau}}}, \quad (21d)$$

with $u_\parallel = u_z^0$ the initial four-velocity component along $\mathbf{E}$, $u_\perp$ the initial four-velocity component perpendicular to $\mathbf{E}$, and $\alpha = \tau_m \omega_E^2$.

For a particle starting at rest, $u_\parallel = u_\perp = 0$ and $\gamma_0 = 1$, the four-velocity simplifies drastically to

$$u^i = c (\cosh(\omega_E \tau), 0, 0, \sinh(\omega_E \tau)). \quad (22)$$

This four-velocity 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 four-position is given by introducing the two complex functions with the help of the hypergeometric functions ${}_2F_1$ (Olver 2010) such that

$$J_1(\tau) = \frac{e^{(\alpha+\omega_E)\tau}}{\alpha + \omega_E} {}_2F_1\left(\frac{1}{2}, 1 + \frac{1}{2\tau}; \frac{3}{2} + \frac{1}{2\tau}; \frac{\gamma^2 e^{2\alpha\tau}}{\gamma^2 - 1}\right), \quad (23a)$$

$$J_2(\tau) = \frac{e^{(\alpha-\omega_E)\tau}}{\alpha - \omega_E} {}_2F_1\left(\frac{1}{2}, 1 - \frac{1}{2\tau}; \frac{3}{2} - \frac{1}{2\tau}; \frac{\gamma^2 e^{2\alpha\tau}}{\gamma^2 - 1}\right). \quad (23b)$$
Then the time and position are given by

$$t = - \frac{\gamma \sqrt{\gamma^2 (e^{2\alpha\tau} - 1) + 1}}{2 (\gamma^2 - 1)} (J_2(\tau) + J_1(\tau)) + C_0, \quad (24a)$$

$$x/c = 0, \quad (24b)$$

$$y/c = \frac{1}{\alpha} \arctan \left( \sqrt{\frac{\gamma^2 (e^{2\alpha\tau} - 1) + 1}{\gamma^2 - 1}} \right) + C_2, \quad (24c)$$

$$z/c = \frac{\gamma \sqrt{\gamma^2 (e^{2\alpha\tau} - 1) + 1}}{2 (\gamma^2 - 1)} (J_2(\tau) - J_1(\tau)) + C_3, \quad (24d)$$

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 $\tau_{acc} \sim 1/\omega_E$. On the other hand, the radiation damping timescale is $\tau_{rad} \sim 1/\alpha$. The ratio of the two timescales is therefore $\tau_{acc}/\tau_{rad} \sim \tau_m \omega_E = b$. As expected, for small damping parameters $b \ll 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 $\tau_{rad} \sim \tau_{acc}/b \gg \tau_{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 $\omega_E$ as normalisation, which leads to a normalised proper time $\tilde{\tau} = \omega_E \tau$. Therefore, the only relevant parameter apart from the initial conditions is $b = \tau_m \omega_E$ and $\alpha \tau = b \tilde{\tau}$. 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_\parallel = 0$ and $u_\perp \neq 0$. We chose different initial Lorentz factors in the set $\log \gamma_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 $\omega_E \tau \sim 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 $\gamma_0 \gg 1$. The time derivative of the Lorentz factor is always negative for $\gamma_0 > 1$ because

$$\frac{d\gamma(\tau)}{d\tau} = -\alpha \gamma_0 (\gamma_0^2 - 1) < 0, \quad (25)$$

meaning that the particle first decelerates due to the radiative friction. At large times, when $\alpha \tau \gg 1$ and $\omega_E \tau \gg 1$, the Lorentz factor behaves as $\gamma(\tau) \approx \cosh(\omega_E \tau)$, losing its information about the initial state. It resembles the motion of a particle starting at rest, independently of $\gamma_0$. This occurs because the perpendicular motion is strongly damped, $\lim_{\tau \to +\infty} u_\perp \to 0$, and only the parallel velocity $u_\parallel$ survives at large times with $\lim_{\tau \to +\infty} u_\parallel \to c \sinh(\omega_E \tau)$. In between, the normalised time remains small and the Lorentz factor can be approximated by

$$\gamma(\tau) \approx \frac{\gamma_0 \cosh(\omega_E \tau)}{\sqrt{1 + 2 \gamma_0^2 \alpha \tau}}. \quad (26)$$

Therefore, before the acceleration phase starts, there is a deceleration step arising at time $\omega_E \tau \sim 1/2 \gamma_0^2 b$. These values agree with the curves in Fig. 1. If $\gamma_0^2 b \lesssim 1$, 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 $\gamma_0 = 10^4$ and $b = 10^{-10}$ (orange dots) or $\gamma_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 $\gamma_0$ because at time $\omega_E \tau \approx 1$ it braked to a momentum

$$u_\perp(\tau) \approx \frac{u_\perp^0}{\sqrt{1 + 2 \gamma_0^2 b}}. \quad (27)$$

Thus, the radiation reaction again impacts the motion if $\gamma_0^2 b \gtrsim 1$. Consequently, it is the combination $\gamma_0^2 b$ that controls the damping efficiency, not $b$ alone found from the simple arguments above.

Because the four-position 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 $\Delta\tau$ for $\log \gamma_0 = 4$ and $\log b = -5$ (in blue, orange, green, and red, respectively). The second-order 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 velocity-Verlet algorithm shown in Sect. 2.

## 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 four-velocity in a magnetic field directed along the $z$-axis with $\mathbf{B} = B \mathbf{e}_z$ is given by

$$\frac{u'}{c} = \gamma(\tau) = \frac{\gamma_0 c}{\sqrt{\gamma_0^2 c^2 - u_\parallel^2 - u_\perp^2 e^{-2\alpha\tau}}}, \quad (28a)$$

$$\frac{u^x}{c} = \frac{u_x^0 \cos(\omega_B \tau) + u_y^0 \sin(\omega_B \tau)}{\sqrt{(\gamma_0^2 c^2 - u_\parallel^2) e^{2\alpha\tau} - u_\perp^2}}, \quad (28b)$$

$$\frac{u^y}{c} = \frac{u_y^0 \cos(\omega_B \tau) - u_x^0 \sin(\omega_B \tau)}{\sqrt{(\gamma_0^2 c^2 - u_\parallel^2) e^{2\alpha\tau} - u_\perp^2}}, \quad (28c)$$

$$\frac{u^z}{c} = \frac{u_z^0}{\sqrt{\gamma_0^2 c^2 - u_\parallel^2 - u_\perp^2 e^{-2\alpha\tau}}}, \quad (28d)$$

with $u_\parallel = u_z^0$ the initial four-velocity component along $\mathbf{B}$, $u_\perp$ the initial four-velocity component perpendicular to $\mathbf{B}$, and $\alpha = \tau_m \omega_B^2$. Apart from the change in the gyro-frequency, the magnetic field strength impacts only the timescale for the decay via the exponential terms of arguments $2\alpha\tau$.

To perform simulations we use the characteristic frequency $\omega_B$ as normalisation with $\tilde{\tau} = \omega_B \tau$. Therefore, the only relevant parameters, apart from the initial conditions, are $b = \tau_m \omega_B$ and $\alpha \tau = b \tilde{\tau}$. The length scale is therefore given in units of the non-relativistic Larmor radius

$$r_B = \frac{c}{\omega_B}. \quad (29)$$

Integrating the four-velocity vector, an exact analytical expression for the four-position is computed with the help of the
<figure><img src="image_2.png" /></figure>

**Fig. 1.** Increase in the Lorentz factor due to radiation reaction for different initial Lorentz factors $\log \gamma_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.

<figure><img src="image_3.png" /></figure>

**Fig. 2.** Relative error of the position y, z and time t shown in the legend. The error decreases with second order in $\Delta\tau$ as given by the green line $\Delta\tau^2$ for $\log \gamma_0 = 4$ and $\log b = -5$. The t and z errors overlap and are indistinguishable.

hypergeometric functions ${}_2F_1$. Introducing the complex functions

$$H_1(\tau) = e^{+i\omega_B\tau} {}_2F_1\left(\frac{1}{2}, +\frac{i}{2b}; 1 + \frac{i}{2b}; \frac{\gamma^2 e^{-2\alpha\tau}}{\gamma^2 - 1}\right), \quad (30a)$$

$$H_2(\tau) = e^{-i\omega_B\tau} {}_2F_1\left(\frac{1}{2}, -\frac{i}{2b}; 1 - \frac{i}{2b}; \frac{\gamma^2 e^{-2\alpha\tau}}{\gamma^2 - 1}\right). \quad (30b)$$

the solution reads

$$t = \frac{1}{\alpha} \tanh^{-1} \left( \frac{\gamma e^{\alpha \tau}}{\sqrt{\gamma^2 (e^{2\alpha \tau} - 1) + 1}} \right) + C_0, \quad (31a)$$

$$x/r_B = \frac{H_1(\tau) + H_2(\tau)}{2i} + C_1, \quad (31b)$$

$$y/r_B = \frac{H_1(\tau) - H_2(\tau)}{2} + C_2, \quad (31c)$$

$$z/r_B = 0, \quad (31d)$$

where the $C_i$ with $i \in [0...2]$ are complex constants of integration to satisfy the initial conditions.

<figure><img src="image_4.png" /></figure>

**Fig. 3.** Particle orbit in a uniform and constant magnetic field and subject to radiation reaction. The initial Lorentz factor is $\log \gamma_0 = 4$. The inset shows the strong damped motion in green and even stronger damping in red where the spiralling is not seen.

<figure><img src="image_5.png" /></figure>

**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.

The particle trajectory follows a spiral, as shown in Fig. 3. The particle comes to rest after a typical time $\omega_B \tau_\infty \gg 1/b$.

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 $2\alpha\gamma_0^2\tau \approx 1$. 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.

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 $\Delta\tau$. Here the method is also second order in time, as expected.
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Comparison between the analytical solution, Eq. (31) (red solid line) and the numerical simulation (blue dots) for $\log \gamma_0 = 4$ and $\log b = -5$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Relative error of the position x and y as shown (see inset for legend). The error decreases with second order in $\Delta\tau$, as given by the violet line $\Delta\tau^2$ for $\log \gamma_0 = 4$ and $\log b = \{-5, -10\}$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**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.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Orbit in the electric drift frame $(x', y')$ for $\log \gamma_0 = 4$ and $\log b = -5$ for different proper time steps $\log(\omega_B \Delta\tau) = \{-1, -2\}$.</figcaption></figure>

### 3.3. Cross electric and magnetic field

The cross electric and magnetic field configuration is a stringent test for an ultra-relativistic particle pusher. If the electric field strength is less than the magnetic field strength $E < cB$, 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 \wedge B/B^2$. Therefore, the velocity is $\beta_E = v_E/c = E/cB$ and the corresponding final Lorentz factor $\gamma_\infty = (1 - \beta_E^2)^{-1/2}$.

We performed simulations with $E/cB = 0.999$ and initial Lorentz factors $\log \gamma_0 = \{0, 4, 8\}$. The final Lorentz factor is $\gamma_\infty \approx 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 ultra-relativistic motion with high precision.

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 \gamma_0 = 4$ and $\log b = -5$, using different proper time steps such as $\log(\omega_B \Delta\tau) = \{-1, -2\}$. The analytical solution found from the appropriate parameters in Eq. (31) is overlapped in red solid lines.

Figure 9 shows the relative error in the $x'$ and $y'$ position depending on the proper time step $\Delta\tau$. The scheme converges to second order in proper time step.
<figure><img src="image_10.png" /><figcaption>Fig. 9. Relative error of the particle position associated with the electric drift motion shown in Fig. 8.</figcaption></figure>

### 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 $I_2 \neq 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 four-velocity reads

$$ \frac{u^t}{c} = \gamma(\tau) = \frac{\gamma_0 c \cosh(\omega_E \tau) + u_z^0 \sinh(\omega_E \tau)}{\sqrt{\gamma_0^2 c^2 - u_\parallel^2 - u_\perp^2 e^{-2\alpha\tau}}}, \quad (32a) $$

$$ \frac{u^x}{c} = \frac{u_x^0 \cos(\omega_B \tau) + u_y^0 \sin(\omega_B \tau)}{\sqrt{(\gamma_0^2 c^2 - u_\parallel^2) e^{2\alpha\tau} - u_\perp^2}}, \quad (32b) $$

$$ \frac{u^y}{c} = \frac{u_y^0 \cos(\omega_B \tau) - u_x^0 \sin(\omega_B \tau)}{\sqrt{(\gamma_0^2 c^2 - u_\parallel^2) e^{2\alpha\tau} - u_\perp^2}}, \quad (32c) $$

$$ \frac{u^z}{c} = \frac{u_z^0 \cosh(\omega_E \tau) + \gamma_0 c \sinh(\omega_E \tau)}{\sqrt{\gamma_0^2 c^2 - u_\parallel^2 - u_\perp^2 e^{-2\alpha\tau}}}. \quad (32d) $$

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 $\alpha = \tau_m (\lambda_E^2 + \lambda_B^2)$, including non-vanishing electric and magnetic contributions.

After a transition time, the gyro-motion 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 $\lambda_E$ and $\lambda_B$.

Figure 10 shows the analytic evolution of the Lorentz factor for a particle starting with only a perpendicular velocity component such that $\log \gamma_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

<figure><img src="image_11.png" /><figcaption>Fig. 10. Analytic evolution of the Lorentz factor with initial condition $\log \gamma_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.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 11. Evolution of the Lorentz factor for a parallel electromagnetic field configuration with initial Lorentz factors $\log \gamma_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.</figcaption></figure>

field $E_0 \ll 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 ultra-relativistic speeds on a timescale $1/\omega_E \gg 1/\omega_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 \gg 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 force-free neutron star magnetospheres, and initial Lorentz factors $\log \gamma_0 = \{0, 4, 8\}$ and $\log b = \{0, -5, -10, -15\}$.

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 four-velocity component becomes comparable to the perpendicular four-velocity component for an initial velocity perpendicular to the magnetic
<figure><img src="image_13.png" /><figcaption>**Fig. 12.** Alignment of the velocity vector with the radiation reaction limit direction given by $\theta = 0^\circ$ for different strengths of the parallel electric field component $E_\parallel$ compared to the perpendicular component $E_\perp$ for $E_\perp/cB = 0.1$ (dashed lines) and $E_\perp/cB = 0.999$ (solid lines). The damping parameter is $\log b = -5$ and the initial Lorentz factor $\log \gamma_0 = 4$.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. 13.** Decrease in the Lorentz factor due to radiation reaction when drifting in a dipole magnetic field with initial Lorentz factor $\log \gamma_0 = 4$ and different damping constants.</figcaption></figure>

field line. This represents the worst case, useful for comparison to the radiation reaction limit regime.

### 3.5. Almost cross field

In a plasma-filled 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_\parallel \ll E_\perp$. As an application towards this configuration, we computed the motion of a particle in an almost cross electric and magnetic field with $\log(E_\parallel \ll E_\perp) = \{-1, -2, -3, -4\}$ and different ratios $E_\perp/cB = \{0.1, 0.999\}$. The particle starts with an initial velocity in a plane perpendicular to **B** and Lorentz factor $\log \gamma_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 $\theta$ 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_\perp/cB$, but depends strongly on the ratio $E_\parallel/E_\perp$. 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.

<figure><img src="image_15.png" /><figcaption>**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.</figcaption></figure>

### 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 convergence 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 $\omega_B$ and the particle initial Lorentz factor is $\gamma_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.

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.
<figure><img src="image_16.png" /><figcaption>Fig. 15. Decrease in the Lorentz factor of particles subject to the mirror effect.</figcaption></figure>

### 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(\tau_m \omega_B) \lesssim -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(\tau_m \omega_B) = -5$, the cyclotron motion is rapidly damped and the particle falls onto the star (green trajectory). For $log(\tau_m \omega_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).

## 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 $\omega = \Omega_*$ 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

<figure><img src="image_17.png" /><figcaption>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.</figcaption></figure>

the stellar rotation, which is many orders of magnitude smaller than the cyclotron frequency.

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 $\gamma_0 = 10^3$ or starting at rest. The neutron star obliquity is set to $\chi = \{0^\circ, 30^\circ, 60^\circ, 90^\circ, 120^\circ, 150^\circ, 180^\circ\}$. 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 half-way between the surface and the light-cylinder (a geometric average); and at the light cylinder, thus $r_0 = \{R_*, \sqrt{R_*} r_L, r_L\}$. 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
<table><caption>**Table 1.** Typical period and surface magnetic field strength of millisecond pulsars, young pulsars, and magnetars.</caption><thead><tr><th>Neutron star</th><th>P<sub>*</sub> (s)</th><th>log B<sub>*</sub> (T)</th><th>log a<sub>B</sub></th><th>log a<sub>E</sub></th><th>-log b</th></tr></thead><tbody><tr><td>Millisecond</td><td>0.005</td><td>5</td><td>13.1 / 9.9 / 9.5</td><td>11.8 / 8.6 / 8.2</td><td>20.1 / 23.4 / 22.5</td></tr><tr><td>Young</td><td>1</td><td>8</td><td>18.4 / 15.2 / 14.8</td><td>14.8 / 11.6 / 11.2</td><td>22.4 / 25.7 / 24.8</td></tr><tr><td>Magnetar</td><td>10</td><td>10</td><td>21.4 / 18.2 / 17.8</td><td>16.8 / 13.6 / 13.2</td><td>23.4 / 26.7 / 25.8</td></tr></tbody></table>

**Notes.** The relevant parameters without dimension are given by the strength parameters for the magnetic field $a_B$ and for the electric field $a_E$ and the damping parameter $b$ for electrons / protons / iron nuclei.

<table><caption>**Table 2.** Maximum Lorentz factor orders of magnitude from conservative arguments about neutron star magnetospheres.</caption><thead><tr><th rowspan="2">Neutron star</th><th colspan="3">log γ<sup>full</sup><sub>max</sub> / log γ<sup>pc</sup><sub>max</sub></th></tr><tr><th>Electron</th><th>Proton</th><th>Iron</th></tr></thead><tbody><tr><td>Millisecond</td><td>10.5 / 9.3</td><td>7.3 / 6.0</td><td>7.0 / 5.7</td></tr><tr><td>Young</td><td>11.2 / 7.7</td><td>8.0 / 4.4</td><td>7.7 / 4.1</td></tr><tr><td>Magnetar</td><td>12.2 / 7.7</td><td>9.0 / 4.4</td><td>8.7 / 4.1</td></tr></tbody></table>

**Notes.** Values for full potential drops are given to the left of the slash (/) and for polar cap potential drops to the right in logarithmic scale.

accelerating field, and thus

$$ \gamma_{\max}^{\text{full}} \approx \frac{q \Omega_* B_* R_*^2}{m c^2} = \frac{R_*}{r_L} \frac{R_*}{r_B}, \qquad (33) $$

where $r_B = c/\omega_B$ is the non-relativistic Larmor radius. If the accelerating potential is only available across the polar caps as expected from nearly force-free magnetosphere models, the maximum energy corresponds to

$$ \gamma_{\max}^{\text{pc}} \approx \frac{q \Omega_*^2 B_* R_*^3}{m c^3} = \left(\frac{R_*}{r_L}\right)^2 \frac{R_*}{r_B} \approx \frac{R_*}{r_L} \gamma_{\max}^{\text{full}}, \qquad (34) $$

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 $\gamma_{\max}^{\text{full}}$ 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.

### 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 $\chi$. 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 $\gamma_f \sim 10^9$ in the LLR approximation for millisecond pulsars. This is, however, two orders of magnitude less than without the radiation reaction where $\gamma_f \sim 10^{11}$. Second, as expected, protons and iron nuclei acquire much less energy, only about $\gamma_f \sim 10^6$ for millisecond pulsars, whether LLR is used or not. For young pulsars, electrons also reach $\gamma_f \sim 10^9$ in the LLR regime instead of $\gamma_f \sim 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 $\gamma_f \sim 10^4-10^{4.5}$. For magnetars, the radiation reaction remains negligible irrespective of the nature of each species. Electrons reach energies up to $\gamma_f \sim 10^{7.5}$, whereas protons and iron nuclei $\gamma_f \sim 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.

### 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 $\gamma_f \sim 10^{11}$ for magnetars, but only $\gamma_f \sim 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 $\gamma_f \sim 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 $\gamma_f \sim 10^7$, whereas young pulsars and magnetars one decade more up to $\gamma_f \sim 10^8$, whether the radiation reaction is included or not. Electrons are trapped with similar Lorentz factors, although slightly less for millisecond pulsars.
<figure><img src="image_18.png" /><figcaption>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.</figcaption></figure>

## 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

$$ \gamma_{\max} \approx 2 (a_{B_L}/\pi)^{2/3}, \qquad (35) $$

where $a_{B_L}$ 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 $\gamma_f \approx 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
<figure><img src="image_19.png" /><figcaption>Fig. 18. Same as Fig. 17, but for crashed particles.</figcaption></figure>

with Lorentz factor higher than $10^{9.1}$ within the magnetosphere. Equation (35) is satisfied for non-null 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

$$ \gamma_{\max} \approx a_{BL}. \qquad (36) $$

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
<figure><img src="image_20.png" /><figcaption>**Fig. 19.** Same as Fig. 17, but for trapped particles.</figcaption></figure>

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 $\gamma_f \lesssim 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 $\gamma_f \lesssim 10^9$. The flaw in these studies is that particles evolve in a prescribed external field without
<table><caption>**Table 3.** Maximum Lorentz factors $\gamma_{\text{max}}$ for the three kinds of particles: electrons / protons / iron nuclei.</caption><thead><tr><th rowspan="2">Neutron star</th><th rowspan="2">log a<sub>B<sub>L</sub></sub></th><th colspan="3">log γ<sub>max</sub></th></tr><tr><th>Crashed</th><th>Trapped</th><th>Escaped</th></tr></thead><tbody><tr><td>Millisecond</td><td>9.3 / 6.0 / 5.7</td><td>7.1 / 7.1 / 6.8</td><td>6.5 / 7.1 / 6.8</td><td>9.1 / 6.0 / 5.6</td></tr><tr><td>Young</td><td>7.6 / 4.4 / 4.1</td><td>6.0 / 7.8 / 7.5</td><td>7.8 / 7.8 / 7.5</td><td>8.0 / 4.1 / 3.7</td></tr><tr><td>Magnetar</td><td>7.6 / 4.4 / 4.1</td><td>6.4 / 8.1 / 7.7</td><td>8.2 / 8.1 / 7.7</td><td>7.1 / 3.2 / 2.9</td></tr></tbody></table>

Notes. The value of the strength parameter at the light cylinder is also given.

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 $\gamma_f \le 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 that were one order of magnitude lower with $\gamma_f \le 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 $\gamma_f \le 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 non-linear in a significant radiation reaction regime. General relativity does not significantly change these conclusions, as shown by Philippov & Spitkovsky (2018), who included frame-dragging effects and found $\gamma_f \le 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 \times 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 understand the origin of ultra-high-energy 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 $\gamma_f \le 700$ and for protons of about $\gamma_f \le 40$. These state-of-the-art results emphasise the difficulty in moving towards a realistic and self-consistent description of neutron star magnetospheres. The main bottleneck is the particle pusher, which requires us to temporally resolve the gyro-motion. 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 gyro-frequency 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 gyro-frequency, 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_\pm = c$, in the LLR approach particles are kicked with high initial Lorentz factors $\gamma_0 = 10^3$ and a velocity parallel to $v_\pm$ 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 $\sqrt{R_*} r_L$ (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.

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
<figure><img src="image_21.png" /></figure>

**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 $\sqrt{R_*} r_L$ in the *middle row*, and at a distance $r = r_L$ in the *bottom row*.

**Fig. 21.** Same as Fig. 20, but for protons.

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 time-consuming full integration of the equation of motion in LLR and an artificial and unrealistic down-scaling 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 $\Delta t^2$ behaviour is shown as blue filled circles. In the previous simulations we fixed the time step to $\Delta t \approx 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

<figure><img src="image_22.png" /></figure>

**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 $\Delta t^2$ decrease is shown as blue filled circles.

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 three-velocity vector. However, the 3D version of the LLR equation also
<figure><img src="image_23.png" /><figcaption>Fig. 23. Comparison between the radiation reaction limit (dotted symbols), the LLR (thick solid lines), and the approximated LLR motion (thin dashed lines).</figcaption></figure>

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 ultra-relativistic speeds). The expressions for the four-velocity then become

$$u_E = \lambda_B \frac{u_E^0 \cosh(\lambda_E \tau) + \tilde{F} u_E^0 \sinh(\lambda_E \tau) / \lambda_E}{\sqrt{(\lambda_E^2 + \lambda_B^2) |u_E^0|^2 + (\lambda_B^2 - (\lambda_E^2 + \lambda_B^2) |u_E^0|^2) e^{-2 \lambda_B^2 \tau_0 \tau}}}, \quad (37a)$$

$$u_B = \lambda_E \frac{u_B^0 \cos(\lambda_B \tau) + \tilde{F} u_B^0 \sin(\lambda_B \tau) / \lambda_B}{\sqrt{(\lambda_E^2 + \lambda_B^2) |u_B^0|^2 + (\lambda_E^2 - (\lambda_E^2 + \lambda_B^2) |u_B^0|^2) e^{2 \lambda_B^2 \tau_0 \tau}}}, \quad (37b)$$

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.

## 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 gyro-motion, 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.

*Acknowledgements.* I am grateful to the referee for helpful comments and suggestions that helped to improve this work. This work has been supported by the CEFIPRA grant IFC/F5904-B/2018 and ANR-20-CE31-0010. 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 Alsac calcul/Big Data.

## References

* Abraham, M. 1902, Ann. Phys., 315, 105

* Abraham, M. 1904, Ann. Phys., 319, 236

* Boghosian, B. M. 1987, PhD Thesis, University of California

* Boris, J. 1970, Proceeding of Fourth Conference on Numerical Simulations of Plasmas

* Cerutti, B., Philippov, A., Parfrey, K., & Spitkovsky, A. 2015, MNRAS, 448, 606

* Cerutti, B., Philippov, A. A., & Spitkovsky, A. 2016, MNRAS, 457, 2401

* Cerutti, B., Philippov, A. A., & Dubus, G. 2020, A&A, 642, A204

* Deutsch, A. J. 1955, Ann. d'Astrophys., 18, 1

* Dirac, P. A. M. 1938, Proc. R. Soc. Lond. Ser. A, 167, 148

* Elkina, N. V., Fedotov, A. M., Herzing, C., & Ruhl, H. 2014, Phys. Rev. E, 89, 053315

* Ferrari, A., & Trussoni, E. 1974, A&A, 36, 267

* Finkbeiner, B., Herold, H., Ertl, T., & Ruder, H. 1989, A&A, 225, 479

* Finkbeiner, B., Herold, H., & Ruder, H. 1990, A&A, 238, 462

* Fulton, T., & Rohrlich, F. 1960, Ann. Phys., 9, 499

* Gordon, D. F., & Hafizi, B. 2021, Comput. Phys. Commun., 258, 107628

* Gordon, D. F., Hafizi, B., & Palastro, J. 2017a, AIP Conf. Proc., 1812, 050002

* Gordon, D. F., Palastro, J. P., & Hafizi, B. 2017b, Phys. Rev. A, 95, 033403

* Guépin, C., Cerutti, B., & Kotera, K. 2020, A&A, 635, A138
Hadad, Y., Labun, L., Rafelski, J., et al. 2010, Phys. Rev. D, 82, 096012

Heintzmann, H., & Schrüfer, E. 1973, Phys. Lett. A, 43, 287

Herold, H., Ertl, T., & Ruder, H. 1985, Mitteilungen der Astronomischen Gesellschaft Hamburg, 63, 174

Kalapotharakos, C., Brambilla, G., Timokhin, A., Harding, A. K., & Kazanas, D. 2018, ApJ, 857, 44

Kelner, S. R., Prosekin, A. Y., & Aharonian, F. A. 2015, AJ, 149, 33

Kulsrud, R. M. 1972, ApJ, 174, L25

Landau, L., & Lifschitz, E. 1989, Physique théorique: Tome 2, Théorie des champs (Moscou: Mir)

Laue, H., & Thielheim, K. O. 1986, ApJS, 61, 465

Li, F., Decyk, V. K., Miller, K. G., et al. 2021, J. Comput. Phys., 438, 110367

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)

Mestel, L., Robertson, J. A., Wang, Y. M., & Westfold, K. C. 1985, MNRAS, 217, 443

Michel, F., & Li, H. 1999, Phys. Rep., 318, 227

Olver, F. W. J. 2010, NIST Handbook of Mathematical Functions (Cambridge: Cambridge University Press)

Philippov, A. A., & Spitkovsky, A. 2018, ApJ, 855, 94

Philippov, A. A., Spitkovsky, A., & Cerutti, B. 2015, ApJ, 801, L19

Piazza, A. D. 2008, Lett. Math. Phys., 83, 305

Pétri, J. 2020, J. Plasma Phys., 86, 825860402

Pétri, J. 2021, MNRAS, 503, 2123

Rohrlich, F. 2007, Classical Charged Particles, 3rd edn. (Singapore: World Scientific Pub Co Inc)

Tomczak, I., & Pétri, J. 2020, J. Plasma Phys., 86, 825860401

Vranic, M., Martins, J. L., Fonseca, R. A., & Silva, L. O. 2016, Comput. Phys. Commun., 204, 141