© J. Pétri 2022

1. Introduction

Neutron stars are known to harbour ultra-strong magnetic fields close to or even above the quantum critical limit of B_c ≈ 4.4 × 10⁹ 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⁹. 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

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{u}^i}{\mathrm{d}\tau }& =\frac{q}{m}{F}^{\mathit{ik}}{u}_k+\frac{q{\tau }_{\mathrm{m}}}{m}{g}^i,\hfill \end{array}$$(1a)

$$\begin{array}{cc}\hfill g^i& ={\partial }_{\ell }{F}^{\mathit{ik}}{u}_k{u}^{\ell }+\frac{q}{m}(F^{\mathit{ik}}{F}_{k\ell }{u}^{\ell }+(F^{\ell m}{u}_m)(F_{\ell k}{u}^k)\frac{u^i}{c^2}),\hfill \end{array}$$(1b)

where q and m are the particle charge and rest mass, uⁱ its four-velocity, τ 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)

$$\begin{array}{c}\hfill {\tau }_{\mathrm{m}}=\frac{q^2}{6\pi {\epsilon }_{0}m{c}^3}.\end{array}$$(2)

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

$$\begin{array}{c}\hfill {\tau }_{\mathrm{e}}=\frac{2}{3}\frac{r_{\mathrm{e}}}{c}=6.26\times {10}^{-24}\mathrm{s}.\end{array}$$(3)

The typical timescale for the radiation reaction is therefore

$$\begin{array}{c}\hfill {\tau }_{\mathrm{m}}=\frac{2}{3}\frac{r_{\mathrm{m}}}{c}=\left(\frac{q^2/e^2}{m/m_{\mathrm{e}}}\right){\tau }_{\mathrm{e}}.\end{array}$$(4)

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

$$\begin{array}{c}\hfill {\tau }_{\mathrm{p}}=\frac{m_{\mathrm{e}}}{m_{\mathrm{p}}}{\tau }_{\mathrm{e}}=3.41\times {10}^{-27}\mathrm{s}.\end{array}$$(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_k^i$. 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

$$\begin{array}{c}\hfill \frac{\mathrm{d}u}{\mathrm{d}\tau }=Gu,\end{array}$$(6)

where the electromagnetic tensor F has been replaced by G = q F/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_k^i$ 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

$$\begin{array}{cc}& G{u}_E=\pm {\lambda }_E{u}_E,\hfill \end{array}$$(7a)

$$\begin{array}{cc}& G{u}_B=\pm i{\lambda }_B{u}_B.\hfill \end{array}$$(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)

$$\begin{array}{cc}& P=\frac{{\lambda }_B^{2}I+G^2}{{\lambda }_E^2+{\lambda }_B^2},\hfill \end{array}$$(8a)

$$\begin{array}{cc}& Q=\frac{{\lambda }_E^{2}I-G^2}{{\lambda }_E^2+{\lambda }_B^2},\hfill \end{array}$$(8b)

where I is the identity matrix. These operators are well defined only if ${\lambda }_E^2+{\lambda }_B^2\ne 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

$$\begin{array}{cc}& u_E=Pu,\hfill \end{array}$$(9a)

$$\begin{array}{cc}& u_B=Qu.\hfill \end{array}$$(9b)

The equation of motion decouples into two parts given by

$$\begin{array}{c}\hfill \frac{{\mathrm{d}}^2{u}_E}{\mathrm{d}{\tau }^2}=+{\lambda }_E^2{u}_E,\end{array}$$(10a)

$$\begin{array}{c}\hfill \frac{{\mathrm{d}}^2{u}_B}{\mathrm{d}{\tau }^2}=-{\lambda }_B^2{u}_B.\end{array}$$(10b)

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

$$\begin{array}{cc}& u_E(\tau )=u_E^{0}cosh({\lambda }_E\tau )+G{u}_E^0\frac{sinh({\lambda }_E\tau )}{{\lambda }_E},\hfill \end{array}$$(11a)

$$\begin{array}{cc}& u_B(\tau )=u_B^{0}cos({\lambda }_B\tau )+G{u}_B^0\frac{sin({\lambda }_B\tau )}{{\lambda }_B}.\hfill \end{array}$$(11b)

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

$$\begin{array}{cc}& \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 }}},\hfill \end{array}$$(12a)

