Open Access
Issue
A&A
Volume 666, October 2022
Article Number A5
Number of page(s) 18
Section Astrophysical processes
DOI https://doi.org/10.1051/0004-6361/202243634
Published online 27 September 2022

© J. Pétri 2022

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe-to-Open model. Subscribe to A&A to support open access publication.

1. Introduction

Neutron stars are known to harbour ultra-strong magnetic fields close to or even above the quantum critical limit of Bc ≈ 4.4 × 109 T. The subclass of magnetars usually sustains field strengths well above this value of Bc. These stars are therefore able to accelerate leptons and hadrons to extremely relativistic regimes of very high Lorentz factors γ ≈ 109. 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

(1a)

(1b)

where q and m are the particle charge and rest mass, ui its four-velocity, τ its proper time, Fik the electromagnetic or Faraday tensor, c the speed of light, and τm the light crossing time across the particle classical radius rm (within a factor unity)

(2)

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

(3)

The typical timescale for the radiation reaction is therefore

(4)

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

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

(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 uB and uE, such that u = uE + uB. The real eigenvalues of are ±λE, whereas the imaginary eigenvalues are ±iλB, because λE and λB are real and positive numbers, with dimensions similar to pulsation thus in 1/s. Then each vector uE and uB remains in an eigensubspace satisfying

(7a)

(7b)

The vector components uE and uB are obtained by defining the projection operators onto the subspaces E and B by (Boghosian 1987)

(8a)

(8b)

where I is the identity matrix. These operators are well defined only if . 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

(9a)

(9b)

The equation of motion decouples into two parts given by

(10a)

(10b)

The exact analytical solutions with initial conditions and are

(11a)

(11b)

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

(12a)

(12b)

with . These expressions are similar to the original formulas found by Heintzmann & Schrüfer (1973). The radiation reaction effect becomes perceptible after a time τ ≈ 1/α. The component uE is associated with the accelerating motion induced by the electric field, whereas the uB component is related to the gyro-motion in the magnetic field. When α vanishes, the radiation reaction effect disappears. The denominators in uE and uB 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 L0 = c/ω, the timescale T0 = 1/ω, the magnetic field strength B0 = mω/q, and the electric field strength E0 = cB0. Normalised quantities will be overlaid with a tilde symbol.

The two important parameters defining the family of solutions are the field strength parameters aB and aE and the radiation reaction efficiency ωτm according to the following definitions:

(13a)

(13b)

(13c)

Introducing the weighted and normalised electromagnetic field tensor by and a normalised time , the Landau-Lifshitz Eq. (1) is rewritten without dimensions as

(14a)

(14b)

The normalised and reduced Landau–Lifshitz equation reads

(15)

The particle four-velocity depends only on the strength parameters aB and aE and on the radiation reaction strength parameter b. Therefore, it is unnecessary to compute trajectories for different particles possessing the same numbers aB, aE, b. The only differences are reflected in the physical timescale and space scales involved.

Generally speaking, we admit that the radiation reaction is negligible whenever the timescale of damping, given by 1/α, becomes larger than the characteristic timescale of our system, which is 1/ω. Expressed in quantities without dimension, we get . Therefore, the relevant parameter to quantify radiation reaction is not b, but the combination of b and the strength parameters aB and aE. 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 (ex, ey, ez).

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 2F1 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:

(16a)

(16b)

(16c)

The superscript n refers to the proper time τn = 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

(17)

which corresponds to particles moving exactly at the speed of light. Here E0 and B0 are the strength of the electric and magnetic field in the frame where they are aligned. They are obtained from the electromagnetic invariants and ℐ2 = cE ⋅ B = cE0B0. Imposing E0 ≥ 0 we find

(18a)

(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, |ℐ2|≪|ℐ1| and ℐ1 < 0. The velocity then reduces to

(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, E0 = 0, the radiated power also vanishes, and the particle moves along the field lines with v± = ±cB/B. Second, if the electric field is orthogonal to the magnetic field, E ⋅ B = 0 and E < cB, the particle motion is decomposed into an electric drift and a motion along B such that

(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 aB ≈ 1020 and γ ≈ 1010. 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 = Eez, we obtain

(21a)

(21b)

(21c)

(21d)

with the initial four-velocity component along E, u the initial four-velocity component perpendicular to E, and .

For a particle starting at rest, u = u = 0 and γ0 = 1, the four-velocity simplifies drastically to

(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 2F1 (Olver 2010) such that

(23a)

(23b)

Then the time and position are given by

(24a)

(24b)

(24c)

(24d)

with C0, C2, C3 complex constants of integration to satisfy the initial conditions.

Returning to Eq. (21), the typical electric acceleration timescale is τacc ∼ 1/ωE. On the other hand, the radiation damping timescale is τrad ∼ 1/α. The ratio of the two timescales is therefore τacc/τrad ∼ τmωE = b. As expected, for small damping parameters b ≪ 1, the acceleration time is much shorter than the radiative damping and the particle is accelerated as if it would not radiate, until the time τrad ∼ τacc/b ≫ τacc. We note that this rough estimate needs to be corrected by taking into account the initial Lorentz factor, as discussed below.

To perform the simulations we use the characteristic frequency ωE as normalisation, which leads to a normalised proper time . Therefore, the only relevant parameter apart from the initial conditions is b = τmωE and . Particles starting at rest or possessing an initial velocity directed along the electric field do not suffer from the radiative force. Consequently, as a typical example, particles start with an initial velocity perpendicular to the electric field, meaning u = 0 and u ≠ 0. We chose different initial Lorentz factors in the set log γ0 = {0, 4, 8}. The damping factor is given by log b = {0, −5, −10, −15}. As output of the simulations, we plot the Lorentz factor increasing according to Eq. (21a) and shown in Fig. 1. For a particle starting at rest, whatever the damping parameter b, the solution is always equal to Eq. (22) with an acceleration arising around the time ωEτ ∼ 1. This configuration is particular and is not impacted by the radiation reaction. More interesting are the particles starting with a substantial kick velocity and high Lorentz factors γ0 ≫ 1. The time derivative of the Lorentz factor is always negative for γ0 > 1 because

(25)

thumbnail Fig. 1.

Increase in the Lorentz factor due to radiation reaction for different initial Lorentz factors log γ0 = {0, 4, 8} and different damping factors log b = {0, −5, −10, −15}. The vertical lines show the time when the damping sets in before the electric acceleration phase starts. The coloured points show the simulation results and the black solid lines correspond to the exact analytical solutions.

meaning that the particle first decelerates due to the radiative friction. At large times, when ατ ≫ 1 and ωEτ ≫ 1, the Lorentz factor behaves as γ(τ)≈cosh(ωEτ), losing its information about the initial state. It resembles the motion of a particle starting at rest, independently of γ0. This occurs because the perpendicular motion is strongly damped, , and only the parallel velocity u survives at large times with . In between, the normalised time remains small and the Lorentz factor can be approximated by

(26)

Therefore, before the acceleration phase starts, there is a deceleration step arising at time . These values agree with the curves in Fig. 1. If , the radiation reaction force has no time to set in and the motion tends to a purely accelerated regime given by Eq. (22). This is the case with γ0 = 104 and b = 10−10 (orange dots) or γ0 = 108 and b = 10−15, which is just on the edge of this condition, showing a weak deceleration right before the electric boost (blue dots). The perpendicular momentum decrease is not necessarily significant before the acceleration; it is controlled by b and γ0 because at time ωEτ ≈ 1 it braked to a momentum

(27)

Thus, the radiation reaction again impacts the motion if . Consequently, it is the combination that controls the damping efficiency, not b alone found from the simple arguments above.

Because the 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 γ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.

thumbnail 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 = Bez is given by

(28a)

(28b)

(28c)

(28d)

with the initial four-velocity component along B, u the initial four-velocity component perpendicular to B, and . 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 . Therefore, the only relevant parameters, apart from the initial conditions, are b = τmωB and . The length scale is therefore given in units of the non-relativistic Larmor radius

(29)

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

(30a)

(30b)

the solution reads

(31a)

(31b)

(31c)

(31d)

where the Ci 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.

thumbnail Fig. 3.

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

The corresponding Lorentz factor decreases according to Eq. (28a) and is shown in Fig. 4. The time when damping sets in is given approximately by . These times are shown as coloured vertical lines in Fig. 4. If the particle moves along the field line, it experiences no damping and keeps a uniform motion.

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

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

thumbnail 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 < 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 vE = E ∧ B/B2. Therefore, the velocity is βE = vE/c = E/cB and the corresponding final Lorentz factor .

We performed simulations with E/cB = 0.999 and initial Lorentz factors log γ0 = {0, 4, 8}. The final Lorentz factor is γ ≈ 22.3. Figure 7 shows the Lorentz factor (coloured dots) compared to the analytical expression (solid black lines). The agreement is excellent and demonstrates the high efficiency of our algorithm to capture ultra-relativistic motion with high precision.

thumbnail Fig. 7.

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

The quantitative agreement is checked by transforming the trajectory to the electric drift frame denoted by the coordinates (x′,y′). In this frame the orbital radius is decreasing, as shown in Fig. 8 for log γ0 = 4 and log b = −5, using different proper time steps such as log(ωB Δτ) = {−1, −2}. The analytical solution found from the appropriate parameters in Eq. (31) is overlapped in red solid lines.

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

thumbnail Fig. 9.

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

3.4. Parallel electric and magnetic field

Another interesting configuration not reducible to any of the previous configurations are the parallel electric and magnetic fields. In this case the second electromagnetic invariant does not vanish ℐ2 ≠ 0. Consequently, there is no frame where either the electric or magnetic field vanishes. The electric and magnetic velocity components uE and uB decouple into an acceleration along the common direction and a gyration in the same direction. Assuming this direction to be ez, the four-velocity reads

(32a)

(32b)

(32c)

(32d)

We recognise the special cases of a pure electric field for ut, uz and a pure magnetic field for ux, uy, the only difference being the value of , 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 γ0 = 8. The damping factor is set to log b = {0, −5, −10, −15} and the electric field strength E0 is varied relative to B0 such that log(E0/cB0) = {−2, 0, 2} (solid lines, dashed lines, and dotted lines, respectively). For a weak electric field E0 ≪ cB0 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 (ez). Increasing E0 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 E0 ≫ cB0 electric acceleration quickly sets in. Figure 11 shows some results of numerical simulations with a weak electric field log(E0/cB0) = − 2, pertinent for almost force-free neutron star magnetospheres, and initial Lorentz factors log γ0 = {0, 4, 8} and log b = {0, −5, −10, −15}.

thumbnail 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/cB0) = {−2, 0, 2} in solid lines, dashed lines, and dotted lines, respectively.

thumbnail 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/cB0) = − 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/cB = {0.1, 0.999}. The particle starts with an initial velocity in a plane perpendicular to B and Lorentz factor log γ0 = 4 in a field with log b = −5. Figure 12 shows the evolution to alignment of the velocity vector with the radiation reaction limit direction for a weak parallel electric field component (see the legend). We note that the angle θ should not be interpreted as the angle between the velocity vector and the magnetic field direction because the velocity in Eq. (17) is not along B. We observe that the time required for alignment is insensitive to the ratio E/cB, 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.

thumbnail Fig. 12.

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

3.6. Dipole magnetic field

As a step towards realistic configurations, we also investigate particle motion in a static magnetic dipole. Unfortunately, no simple exact analytical expressions are available for checking the algorithm; therefore, no quantitative accurate convergance test can be performed. First, we study the magnetic drift in the equatorial plane of a dipole field. Next, we look at trapped particles due to the mirror effect.

3.6.1. Magnetic drift

The magnetic drift motion in the equatorial plane of a dipole field is an interesting example to test our algorithm. The characteristic frequency is again ωB and the particle initial Lorentz factor is γ0 = 104. 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.

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

thumbnail 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).

thumbnail Fig. 15.

Decrease in the Lorentz factor of particles subject to the mirror effect.

thumbnail 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* = 108 T; millisecond pulsars with period P* = 5 ms and B* = 105 T; and magnetars with period P* = 10 s and B* = 1010 T. These quantities and their associated normalised strengths and damping parameters aB, aE, 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.

Table 1.

Typical period and surface magnetic field strength of millisecond pulsars, young pulsars, and magnetars.

We distinguish three kind of particles: a first group crashing onto the stellar surface, a second group of trapped particles, and a third group of escaping particles, all accelerated to high energies. Particles are considered trapped when they still have not crashed onto the surface or have not yet escaped the light cylinder. Particles are placed regularly within the light cylinder, starting at rest or with an initial velocity vector oriented along the magnetic field line, directed towards the star or towards infinity, with a Lorentz factor equal to γ0 = 103 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 r0: at the surface R*; approximately half-way between the surface and the light-cylinder (a geometric average); and at the light cylinder, thus . For comparison, we performed simulations with and without the radiation reaction.

4.2. Orders of magnitude

Before presenting the accurate numerical simulations of particle trajectories and their radiation reaction in the vicinity of neutron stars, we recall the orders of magnitude of the maximum Lorentz factors expected when charges are accelerated in the electric potential produced by a rotating magnetised perfectly conducting star. The most optimistic view adopts the full potential drop between the pole and the equator as an estimate of the accelerating field, and thus

(33)

where rB = 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

(34)

which is a factor R*/rL smaller than for the former case. Table 2 summarises the maximum Lorentz factors for electrons, protons, and iron around millisecond pulsars, young pulsars, and magnetars. The values reported in this table for are at best upper limits for the vacuum case. Only an accurate numerical integration of the equation of motion gives robust results, as we now show.

Table 2.

Maximum Lorentz factor orders of magnitude from conservative arguments about neutron star magnetospheres.

4.3. Escaping particles

Particles reaching distances larger than 10 rL 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 ∼ 109 in the LLR approximation for millisecond pulsars. This is, however, two orders of magnitude less than without the radiation reaction where γf ∼ 1011. Second, as expected, protons and iron nuclei acquire much less energy, only about γf ∼ 106 for millisecond pulsars, whether LLR is used or not. For young pulsars, electrons also reach γf ∼ 109 in the LLR regime instead of γf ∼ 1011 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 ∼ 104 − 104.5. For magnetars, the radiation reaction remains negligible irrespective of the nature of each species. Electrons reach energies up to γf ∼ 107.5, whereas protons and iron nuclei γf ∼ 103 − 103.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.

thumbnail 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 ∼ 1011 for magnetars, but only γf ∼ 108.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 ∼ 107.5 − 108.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.

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

thumbnail 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

(35)

where aBL 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 ≈ 109.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 109.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 aBL and the Lorentz factor such that

(36)

Table 3.

Maximum Lorentz factors γmax for the three kinds of particles: electrons / protons / iron nuclei.

This increase in the acceleration efficiency is imputed to the presence of a still strong radial component of the electric field, which was absent in the study of Pétri (2021). The simulations performed in this section only followed a small number of particles, due to the stringent computation time required to accurately evolve the particle velocity and position. Describing the plasma feedback onto the electromagnetic field would require a much larger number of particles coupled to the evolution of the electromagnetic field via Maxwell equations, leading to a PIC code. To date, although PIC codes exist and have been adapted to simulate neutron star magnetospheres, none has yet been able to handle the parameter space explored in the present work. Therefore, let us contrast our results in the light of the existing kinetic descriptions of the magnetosphere.

4.7. Comparison to previous works

Several investigations of particle acceleration and radiation reaction around neutron stars have been attempted in the literature. However, due to severe numerical limitations, studies employing realistic field strengths for the neutron star are very rare. We mention, however, the pioneering work of Finkbeiner et al. (1989) and Finkbeiner et al. (1990), who employed a single test particle approach with the radiation reaction and found acceleration around the neutron star up to Lorentz factors of about γf ≲ 109 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 ≲ 109. 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 ≲ 103 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 ≲ 102. 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 ≲ 103 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 108 and slightly larger for protons (almost 109), but for faster rotations in a field of an oblique rotating dipole with strengths of 5 × 106 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 109 and for protons of about 106.

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 γ0 = 103 and a velocity parallel to v± in order to have comparable initial conditions for both sets of runs.

Figure 20 shows a sample of electron trajectories and demonstrates the reasonable results obtained by this asymptotic regime for a millisecond pulsar. However, the precision depends on the particle initial position. For motions starting at the surface (upper row in Fig. 20) some trajectories are well reproduced by the radiation reaction limit regime. The accuracy is less good for the brown and green paths, although the general trend is conserved. When starting at larger distances from the surface, for example at (middle row of Fig. 20), we observe better agreement between the two regimes. The best results are obtained for particles well away from the surface, starting at r = rL (bottom row of Fig. 20). All trajectories computed in the radiation reaction limit regime overlap with the LLR integration.

thumbnail Fig. 20.

Sample of trajectories for electrons obtained with LLR in the field of a millisecond pulsar (solid lines) and in the radiation reaction limit (dotted symbols). The trajectories are projected along the xy plane in the left panel and on the xz plane in the right panel. Particles are launched from the stellar surface in the top row, at a distance in the middle row, and at a distance r = 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 = rL 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.

thumbnail 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 r0/rL in the legend. For reference the Δt2 behaviour is shown as blue filled circles. In the previous simulations we fixed the time step to Δt ≈ 10−4, and thus we expect a precision better than three digits in all cases. The discrepancy between the Landau–Lifshitz and radiation reaction regimes can therefore not be explained by a discretisation effect. We also checked that the initial condition of the velocity does not impact the trajectory in Landau–Lifshitz. The explanation must be searched for in the deficiency of the radiation reaction regime to satisfactorily account for all possible trajectories. This regime assumes a radiative friction force opposite to the 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

(37a)

(37b)

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

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

  1. Abraham, M. 1902, Ann. Phys., 315, 105 [NASA ADS] [CrossRef] [Google Scholar]
  2. Abraham, M. 1904, Ann. Phys., 319, 236 [NASA ADS] [CrossRef] [Google Scholar]
  3. Boghosian, B. M. 1987, PhD Thesis, University of California [Google Scholar]
  4. Boris, J. 1970, Proceeding of Fourth Conference on Numerical Simulations of Plasmas [Google Scholar]
  5. Cerutti, B., Philippov, A., Parfrey, K., & Spitkovsky, A. 2015, MNRAS, 448, 606 [Google Scholar]
  6. Cerutti, B., Philippov, A. A., & Spitkovsky, A. 2016, MNRAS, 457, 2401 [Google Scholar]
  7. Cerutti, B., Philippov, A. A., & Dubus, G. 2020, A&A, 642, A204 [Google Scholar]
  8. Deutsch, A. J. 1955, Ann. d’Astrophys., 18, 1 [Google Scholar]
  9. Dirac, P. A. M. 1938, Proc. R. Soc. Lond. Ser. A, 167, 148 [NASA ADS] [CrossRef] [Google Scholar]
  10. Elkina, N. V., Fedotov, A. M., Herzing, C., & Ruhl, H. 2014, Phys. Rev. E, 89, 053315 [NASA ADS] [CrossRef] [Google Scholar]
  11. Ferrari, A., & Trussoni, E. 1974, A&A, 36, 267 [NASA ADS] [Google Scholar]
  12. Finkbeiner, B., Herold, H., Ertl, T., & Ruder, H. 1989, A&A, 225, 479 [NASA ADS] [Google Scholar]
  13. Finkbeiner, B., Herold, H., & Ruder, H. 1990, A&A, 238, 462 [NASA ADS] [Google Scholar]
  14. Fulton, T., & Rohrlich, F. 1960, Ann. Phys., 9, 499 [NASA ADS] [CrossRef] [Google Scholar]
  15. Gordon, D. F., & Hafizi, B. 2021, Comput. Phys. Commun., 258, 107628 [NASA ADS] [CrossRef] [Google Scholar]
  16. Gordon, D. F., Hafizi, B., & Palastro, J. 2017a, AIP Conf. Proc., 1812, 050002 [CrossRef] [Google Scholar]
  17. Gordon, D. F., Palastro, J. P., & Hafizi, B. 2017b, Phys. Rev. A, 95, 033403 [NASA ADS] [CrossRef] [Google Scholar]
  18. Guépin, C., Cerutti, B., & Kotera, K. 2020, A&A, 635, A138 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  19. Hadad, Y., Labun, L., Rafelski, J., et al. 2010, Phys. Rev. D, 82, 096012 [NASA ADS] [CrossRef] [Google Scholar]
  20. Heintzmann, H., & Schrüfer, E. 1973, Phys. Lett. A, 43, 287 [NASA ADS] [CrossRef] [Google Scholar]
  21. Herold, H., Ertl, T., & Ruder, H. 1985, Mitteilungen der Astronomischen Gesellschaft Hamburg, 63, 174 [Google Scholar]
  22. Kalapotharakos, C., Brambilla, G., Timokhin, A., Harding, A. K., & Kazanas, D. 2018, ApJ, 857, 44 [Google Scholar]
  23. Kelner, S. R., Prosekin, A. Y., & Aharonian, F. A. 2015, AJ, 149, 33 [Google Scholar]
  24. Kulsrud, R. M. 1972, ApJ, 174, L25 [NASA ADS] [CrossRef] [Google Scholar]
  25. Landau, L., & Lifchitz, E. 1989, Physique théorique: Tome 2, Théorie des champs (Moscou: Mir) [Google Scholar]
  26. Laue, H., & Thielheim, K. O. 1986, ApJS, 61, 465 [NASA ADS] [CrossRef] [Google Scholar]
  27. Li, F., Decyk, V. K., Miller, K. G., et al. 2021, J. Comput. Phys., 438, 110367 [NASA ADS] [CrossRef] [Google Scholar]
  28. Lorentz, H. A. H. A. 1916, The Theory of Electrons and its Applications to the Phenomena of Light and Radiant Heat (Leipzig: G.E. Stechert) [Google Scholar]
  29. Mestel, L., Robertson, J. A., Wang, Y. M., & Westfold, K. C. 1985, MNRAS, 217, 443 [Google Scholar]
  30. Michel, F., & Li, H. 1999, Phys. Rep., 318, 227 [NASA ADS] [CrossRef] [Google Scholar]
  31. Olver, F. W. J. 2010, NIST Handbook of Mathematical Functions (Cambridge: Cambridge University Press) [Google Scholar]
  32. Philippov, A. A., & Spitkovsky, A. 2018, ApJ, 855, 94 [NASA ADS] [CrossRef] [Google Scholar]
  33. Philippov, A. A., Spitkovsky, A., & Cerutti, B. 2015, ApJ, 801, L19 [NASA ADS] [CrossRef] [Google Scholar]
  34. Piazza, A. D. 2008, Lett. Math. Phys., 83, 305 [NASA ADS] [CrossRef] [Google Scholar]
  35. Pétri, J. 2020, J. Plasma Phys., 86, 825860402 [CrossRef] [Google Scholar]
  36. Pétri, J. 2021, MNRAS, 503, 2123 [Google Scholar]
  37. Rohrlich, F. 2007, Classical Charged Particles, 3rd edn. (Singapore: World Scientific Pub Co Inc) [CrossRef] [Google Scholar]
  38. Tomczak, I., & Pétri, J. 2020, J. Plasma Phys., 86, 825860401 [NASA ADS] [CrossRef] [Google Scholar]
  39. Vranic, M., Martins, J. L., Fonseca, R. A., & Silva, L. O. 2016, Comput. Phys. Commun., 204, 141 [NASA ADS] [CrossRef] [Google Scholar]

All Tables

Table 1.

Typical period and surface magnetic field strength of millisecond pulsars, young pulsars, and magnetars.

Table 2.

Maximum Lorentz factor orders of magnitude from conservative arguments about neutron star magnetospheres.

Table 3.

Maximum Lorentz factors γmax for the three kinds of particles: electrons / protons / iron nuclei.

All Figures

thumbnail 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
thumbnail 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
thumbnail 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
thumbnail 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
thumbnail 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
thumbnail 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
thumbnail 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
thumbnail 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
thumbnail Fig. 9.

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

In the text
thumbnail 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/cB0) = {−2, 0, 2} in solid lines, dashed lines, and dotted lines, respectively.

In the text
thumbnail 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/cB0) = − 2, and log b = {0, −5, −10, −15}. The black solid lines correspond to the exact analytical solutions.

In the text
thumbnail 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
thumbnail 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
thumbnail 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
thumbnail Fig. 15.

Decrease in the Lorentz factor of particles subject to the mirror effect.

In the text
thumbnail 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
thumbnail 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
thumbnail Fig. 18.

Same as Fig. 17, but for crashed particles.

In the text
thumbnail Fig. 19.

Same as Fig. 17, but for trapped particles.

In the text
thumbnail Fig. 20.

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

In the text
thumbnail Fig. 21.

Same as Fig. 20, but for protons.

In the text
thumbnail 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
thumbnail Fig. 23.

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

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.