$$\begin{array}{cc}& \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 }}},\hfill \end{array}$$(12b)

with $\alpha ={\tau }_{\text{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 τ ≈ 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 gyro-motion 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 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 ω 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₀ = c/ω, the timescale T₀ = 1/ω, the magnetic field strength B₀ = m ω/q, and the electric field strength E₀ = c B₀. 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:

$$\begin{array}{cc}\hfill a_B& =\frac{B}{B_0}=\frac{{\omega }_B}{\omega },\hfill \end{array}$$(13a)

$$\begin{array}{cc}\hfill a_E& =\frac{E}{E_0}=\frac{{\omega }_E}{\omega },\hfill \end{array}$$(13b)

$$\begin{array}{cc}\hfill b& =\omega {\tau }_{\mathrm{m}}.\hfill \end{array}$$(13c)

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

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{\stackrel{\sim }{u}}^i}{\mathrm{d}\stackrel{\sim }{\tau }}& ={\stackrel{\sim }{F}}^{\mathit{ik}}{\stackrel{\sim }{u}}_k+b{\stackrel{\sim }{g}}^i,\hfill \end{array}$$(14a)

$$\begin{array}{cc}\hfill {\stackrel{\sim }{g}}^i& ={\stackrel{\sim }{\partial }}_{\ell }{\stackrel{\sim }{F}}^{\mathit{ik}}{\stackrel{\sim }{u}}_k{\stackrel{\sim }{u}}^{\ell }+({\stackrel{\sim }{F}}^{\mathit{ik}}{\stackrel{\sim }{F}}_{k\ell }{\stackrel{\sim }{u}}^{\ell }+({\stackrel{\sim }{F}}^{\ell m}{\stackrel{\sim }{u}}_m)({\stackrel{\sim }{F}}_{\ell k}{\stackrel{\sim }{u}}^k){\stackrel{\sim }{u}}^i).\hfill \end{array}$$(14b)

The normalised and reduced Landau–Lifshitz equation reads

$$\begin{array}{c}\hfill \frac{\mathrm{d}{\stackrel{\sim }{u}}^i}{\mathrm{d}\stackrel{\sim }{\tau }}={\stackrel{\sim }{F}}^{\mathit{ik}}{\stackrel{\sim }{u}}_k+b({\stackrel{\sim }{F}}^{\mathit{ik}}{\stackrel{\sim }{F}}_{k\ell }{\stackrel{\sim }{u}}^{\ell }+({\stackrel{\sim }{F}}^{\ell m}{\stackrel{\sim }{u}}_m)({\stackrel{\sim }{F}}_{\ell k}{\stackrel{\sim }{u}}^k){\stackrel{\sim }{u}}^i).\end{array}$$(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/α, becomes larger than the characteristic timescale of our system, which is 1/ω. Expressed in quantities without dimension, we get ${\tau }_{\text{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 ₂F₁ 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:

$$\begin{array}{cc}& {\mathit{u}}^{n+1/2}=\mathit{L}(\mathrm{\Delta }\tau /2,{\mathit{u}}^n,\mathit{E}({\mathit{x}}^n),\mathit{B}({\mathit{x}}^n)),\hfill \end{array}$$(16a)

$$\begin{array}{cc}& {\mathit{x}}^{n+1}={\mathit{x}}^n+{\mathit{u}}^{n+1/2}\mathrm{\Delta }\tau ,\hfill \end{array}$$(16b)

$$\begin{array}{cc}& {\mathit{u}}^{n+1}=\mathit{L}(\mathrm{\Delta }\tau /2,{\mathit{u}}^{n+1/2},\mathit{E}({\mathit{x}}^{n+1}),\mathit{B}({\mathit{x}}^{n+1})).\hfill \end{array}$$(16c)

The superscript n refers to the proper time τⁿ = n Δτ, and the same for the half-integer 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 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

$$\begin{array}{c}\hfill {\mathit{v}}_{\pm }=\frac{\mathbf{E}\wedge \mathbf{B}\pm ({\mathbf{E}}_{\mathbf{0}}\mathbf{E}/\mathbf{c}+\mathbf{c}{\mathbf{B}}_{\mathbf{0}}\mathbf{B})}{E_0^2/c^2+B^2},\end{array}$$(17)

which corresponds to particles moving exactly at the speed of light. Here E₀ and B₀ are the strength of the electric and magnetic field in the frame where they are aligned. They are obtained from the electromagnetic invariants ${\mathcal{I}}_1={\mathbf{E}}^2-c^2{\mathbf{B}}^2=E_0^2-c^2{B}_0^2$ and ℐ₂ = c E ⋅ B = c E₀ B₀. Imposing E₀ ≥ 0 we find

$$\begin{array}{cc}& E_0^2=\frac{1}{2}({\mathcal{I}}_1+\sqrt{{\mathcal{I}}_1^2+4{\mathcal{I}}_2^2}),\hfill \end{array}$$(18a)

$$\begin{array}{cc}& c{B}_0=\text{sign}({\mathcal{I}}_2)\sqrt{E_0^2-{\mathcal{I}}_1}.\hfill \end{array}$$(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, |ℐ₂|≪|ℐ₁| and ℐ₁ < 0. The velocity then reduces to

$$\begin{array}{c}\hfill {\mathit{v}}_{\pm }\approx \frac{\mathbf{E}\wedge \mathbf{B}}{B^2}\pm \text{sign}(\mathbf{E}·\mathbf{B})\frac{\sqrt{c^2{B}^2-E^2}}{B^2}\mathbf{B}.\end{array}$$(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 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, 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

$$\begin{array}{c}\hfill {\mathit{v}}_{\pm }=\frac{\mathbf{E}\wedge \mathbf{B}}{B^2}\pm \frac{\sqrt{c^2{B}^2-E^2}}{B^2}\mathbf{B}.\end{array}$$(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 ≈ 10²⁰ and γ ≈ 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 E = E e_z, we obtain

$$\begin{array}{cc}& \frac{u^t}{c}=\gamma (\tau )=\frac{{\gamma }_{0}ccosh({\omega }_E\tau )+u_z^{0}sinh({\omega }_{\mathrm{E}}\tau )}{\sqrt{{\gamma }_0^2{c}^2-u_{\Vert }^2-u_{\perp }^2{e}^{-2\alpha \tau }}},\hfill \end{array}$$(21a)

$$\begin{array}{cc}& \frac{u^x}{c}=\frac{u_x^0}{\sqrt{({\gamma }_0^2{c}^2-u_{\Vert }^2)e^{2\alpha \tau }-u_{\perp }^2}},\hfill \end{array}$$(21b)

$$\begin{array}{cc}& \frac{u^y}{c}=\frac{u_y^0}{\sqrt{({\gamma }_0^2{c}^2-u_{\Vert }^2)e^{2\alpha \tau }-u_{\perp }^2}},\hfill \end{array}$$(21c)

$$\begin{array}{cc}& \frac{u^z}{c}=\frac{u_z^{0}cosh({\omega }_E\tau )+{\gamma }_{0}csinh({\omega }_E\tau )}{\sqrt{{\gamma }_0^2{c}^2-u_{\Vert }^2-u_{\perp }^2{e}^{-2\alpha \tau }}},\hfill \end{array}$$(21d)

with $u_{\parallel }=u_z^0$ the initial four-velocity component along E, u_⊥ the initial four-velocity component perpendicular to E, and $\alpha ={\tau }_{\text{m}}{\omega }_{\text{E}}^2$.

For a particle starting at rest, u_∥ = u_⊥ = 0 and γ₀ = 1, the four-velocity simplifies drastically to

$$\begin{array}{c}\hfill u^i=c(cosh({\omega }_E\tau ),0,0,sinh({\omega }_{\mathrm{E}}\tau )).\end{array}$$(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 ₂F₁ (Olver 2010) such that

$$\begin{array}{cc}& J_1(\tau )=\frac{e^{(\alpha +{\omega }_E)\tau }}{\alpha +{\omega }_E}{}_2{F}_1(\frac{1}{2},1+\frac{1}{2b};\frac{3}{2}+\frac{1}{2b};\frac{{\gamma }^2{e}^{2\alpha \tau }}{{\gamma }^2-1}),\hfill \end{array}$$(23a)

$$\begin{array}{cc}& J_2(\tau )=\frac{e^{(\alpha -{\omega }_E)\tau }}{\alpha -{\omega }_E}{}_2{F}_1(\frac{1}{2},1-\frac{1}{2b};\frac{3}{2}-\frac{1}{2b};\frac{{\gamma }^2{e}^{2\alpha \tau }}{{\gamma }^2-1}).\hfill \end{array}$$(23b)

Then the time and position are given by

$$\begin{array}{cc}\hfill t& =-\frac{\gamma \sqrt{{\gamma }^2(e^{2\alpha \tau }-1)+1}}{2({\gamma }^2-1)}(J_2(\tau )+J_1(\tau ))+C_0,\hfill \end{array}$$(24a)

$$\begin{array}{cc}\hfill x/c& =0,\hfill \end{array}$$(24b)

$$\begin{array}{cc}\hfill y/c& =\frac{1}{\alpha }arctan\left(\sqrt{\frac{{\gamma }^2(e^{2\alpha \tau }-1)+1}{{\gamma }^2-1}}\right)+C_2,\hfill \end{array}$$(24c)

$$\begin{array}{cc}\hfill 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,\hfill \end{array}$$(24d)

with C₀, C₂, C₃ 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 $\stackrel{\sim }{\tau }={\omega }_E\tau$. Therefore, the only relevant parameter apart from the initial conditions is b = τ_m ω_E and $\alpha \tau =b\stackrel{\sim }{\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_∥ = 0 and u_⊥ ≠ 0. We chose different initial Lorentz factors in the set log γ₀ = {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 γ₀ ≫ 1. The time derivative of the Lorentz factor is always negative for γ₀ > 1 because

$$\begin{array}{c}\hfill \frac{\mathrm{d}\gamma (\tau )}{\mathrm{d}\tau }=-\alpha {\gamma }_0({\gamma }_0^2-1)<0,\end{array}$$(25)

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 γ₀. This occurs because the perpendicular motion is strongly damped, $\underset{\tau \to +\infty }{lim}{u}_{\perp }\to 0$, and only the parallel velocity u_∥ survives at large times with $\underset{\tau \to +\infty }{lim}{u}_{\Vert }\to csinh({\omega }_E\tau )$. In between, the normalised time remains small and the Lorentz factor can be approximated by

$$\begin{array}{c}\hfill \gamma (\tau )\approx \frac{{\gamma }_{0}cosh({\omega }_E\tau )}{\sqrt{1+2{\gamma }_0^2\alpha \tau }}.\end{array}$$(26)

Therefore, before the acceleration phase starts, there is a deceleration step arising at time ${\omega }_{\text{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 γ₀ = 10⁴ and b = 10⁻¹⁰ (orange dots) or γ₀ = 10⁸ and b = 10⁻¹⁵, 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 γ₀ because at time ω_E τ ≈ 1 it braked to a momentum

$$\begin{array}{c}\hfill u_{\perp }(\tau )\approx \frac{u_{\perp }^0}{\sqrt{1+2{\gamma }_0^{2}b}}.\end{array}$$(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 Δτ for log γ₀ = 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.

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 four-velocity in a magnetic field directed along the z-axis with B = B e_z is given by

$$\begin{array}{cc}& \frac{u^t}{c}=\gamma (\tau )=\frac{{\gamma }_{0}c}{\sqrt{{\gamma }_0^2{c}^2-u_{\Vert }^2-u_{\perp }^2{e}^{-2\alpha \tau }}},\hfill \end{array}$$(28a)

$$\begin{array}{cc}& \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_{\Vert }^2)e^{2\alpha \tau }-u_{\perp }^2}},\hfill \end{array}$$(28b)

$$\begin{array}{cc}& \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_{\Vert }^2)e^{2\alpha \tau }-u_{\perp }^2}},\hfill \end{array}$$(28c)

$$\begin{array}{cc}& \frac{u^z}{c}=\frac{u_z^0}{\sqrt{{\gamma }_0^2{c}^2-u_{\Vert }^2-u_{\perp }^2{e}^{-2\alpha \tau }}},\hfill \end{array}$$(28d)

with $u_{\parallel }=u_z^0$ the initial four-velocity component along B, u_⊥ the initial four-velocity component perpendicular to B, and $\alpha ={\tau }_{\text{m}}{\omega }_{\text{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 α τ.

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

$$\begin{array}{c}\hfill r_B=\frac{c}{{\omega }_B}.\end{array}$$(29)

Integrating the four-velocity vector, an exact analytical expression for the four-position is computed with the help of the hypergeometric functions ₂F₁. Introducing the complex functions

$$\begin{array}{cc}& H_1(\tau )=e^{+i{\omega }_B\tau }{}_2{F}_1(\frac{1}{2},+\frac{i}{2b};1+\frac{i}{2b};\frac{{\gamma }^2{e}^{-2\alpha \tau }}{{\gamma }^2-1}),\hfill \end{array}$$(30a)

$$\begin{array}{cc}& H_2(\tau )=e^{-i{\omega }_B\tau }{}_2{F}_1(\frac{1}{2},-\frac{i}{2b};1-\frac{i}{2b};\frac{{\gamma }^2{e}^{-2\alpha \tau }}{{\gamma }^2-1}).\hfill \end{array}$$(30b)

the solution reads

$$\begin{array}{cc}\hfill t& =\frac{1}{\alpha }{tanh}^{-1}\left(\frac{\gamma e^{\alpha \tau }}{\sqrt{{\gamma }^2(e^{2\alpha \tau }-1)+1}}\right)+C_0,\hfill \end{array}$$(31a)

$$\begin{array}{cc}\hfill x/r_B& =\frac{H_1(\tau )+H_2(\tau )}{2i}+C_1,\hfill \end{array}$$(31b)

$$\begin{array}{cc}\hfill y/r_B& =\frac{H_1(\tau )-H_2(\tau )}{2}+C_2,\hfill \end{array}$$(31c)

$$\begin{array}{cc}\hfill z/r_B& =0,\hfill \end{array}$$(31d)

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

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 ultra-relativistic 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². Therefore, the velocity is β_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 γ₀ = {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 ultra-relativistic 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 γ₀ = 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 ℐ₂ ≠ 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

$$\begin{array}{cc}& \frac{u^t}{c}=\gamma (\tau )=\frac{{\gamma }_{0}ccosh({\omega }_{\mathrm{E}}\tau )+u_z^{0}sinh({\omega }_{\mathrm{E}}\tau )}{\sqrt{{\gamma }_0^2{c}^2-u_{\Vert }^2-u_{\perp }^2{e}^{-2\alpha \tau }}},\hfill \end{array}$$(32a)

$$\begin{array}{cc}& \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_{\Vert }^2)e^{2\alpha \tau }-u_{\perp }^2}},\hfill \end{array}$$(32b)

$$\begin{array}{cc}& \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_{\Vert }^2)e^{2\alpha \tau }-u_{\perp }^2}},\hfill \end{array}$$(32c)

$$\begin{array}{cc}& \frac{u^z}{c}=\frac{u_z^{0}cosh({\omega }_E\tau )+{\gamma }_{0}csinh({\omega }_{\mathrm{E}}\tau )}{\sqrt{{\gamma }_0^2{c}^2-u_{\Vert }^2-u_{\perp }^2{e}^{-2\alpha \tau }}}.\hfill \end{array}$$(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 }_{\text{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 λ_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 γ₀ = 8. The damping factor is set to log b = {0, −5, −10, −15} and the electric field strength E₀ is varied relative to B₀ such that log(E₀/c B₀) = {−2, 0, 2} (solid lines, dashed lines, and dotted lines, respectively). For a weak electric field E₀ ≪ c B₀ 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/ω_E ≫ 1/ω_B. The particle performs many orbits before being deflected along the parallel direction (e_z). Increasing E₀ 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₀ ≫ c B₀ electric acceleration quickly sets in. Figure 11 shows some results of numerical simulations with a weak electric field log(E₀/c B₀) = − 2, pertinent for almost force-free neutron star magnetospheres, and initial Lorentz factors log γ₀ = {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 E0 relative to B0 such that log(E0/c B0) = {−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(E0/c B0) = − 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 four-velocity component becomes comparable to the perpendicular four-velocity 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 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_∥ ≪ 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 γ₀ = 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 γ₀ = 10⁴. 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⁸ T; millisecond pulsars with period P_* = 5 ms and B_* = 10⁵ T; and magnetars with period P_* = 10 s and B_* = 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.

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 γ₀ = 10³ 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₀: 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_{\ast },\sqrt{R_{\ast }{r}_{\mathrm{L}}},r_{\mathrm{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 accelerating field, and thus

$$\begin{array}{c}\hfill {\gamma }_{\max }^{\mathrm{full}}\approx \frac{q{\mathrm{\Omega }}_{\ast }{B}_{\ast }{R}_{\ast }^2}{m{c}^2}=\frac{R_{\ast }}{r_{\mathrm{L}}}\frac{R_{\ast }}{r_B},\end{array}$$(33)

where r_B = c/ω_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

$$\begin{array}{c}\hfill {\gamma }_{\max }^{\mathrm{pc}}\approx \frac{q{\mathrm{\Omega }}_{\ast }^2{B}_{\ast }{R}_{\ast }^3}{m{c}^3}={\left(\frac{R_{\ast }}{r_{\mathrm{L}}}\right)}^2\frac{R_{\ast }}{r_{\mathrm{B}}}\approx \frac{R_{\ast }}{r_{\mathrm{L}}}{\gamma }_{\max }^{\mathrm{full}},\end{array}$$(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 }^{\mathrm{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 χ. 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⁹ in the LLR approximation for millisecond pulsars. This is, however, two orders of magnitude less than without the radiation reaction where γ_f ∼ 10¹¹. Second, as expected, protons and iron nuclei acquire much less energy, only about γ_f ∼ 10⁶ for millisecond pulsars, whether LLR is used or not. For young pulsars, electrons also reach γ_f ∼ 10⁹ in the LLR regime instead of γ_f ∼ 10¹¹ 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⁴ − 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³ − 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¹¹ 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.

Fig. 18. Same as Fig. 17, but for crashed particles.

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⁷, whereas young pulsars and magnetars one decade more up to γ_f ∼ 10⁸, whether the radiation reaction is included or not. Electrons are trapped with similar Lorentz factors, although slightly less for millisecond pulsars.

Fig. 19. Same as Fig. 17, but for trapped particles.

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

$$\begin{array}{c}\hfill {\gamma }_{\max }\approx 2{(a_{B_{\mathrm{L}}}/\pi )}^{2/3},\end{array}$$(35)

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

$$\begin{array}{c}\hfill {\gamma }_{\max }\approx a_{B_{\mathrm{L}}}.\end{array}$$(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 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⁹ 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⁹. 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³ 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². 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³ 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 γ_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⁸ and slightly larger for protons (almost 10⁹), but for faster rotations in a field of an oblique rotating dipole with strengths of 5 × 10⁶ 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⁹ and for protons of about 10⁶.

Hadron acceleration has been much less discussed in this context, but it is equally important to understanding 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 γ_f ≲ 700 and for protons of about γ_f ≲ 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_± = c, in the LLR approach particles are kicked with high initial Lorentz factors γ₀ = 10³ 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 $\sqrt{R_{\ast }{r}_{\mathrm{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.

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_{\ast }{r}_{\mathrm{L}}}$ in the middle row, and at a distance r = rL 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.

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

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 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₀/r_L in the legend. For reference the Δt² behaviour is shown as blue filled circles. In the previous simulations we fixed the time step to Δt ≈ 10⁻⁴, 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 three-velocity 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 ultra-relativistic speeds). The expressions for the four-velocity then become

$$\begin{array}{cc}& u_E={\lambda }_B\frac{u_E^{0}cosh({\lambda }_E\tau )+\stackrel{\sim }{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 }}},\hfill \end{array}$$(37a)

$$\begin{array}{cc}& u_B={\lambda }_E\frac{u_B^{0}cos({\lambda }_B\tau )+\stackrel{\sim }{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 }_E^2{\tau }_0\tau }}},\hfill \end{array}$$(37b)

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 r0/rL in the legend. The Δt2 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⁹ 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.

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/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 Alsacalcul/Big Data.

References

All Tables

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 E0 relative to B0 such that log(E0/c B0) = {−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(E0/c B0) = − 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 $\sqrt{R_{\ast }{r}_{\mathrm{L}}}$ in the middle row, and at a distance r = rL 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 r0/rL in the legend. The Δt2 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