© The Authors 2023

1. Introduction

Neutron stars are compact stellar remnants left after supernovae explosions. Due to their intense magnetic and electric fields, they are believed to be efficient sources of ultra-relativistic particles. They also emit photons via the synchrotron and curvature radiation mechanisms or inverse Compton, interacting with their surrounding as well as with the interstellar medium.

In this paper, we use a simple model of a neutron star described by only a few parameters, namely the inclination of the neutron star χ, corresponding to the angle between its rotation and magnetic axes, the angular rotation speed of the neutron star Ω, allowing to define the light cylinder radius r_L = c/Ω, which is the distance at which an object in co-rotation with the neutron star would reach the speed of light c, the radius of the neutron star R and the magnetic field strength B at the surface of the neutron star.

The extreme magnetic field of these remnant stars ranges from B ≃ 10⁵ T to B ≃ 10¹⁰ T. Moreover, coupled with an angular speed between Ω ≃ 1000 rad s⁻¹ and Ω ≃ 0.1 rad s⁻¹, they generate an intense electric field E, accelerating particles around the neutron star to ultra-relativistic Lorentz factors. For our purposes, the neutron star mass M is irrelevant because the gravitational force exerted on charged particles F_g is negligible compared to the Lorentz force F_L. Indeed, both forces have a typical intensity of $F_{\mathrm{g}}=G\frac{M{m}_{\mathrm{p}}}{R^2}\simeq 3.17\times {10}^{-15}\text{N}$ for a proton at the surface of the neutron star compared to F_L = q E = q ΩR B ≃ 5.01 × 10⁻⁷ N. The ratio of these forces is

$$\begin{array}{c}\hfill \frac{F_{\mathrm{g}}}{F_{\mathrm{L}}}=\frac{GM{m}_{\mathrm{p}}}{R^{2}qE}=6.33\times {10}^{-9},\end{array}$$

so the electromagnetic force is approximately 10⁹ times stronger than the gravitational force for protons and even a factor m_p/m_e ∼ 2000 larger for electrons and positrons.

Most numerical simulations of neutron star magnetospheres use the Boris (1970) scheme or better the Vay (2008) scheme; however, these algorithms are not well suited for ultra-strong electromagnetic fields, and some authors have consequently lowered the true field strengths in their simulations to unrealistically low values. Although scaling is sometimes applied to obtain results closer to reality, such scaling cannot be straightforwardly extrapolated, for instance, when radiation reaction is included because of the non-linearities introduced by radiative feedback (Vranic et al. 2016). Moreover, Lorentz factors reached by the particles near pulsars hardly exceed γ = 10⁴ with those algorithms (see, for instance, Brambilla et al. 2018; Philippov & Spitkovsky 2018; Guépin et al. 2020, and Kalapotharakos et al. 2018).

This limitation arises from the huge span in timescales, starting from the ultra-high frequency gyro motion ω_B down to the stellar rotation frequency Ω. These extreme values are synthesised by their ratio, which is also known as the strength parameter,

$$\begin{array}{c}\hfill a=\frac{{\omega }_B}{\mathrm{\Omega }}=\frac{qB}{m\mathrm{\Omega }}\simeq {10}^{10}\end{array}$$

for a proton near a millisecond pulsar (Ω = 10³ rad s⁻¹ and B = 10⁵ T). The situation becomes worse for electron-positron pairs and for young pulsars. This ratio corresponds to the number of gyrations made by a particle during the timescale of evolution of the electromagnetic field due to the stellar rotation. It shows that the difference in timescales makes computing the trajectory of particles in a reasonable amount of time almost impossible since billions of time steps are needed in the pulsar period timescale.

To tackle this issue, our aim is to propose a new technique based on analytical solutions of the Lorentz force equation in constant and ultra-strong electromagnetic fields with an acceptable computational time and, most importantly, an approach that avoids the need for scaling, allowing particles to reach high Lorentz factors with realistic fields. Several authors have worked on analytical solutions to the equation of motion including the radiation reaction, such as Gordon & Hafizi (2021), Laue & Thielheim (1986), Li et al. (2021), and Heintzmann & Schrüfer (1973). Our approach is based on these works. For plane waves in a vacuum, exact solutions are also known because of Hadad et al. (2010), Piazza (2008). These solutions were applied in strong electromagnetic waves by Pétri (2021) and around a dipole by Pétri (2022).

Neutron stars are known to act as unipolar inductors, generating huge electric potential drops between the poles and the equator of the order

$$\begin{array}{c}\hfill \mathrm{\Delta }\varphi =\mathrm{\Omega }B{R}^2\approx {10}^{16}V.\end{array}$$(1)

As a consequence, these stars expel electrons (and maybe protons and ions), filling the magnetosphere with charged particles. The typical Lorentz factor for electrons in this static field is therefore

$$\begin{array}{c}\hfill \gamma =\frac{e\mathrm{\Delta }\varphi }{m_{\mathrm{e}}{c}^2}\approx {10}^{10}.\end{array}$$(2)

If the particle injection rate is high enough, this plasma will screen the electric field, drastically mitigating the potential drop Δϕ and the acceleration efficiency. Resistive (Li et al. 2012) as well as PIC simulations (Cerutti et al. 2015) have indeed showed that only a small fraction of the full potential is available. However, for low particle injection rates, the plasma is unable to screen the electric field, and the full potential drop develops. In these cases, the magnetosphere is almost empty and known as an electrosphere (Krause-Polstorff & Michel 1985; Pétri et al. 2002). Such electrospheres are the subject of the present paper. They represent inactive pulsars that are able to accelerate particles to ultra-relativistic speeds. Our aim is to accurately quantify the final Lorentz factor reached by the outflowing plasma in this large-amplitude low-frequency electromagnetic wave. A similar study was performed by Michel & Li (1999), though with a more analytical perspective.

The outline of the paper is as follows. First in Sect. 2, we summarise the principle of the algorithm. Next in Sect. 3, we show some results obtained in fields where an analytical solution is known before discussing the results of the simulations near pulsars in Sect. 4. Some conclusions are drawn in Sect. 5.

2. Description of the numerical algorithm

In this section, we describe the algorithm developed including radiation reaction. It is similar to the one presented in Pétri (2020), which was used the basis for the work of Tomczak & Pétri (2020). This code successively finds analytical solutions to the equation of the motion of particles in an electromagnetic field (E, B) assumed to remain constant within a time step integration. During this time step, the algorithm solves the equation of motion, approximated by the reduced Landau-Lifshitz equation

$$\begin{array}{c}\hfill \frac{\mathrm{d}{u}^{\mu }}{\mathrm{d}\tau }=\frac{q}{m}{F^{\mu }}_{\nu }{u}^{\nu }-\frac{q^4}{6\pi {\epsilon }_0{m}^3{c}^3}[F^{\mu \sigma }{F}_{\lambda \sigma }{u}^{\lambda }+(F^{\sigma \nu }{u}_{\nu }{F}_{\sigma \lambda }{u}^{\lambda })u^{\mu }].\end{array}$$(3)

We note that the four-velocity of the particle (in contravariant form) is u^μ = (u⁰   u¹   u²   u³) = γc(1   β^x   β^y   β^z) = γc(1   β), with the Lorentz factor $\gamma =\frac{1}{\sqrt{1-{\beta }^2}}$; β is the speed of the particle normalised to the speed of light c, and it can be decomposed onto a Cartesian coordinate basis with the usual labels of x, y and z; and ${F^{\mu }}_{\nu }$ is the electromagnetic tensor given in components by

$$\begin{array}{c}\hfill F{{}^{\mu }}_{\nu }=\left(\begin{array}{cccc}0& E_x/c& E_y/c& E_z/c\\ E_x/c& 0& B_z& -B_y\\ E_y/c& -B_z& 0& B_x\\ E_z/c& B_y& -B_x& 0\end{array}\right).\end{array}$$(4)

For this algorithm, we first needed to switch from the observer’s reference frame denoted by ℛ, in which the neutron star only rotates, to a reference frame where E and B are parallel and denoted by ℛ′, where integrating the Lorentz force is easy. In a Cartesian coordinate system, the evolution of the four-velocity when E and B are constant and along the z-axis is

$$\begin{array}{cc}& u^0/c=\gamma =a(\tau ){\gamma }_0[cosh({\omega }_E\tau )+{\beta }_0^{z}sinh({\omega }_E\tau )]\hfill \end{array}$$(5a)

$$\begin{array}{cc}& u^1/c=\gamma {\beta }^x=b(\tau ){\gamma }_0[{\beta }_0^{x}cos({\omega }_B\tau )+{\beta }_0^{y}sin({\omega }_B\tau )]\hfill \end{array}$$(5b)

$$\begin{array}{cc}& u^2/c=\gamma {\beta }^y=b(\tau ){\gamma }_0[-{\beta }_0^{x}sin({\omega }_B\tau )+{\beta }_0^{y}cos({\omega }_B\tau )]\hfill \end{array}$$(5c)

$$\begin{array}{cc}& u^3/c=\gamma {\beta }^z=a(\tau ){\gamma }_0[sinh({\omega }_E\tau )+{\beta }_0^{z}cosh({\omega }_E\tau )].\hfill \end{array}$$(5d)

We introduced the typical electric and magnetic frequencies as ω_E = q E/mc and ω_B = q B/m. The initial Lorentz factor is indicated with γ₀, and $({\beta }_0^x{\beta }_0^y{\beta }_0^z)$ is the initial velocity normalised to the speed of light c. The coefficients a(τ) and b(τ) bring corrections to the velocity induced by the radiation reaction (setting a(τ) = b(τ) = 1 removes radiation reaction). The coefficients are given by

$$\begin{array}{cc}& a(\tau )=\frac{1}{{\gamma }_0\sqrt{(1-{\beta }_z^2)-({\beta }_x^2+{\beta }_y^2)exp(-2{\tau }_0({\omega }_B^2+{\omega }_E^2)\tau )}}\hfill \end{array}$$(6a)

$$\begin{array}{cc}& b(\tau )=a(\tau )exp(-{\tau }_0({\omega }_B^2+{\omega }_E^2)\tau )\hfill \end{array}$$(6b)

with τ₀ = q²/6πε₀mc³, a characteristic timescale of the energy losses. It is unfortunately not possible to integrate the four-position in an analytical manner, so we decided to use the analytical solution for the position without the radiation reaction

$$\begin{array}{cc}& c(t-t_0)=\frac{{\gamma }_{0}c}{{\omega }_E}[sinh({\omega }_E\tau )+{\beta }_0^{z}cosh({\omega }_E\tau )-{\beta }_0^z]\hfill \end{array}$$(7a)

$$\begin{array}{cc}& x-x_0=\frac{{\gamma }_{0}c}{{\omega }_B}[{\beta }_0^{x}sin({\omega }_B\tau )-{\beta }_0^{y}cos({\omega }_B\tau )+{\beta }_0^y]\hfill \end{array}$$(7b)

$$\begin{array}{cc}& y-y_0=\frac{{\gamma }_{0}c}{{\omega }_B}[{\beta }_0^{x}cos({\omega }_B\tau )-{\beta }_0^x+{\beta }_0^{y}sin({\omega }_B\tau )]\hfill \end{array}$$(7c)

$$\begin{array}{cc}& z-z_0=\frac{{\gamma }_{0}c}{{\omega }_E}[cosh({\omega }_E\tau )-1+{\beta }_0^{z}sinh({\omega }_E\tau )],\hfill \end{array}$$(7d)

with the initial position and time being (x₀, y₀, z₀) and t₀. We integrated the position of the particle according to τ, the proper time; however, as we wanted to set the observer’s time step δt = t − t₀ to be constant, we resorted to finding dτ, the proper time step, so as to always obtain the same observer time step dt.

This scheme has some drawbacks: It returns an approximation of the four-position of the particle with radiation reaction, and it does not efficiently take into account field gradients. However, this method is still more efficient than other approximations, such as Euler or Runge-Kutta, most of which apply a constant speed assumption. Indeed, this scheme does not assume a constant speed but only that the radiation reaction has little effect on the position of the particle. Nonetheless, it still takes into account speed variations in terms of norm and direction. In addition, compared to other methods, it allows for longer time steps and multiple gyrations of the particle around a magnetic field line during those longer time steps. Moreover, since the gyration radius is relatively small compared to typical magnetic field scales (followed by the trajectories of the particles), the error related to the shrinking of the Larmor radius is negligible.

The light-like case, where E ⋅ B = 0 and E² = c² B², must be treated separately because there exists no frame where E and B are parallel. Due to the low likelihood of finding such configurations, we kept solutions found in Pétri (2020) as not corrected for radiation reaction. This means that if a particle were to find a light-like field, the algorithm would keep working but the radiation reaction would not be taken into account for just one time step.

3. Tests

In order to check the correctness and accuracy of our code implementation, we simulated the evolution of a particle in a constant and uniform magnetic field. Analytical solutions are known for the four-velocity but also for the spatial position and observer time in this case (Pétri 2022). An electron is kicked into a constant magnetic field of strength B_z with an initial Lorentz factor of γ₀ = 10⁴. The normalised damping parameter is therefore τ₀ ω_B = 10⁻⁵. As Fig. 1 highlights, the algorithm converges and finds the spiral motion described by a charged particle losing energy in a magnetic field. The exact analytical trajectory is shown in orange solid lines and the numerical simulations with blue dots. Thanks to this comparison, we observed that the algorithm is first order in proper time, according to the particle position. The first order convergence is due to the fact that the four-position was updated according to the pure Lorentz force, neglecting the radiation reaction corrections.

Fig. 1. Particle trajectory in a constant magnetic field with radiation reaction. Panel a: trajectory of the simulated particle (in blue) compared to theoretical positions (in orange) for τ0 ωB = 10−5, γ0 = 104, initial speed along y. Panel b: relative error as a function of the proper time step dτ for the simulations in blue points, compared to the first order error expectations in red solid line.

Being confident about the convergence and accuracy of our algorithm, we then simulated the motion of charged particles in the electromagnetic field of a rotating neutron star.

4. Simulations in the Deutsch field

As an application in the astrophysical context, we explored particle acceleration and its radiation reaction in a strongly magnetised rotating magnetic dipole such as that expected around rotating neutron stars. If the star is surrounded by vacuum, simple analytical expressions are known for the electromagnetic field and given by Deutsch (1955).

4.1. Neutron star settings

Therefore, we injected particles in the field of a rotating neutron star in a vacuum, namely, the Deutsch field (Deutsch 1955), decomposed in spherical coordinates and using the complex form for the magnetic field as

$$\begin{array}{cc}& B_r(\mathit{r},t)=2B[\frac{R^3}{r^3}cos\chi cos\vartheta +\frac{R}{r}\frac{h_1^{(1)}(kr)}{h_1^{(1)}(kR)}sin\chi sin\vartheta e^{i\psi }]\hfill \end{array}$$(8a)

$$\begin{array}{cc}& B_{\vartheta }(\mathit{r},t)=B[\frac{R^3}{r^3}cos\chi sin\vartheta \hfill \\ \hfill & +(\frac{R}{r}\frac{\frac{\mathrm{d}}{\mathrm{d}r}\left(r{h}_1^{(1)}(kr)\right)}{h_1^{(1)}(kR)}+\frac{R^2}{r_{\mathrm{L}}^2}\frac{h_2^{(1)}(kr)}{\frac{\mathrm{d}}{\mathrm{d}r}\left(r{h}_2^{(1)}(kr)\right){|}_R})sin\chi cos\vartheta e^{i\psi }]\hfill \end{array}$$(8b)

$$\begin{array}{cc}& B_{\phi }(\mathit{r},t)=B[\frac{R}{r}\frac{\frac{\mathrm{d}}{\mathrm{d}r}(r{h}_1^{(1)}(kr))}{h_1^{(1)}(kR)}+\frac{R^2}{r_{\mathrm{L}}^2}\frac{h_2^{(1)}(kr)}{\frac{\mathrm{d}}{\mathrm{d}r}\left(r{h}_2^{(1)}(kr)\right){|}_R}cos2\vartheta ]isin\chi e^{i\psi }\hfill \end{array}$$(8c)

and for the electric field as

$$\begin{array}{cc}& E_r(\mathit{r},t)=\mathrm{\Omega }BR[(\frac{2}{3}-\frac{R^2}{r^2}(3{cos}^2\vartheta -1))\frac{R^2}{r^2}cos\chi \hfill \\ \hfill & +3sin\chi sin2\vartheta e^{i\psi }\frac{R}{r}\frac{h_2^{(1)}(kr)}{\frac{\mathrm{d}}{\mathrm{d}r}\left(r{h}_2^{(1)}(kr)\right){|}_R}]\hfill \end{array}$$(9a)

$$\begin{array}{cc}& E_{\vartheta }(\mathit{r},t)=\mathrm{\Omega }BR[-\frac{R^4}{r^4}sin2\vartheta cos\chi \hfill \\ \hfill & +sin\chi e^{i\psi }(\frac{R}{r}\frac{\frac{\mathrm{d}}{\mathrm{d}r}\left(r{h}_2^{(1)}(kr)\right)}{\frac{\mathrm{d}}{\mathrm{d}r}\left(r{h}_2^{(1)}(kr)\right){|}_R}cos2\vartheta -\frac{h_1^{(1)}(kr)}{h_1^{(1)}(kR)})]\hfill \end{array}$$(9b)

$$\begin{array}{cc}& E_{\phi }(\mathit{r},t)=\mathrm{\Omega }BR[\frac{R}{r}\frac{\frac{\mathrm{d}}{\mathrm{d}r}\left(r{h}_2^{(1)}(kr)\right)}{\frac{\mathrm{d}}{\mathrm{d}r}\left(r{h}_2^{(1)}(kr)\right){|}_R}-\frac{h_1^{(1)}(kr)}{h_1^{(1)}(kR)}]isin\chi cos\vartheta e^{i\psi },\hfill \end{array}$$(9c)

where $h_{\mathcal{l}}^{(1)}$ represents the spherical Hankel functions for outgoing waves (Arfken & Weber 2005), k r_L = 1, and the phase is ψ = φ − Ω t, with φ as the phase at t = 0. The physical components of the electromagnetic field are the real parts of the above expressions.

In our simulations, the magnetic field at the surface of the star was set to B = 10⁵ T, which corresponds to the typical values for a millisecond pulsar. We worked in normalised units, setting the reference rotation speed of the neutron star to Ω, the reference speed to the speed of light c, and therefore the light cylinder radius to r_L = c/Ω. By choosing the normalised neutron star radius to be R = 0.1 r_L, we fixed the size of the light cylinder and the stellar rotation speed. Indeed, the neutron star radius is about R = 12 km (Nättilä et al. 2017), which means that the light cylinder radius is about r_L = 120 km. Moreover, since r_L = c/Ω, the angular velocity is Ω = 2500 rad s⁻¹. Thus, a real period of rotation of P is 2.5 ms. We needed the observer’s time step to be small relative to the time of evolution of the magnetic field (which is of the order of the rotation period), so we chose dt/P = 0.0001, meaning dt = 250 ns. Concerning the inclination of the neutron star, the simulations were carried out with an inclination from the following set: χ ∈ {0 ° ;30 ° ;60 ° ;90 ° ;120 ° ;150 ° ;180 ° }.

Particle injection was made using a rejection method in order to obtain a uniform and isotropic distribution of particles around the neutron star. We generated three random numbers, each following an independent uniform distribution law in the interval [ − 0.9; 0.9] and corresponding to the Cartesian coordinates (x, y, z). We then defined the radius at which a particle is injected, $r=\sqrt{x^2+{\mathit{y}}^2+z^2}$, and if r/r_L ≤ 0.1 or r/r_L ≥ 0.9, we removed that particle and generated it again. Otherwise, we kept the particle and injected the next one.

Particles were injected at rest and evolved in time up to a final time of t_f/P = 15. Sometimes particles crashed onto the surface, and the integration was stopped earlier. As for particle species, we considered electrons and protons. In order to obtain reasonable statistics, we simulated 8192 particles per configuration. In the following sections, we describe the final particle properties, including their distribution in space and energy.

4.2. Distribution of particles in space

In this sub-section, we discuss the particle positions at the end of the run as well as at the beginning of the run in order to link them to their final Lorentz factor. We note, however, that the cases of aligned (χ = 0°) and anti-aligned (χ = 180°) rotators are treated separately in Sect. 4.4 because the electromagnetic field is static in these configurations.

The coordinates used to characterise the position of a particle were either in the Cartesian coordinate system (x, y, z) or the spherical coordinate system (r, θ, ϕ). In the spherical coordinate system, r defines the distance of the particle relative to the centre of the neutron star, θ is its colatitude (relative to the rotation axis), and ϕ is its azimuth relative to the x axis, knowing that at the beginning of the simulation the magnetic axis lies in the xOz plane. We also defined three possible final states for the particles: ejected, meaning that at the end r/r_L ≥ 1 (in our case those particles have a radial velocity); trapped, meaning that at the end r/r_L ∈ ]0.1; 1[; or crashed onto the neutron star, meaning that the particles should have reached r/r_L ≤ 0.1 or equivalently r < R at one point.

The statistics of the particles according to their final states are summarised in Table 1 for our total number of 8192 particles. It allowed us to notice that for the protons, as χ increases from 30° to 150°, fewer particles impact the surface of the neutron star, while, inversely, more protons are ejected away from it. We also found it interesting to notice that the maximum number of protons trapped close to the pulsar was obtained for χ = 90° and that a symmetrical behaviour could be observed for electrons, namely, as χ increases, more particles impact the surface of the neutron star, fewer electrons are ejected away from it, and we still find that for χ = 90°, most electrons are trapped close to the neutron star in the same proportion as the protons. We however noticed a few differences: Protons for an inclination χ are ejected more easily than electrons for an inclination π − χ, which either become trapped or crash onto the surface more frequently.

Actually, we noticed that when respectively comparing the trajectories of protons and electrons in Figs. 2–9, we found the protons and electrons to possess very similar trajectories. If a proton starting at a position (r₀; θ₀; ϕ₀) ends at position (r_f; θ_f; ϕ_f) for a pulsar of inclination χ, an electron close to a pulsar of inclination π − χ starting at position (r₀; π − θ₀; ϕ₀) is very likely to end the simulation at position (r_f; π − θ_f; ϕ_f).

Fig. 2. Final state of the protons depending on their initial positions around the pulsar. The obliquity and the initial position are χ = 60°, r0 ∈ [0.3; 0.4] for (a) and r0 ∈ [0.8; 0.9] for (b), χ = 90°, r0 ∈ [0.3; 0.4] for (c) and r0 ∈ [0.8; 0.9] for (d), radiation reaction being enabled.

Fig. 3. Final state of the electrons depending on their initial positions around the pulsar. The obliquity and the initial position are χ = 120°, r0 ∈ [0.3; 0.4] for (a) and r0 ∈ [0.8; 0.9] for (b), χ = 90°, r0 ∈ [0.3; 0.4] for (c) and r0 ∈ [0.8; 0.9] for (d), radiation reaction being enabled.

Fig. 4. Map of the impact spots and Lorentz factor of protons on the surface for inclination χ = 30° (a), 60° (b), 90° (c), and 120° (d). Radiation reaction was enabled, and the magnetic axis is in the ϕ = 0 plane.

Fig. 5. Map of the impact spots and Lorentz factor of electrons on the surface for inclination χ = 150° (a), 120° (b), 90° (c), and 60° (d). Radiation reaction was enabled, and the magnetic axis is in the ϕ = 0 plane.

Fig. 6. Map of the final colatitude, azimuth, and radius (colour) of protons around neutron stars of inclination 60° (a), 90° (b), 120° (c), and 150° (d). Radiation reaction was enabled, and the magnetic axis is in the ϕ = 0 plane.

Fig. 7. Map of the final colatitude, azimuth, and radius (colour) of electrons around neutron stars of inclination 60° (a), 90° (b), 120° (c), and 150° (d). Radiation reaction was enabled, and the magnetic axis is in the ϕ = 0 plane.

Fig. 8. Map of the ejection colatitude and azimuth and Lorentz factor of protons on neutron stars of inclination 30° (a), 60° (b), 90° (c), and 120° (d). Radiation reaction was enabled.

Fig. 9. Map of the ejection colatitude and azimuth and Lorentz factor of electrons on neutron stars of inclination 150° (a), 120° (b), 90° (c), and 60° (d). Radiation reaction was enabled.

The symmetrical behaviour found in these simulations reflects the fact that the trajectories are relatively insensitive to the charge over mass ratio q/m of the particles, except for the sign of the charge itself. Indeed, in the ultra-relativistic regime, the mass of the particles becomes negligible compared to their total kinetic energy, and they can be considered as massless particles just like photons, for instance. However, because the radiation reaction force does not scale linearly with this ratio q/m, it is not at all obvious that the trajectories will remain similar. Nevertheless, we found that the radiation reaction impacts the motion of protons similarly to that of electrons. However, due to the difference in mass between protons and electrons, their respective Lorentz factors, although both ultra relativistic, scale like their mass ratios m_e/m_p.

Moreover, the statistics presented in Table 1 can be linked to the initial positions of the particles.

Starting positions. Indeed, the initial position of the particles has an influence on their final state. Figures 2 and 3 respectively show the map of the initial positions of protons and electrons, and their final state is indicated with a colour code. We chose two altitude intervals. The first is close to the surface, with r/r_L ∈ [0.3, 0.4] and χ = 60°, and the second is close to the light cylinder, with r/r_L ∈ [0.8, 0.9] and χ = 90°. The figures also highlight the importance of the neutron star obliquity χ. In addition, we note the figures show regions with clear boundaries and almost no overlap that are prone to ejection, crashing, or trapping of particles, meaning that the particle’s behaviour is well defined according to their starting point.

Comparing Figs. 2–10 allowed us to find either clear distinctions between the Lorentz factor of ejected particles and the other populations or that their speed would be similar no matter if they are ejected or not. Figure 10 also shows the influence of the initial radius on the spread of energy of the particles: as the initial radius rises, proton energies become less spread for χ = 60°, spanning six orders of magnitude for low altitude (r ∈ [0.3, 0.4]) and only less than two orders of magnitude for high altitude (r ∈ [0.8, 0.9]). However, for some other configurations, such as the orthogonal rotator with χ = 90°, the final Lorentz factor is not affected by the initial altitude when the particles are trapped.

Fig. 10. Final Lorentz factor of the particles shown in Fig. 2.

Upon comparing the results with the radiation reaction in Figs. 2 and 10 to those without radiation reaction in Fig. 11, we noticed that the radiation reaction drastically influences the behaviour of the particles regarding their Lorentz factors. The radiation reaction decreases the final Lorentz factor by at least one order of magnitude and sometimes changes the final state of the particles, depending on their initial position.

Fig. 11. Final Lorentz factor (a) and state (b) of the particles near a neutron star of inclination χ = 60°, with r0 ∈ [0.8; 0.9], and final Lorentz factor (c) and state (d) of the particles near a neutron star of inclination χ = 90°, with r0 ∈ [0.8; 0.9]. Radiation reaction was disabled in both cases.

The impact of radiation reaction on the particle trajectory can be drastic. Indeed, a comparison of two trajectories of particles injected at the same location within the magnetosphere where one has radiation reaction enabled, ‘rr’, and the other does not, ‘no rr’, is shown in Fig. 12. The upper row shows a crashed particle, the middle panel shows an ejected particle, and the bottom panel shows a trapped particle. The left column of the figure projects the trajectory on the (x, y) plane, and the middle column projects the trajectory on the (x, z) plane. The particle with radiation reaction enabled is shown with red and the particle without it is blue. The right column shows the time evolution of the Lorentz factor in both cases. The labels crashed, trapped, or ejected refer to the trajectory as observed in the ‘no rr’ case. Simulations were performed for electrons. For the trapped case, the particles start by following a similar path until the one losing energy because of radiation reaction brakes and only drifts, whereas the particle without radiation reaction decelerates and then accelerates again in another direction (see lower panel). Its Lorentz factor oscillates between 10⁸ and 10¹³ in a periodic fashion associated with a bouncing motion during the trapped stage of the ‘no rr’ case. The particle with ‘rr’ loses a large fraction of its initial energy quickly, decreasing the Lorentz factor from 10⁸ to 10⁶. It then experiences violent oscillations to a point it almost rests, and then it accelerates again. For the crashed case, both particles follow similar tracks, although they are slightly different in the (x, y) plane (top panel of Fig. 12). Nevertheless, the Lorentz factor is four to five orders of magnitude lower in the radiation reaction case. While moving towards the star, the electron efficiently radiates its energy gained in the increasingly stronger electromagnetic field. The largest difference in particle trajectories was observed between a particle ejected without radiation reaction that actually crashes when radiation is enabled (middle panel). The asymptotic Lorentz factor in the ejected case is about 10¹³, whereas the Lorentz factor in the crashed case behaves as in the trapped and crashed case shown in the bottom row of Fig. 12.

Fig. 12. Two trajectories of the same particle injected at the same location but with radiation reaction in one case(‘rr’, in red solid lines) and without it in the other (‘no rr’, in dash-dotted blue lines). The upper panel corresponds to a crashed particle, the middle panel to an ejected particle, and the bottom panel to a trapped particle. The left column projects the trajectory on the (x, y) plane, the middle column projects it on the (x, z) plane, and the right column shows the evolution of the Lorentz factor in a log scale.

The properties of the particles can also be viewed according to their final state: crashed, ejected, or trapped. This point of view is explored more deeply in the following sub-section.

Particles impacting the neutron star. We also investigated the final positions of particles when crashing onto the stellar surface. Thus, as shown in Fig. 4, we report the hotspots on the polar caps, which highlights the fact that the protons impact the star in very localised areas concentrated around the magnetic axis. Indeed, those areas are always located around the azimuth ϕ = 0° and ϕ = 180°, and they respectively stay in the northern and southern hemispheres as χ increases. The spot at ϕ = 0° is always found between θ = 0° and θ = 90°, while the spot at ϕ = 180° is found between θ = 90° and θ = 180°. By comparison, if the radiation reaction is disabled, as in Fig. 13, the hotspots move slightly along the meridian as χ changes. When χ = 30°, the hotspots are respectively at θ = 45° and θ = 135°, while for χ = 90°, they are respectively centred around θ = 108° and θ = 72°. We emphasise that there exists no south-north hemisphere symmetry in these impact maps, neither for protons nor for electrons. This is due to the nature of the electromagnetic field, which is a vector field, and it does not produce the same pattern when χ > 90° compared to χ < 90°.

Fig. 13. Same as Fig. 4 but with radiation reaction disabled.

Trapped particles. The structures formed by trapped particles remained close to the neutron star. These are shown in Figs. 6 and 7. Again, we noticed that these particles tend to avoid some regions while populating other well-defined areas. We also note that the order of the figure is for increasing χ for protons but decreasing χ for electrons. This highlights the symmetry between positive and negative charges when switching from χ to π − χ.

As shown in Fig. 14, particles have quite similar positions whether radiation reaction is enabled or not, with the same regions being populated but some particles being farther away from the neutron star when radiation reaction is disabled. We also noticed that despite having fewer particles for the simulations, the particles cover a wider area when radiation reaction is not enabled.

Fig. 14. Map of the final colatitude, azimuth, and radius (colour) of protons around neutron stars of inclination 60° (a), 90° (b), 120° (c), and 150° (d). Radiation reaction was disabled, and the magnetic axis is in the ϕ = 0 plane.

Actually, Fig. 6 shows that as χ increases, the trapped particles are on average repelled farther away from the surface of the neutron star. For χ = 60°, all these particles are extremely close to the surface of the neutron star (r ∈ [0.105; 0.135]), and it is probable that given more time, these particles would eventually crash onto the surface because of the energy losses. We noticed that the particles formed over densities around ϕ = 0, θ = 2π/3 and ϕ = π, θ = π/3. For χ = 90°, the structure has a spiralling tail of particles that seems to trail towards the surface of the neutron star, while most particles are in a more densely populated region at slightly higher altitudes. Again, these particles are concentrated around ϕ = 0, θ = 3π/4 and ϕ = π, θ = π/4.

Particles are reputed to be trapped for the time of the simulation. What happens later could not be guessed. However, this final state could depend on the duration of the trapping before being ejected or crashing onto the surface. In order to check the ‘stability’ of the trapping state, we performed new simulations with half the total time and twice the total time of the fiducial run, respectively t_f/P = 7.5 and t_f/P = 30. In this further analysis, we injected 2048 electrons and let them evolve in an orthogonal rotator (χ = 90°) for half and twice the time of all other simulations. It showed that as time grows, the particles spiral even more, as Fig. 15 highlights. However, despite their proximity to the neutron star surface, the fraction of particles impacting the surface does not vary. The statistics of the evolution of the particles are as follows: Out of the 2048 particles, 6 impact the surface, 1577 are trapped, and 465 are ejected. We did not notice any significant difference between these runs and concluded that the trapping state lasts for a significant time, at least several neutron star rotation periods, which is long enough to have an impact on the magnetosphere electrodynamics if there is any.

Fig. 15. Map of the final colatitude, azimuth, and radius (colour) of electrons around neutron stars with an inclination of 90°. The simulations lasted half (a) and twice (b) the time of other simulations. There were 2048 particles, radiation reaction was enabled, and the magnetic axis is in the ϕ = 0 plane.

For χ = 120°, the structures are still around ϕ = 0 and ϕ = π, but at θ = 4π/5 and θ = π/5 respectively, and around a dense area, one can notice a more diffuse region with fewer particles at low altitude. For χ = 150°, the structure takes the most different form and corresponds to the part of striped wind below the light cylinder radius. The above discussion for protons also applies to electrons if the values for χ are replaced by π − χ.

Ejected particles. Since the striped wind has been mentioned, we checked if the simulations manage to produce such a structure by plotting the final positions of particles outside the light cylinder radius, as done in Fig. 16.

Fig. 16. Map of the final positions of protons in the y, z plane for x ∈ [ − 1; 1] with χ = 120° (a) and χ = 150° (b), and in the x, y plane for z ∈ [ − 1; 1] with χ = 120° (c) and χ = 150° (d). Radiation reaction was enabled.

The distribution of the particles in space is very reminiscent to the striped wind geometry, especially due to the spiral. However, this wind geometry was retrieved only for χ = 120° and χ = 150° for protons (χ = 30° and χ = 60° for electrons). The theoretical equation describing the striped wind is given by Bogovalov (1999) and reads

$$\begin{array}{c}\hfill r_{\mathrm{s}}(t,\theta ,\varphi )=r_{\mathrm{L}}[\pm arccos(cot\theta cot\chi )+\frac{\mathit{ct}}{r_{\mathrm{L}}}-\varphi +2\ell \pi ],\end{array}$$(10)

where ℓ is an integer and r_s the radius of the distance of the particles composing the striped wind. We observed that the striped wind is found only in the angle θ ∈ [π/2 − χ; π/2 + χ] if χ ≤ π/2 and θ ∈ [π − χ; χ] if χ > π/2.

When comparing the proton positions to the expected theoretical wind structure, as in Fig. 17, we found that the spiral fits quite well. However, the opening angle at which the wind is spread is not the one expected. We believe that adding the particle interactions may improve the results by giving more realistic winds.

Fig. 17. Comparison of final proton positions to the expected striped winds in the y, z plane (a) and in the x, y plane (b). Radiation reaction was enabled.

Upon looking at the final colatitude and azimuth of protons in Fig. 8, we noticed that the most energetic particles prefer some directions but also that a given direction of ejection can be more densely populated or, conversely, almost void of particles. Comparing Figs. 8–18, we noticed some similarities regarding the more or less densely populated directions of ejection. However, the positions of the particles and their Lorentz factors are not the same any more than when radiation reaction is enabled.

Fig. 18. Map of the ejection colatitude and azimuth and Lorentz factor of protons on neutron stars of inclination 30° (a), 60° (b), 90° (c) and 120° (d), radiation reaction disabled.

Taking a look at Figs. 2, 4, 6, 8, 10, and 16, we found a symmetrical behaviour for particles injected on one side or the opposite side of the neutron star. The centre of the neutron star is a point of symmetry for the electric field and a point of anti-symmetry for the magnetic field, meaning that E(x, y, z, t) = − E(−x, −y, −z, t) and B(x, y, z, t) = B(−x, −y, −z, t). This property led us to find that particles injected symmetrically relative to the centre of the neutron star have symmetrical trajectories.

We considered two particles injected at x = (x, y, z) and x′=(−x, −y, −z) = − x with speeds of v(t = 0) = − v′(t = 0). We first find the force acting on these particles at the initial time step: the Lorentz force and the radiation reaction (here in classical formulation), which gives the following for the first particle:

$$\begin{array}{c}\hfill \mathit{F}=q(\mathit{E}(x,y,z)+\mathit{v}\times \mathit{B}(x,y,z))+\frac{{\mu }_0{q}^2}{6\pi c}\frac{{\mathrm{d}}^3\mathit{x}}{\mathrm{d}{t}^3},\end{array}$$(11)

and the following for the second particle:

$$\begin{array}{cc}\hfill \mathit{F}\prime & =q(\mathit{E}(-x,-y,-z)+\mathit{v}\prime \times \mathit{B}(-x,-y,-z))+\frac{{\mu }_0{q}^2}{6\pi c}\frac{{\mathrm{d}}^3\mathit{x}\prime }{\mathrm{d}{t}^3}\hfill \\ \hfill & =-q(\mathit{E}(x,y,z)+\mathit{v}\times \mathit{B}(x,y,z))-\frac{{\mu }_0{q}^2}{6\pi c}\frac{{\mathrm{d}}^3\mathit{x}}{\mathrm{d}{t}^3}.\hfill \end{array}$$(12)

Since the Lorentz force acting on the first particle F_L is the opposite of that acting on the second particle ${\mathit{F}\prime }_{\mathrm{L}}=-{\mathit{F}}_{\mathrm{L}}$, we found that $\frac{{\mathrm{d}}^3\mathit{x}}{\mathrm{d}{t}^3}=-\frac{{\mathrm{d}}^3\mathit{x}\prime }{\mathrm{d}{t}^3}$, meaning in the end that F = −F′. When we integrated the force, we found that the particles have symmetrical speeds relative to the centre of the neutron star: v = −v′, and if we integrated this speed, we found that x′= − x holds regardless of the time of integration.

4.3. Lorentz factor distribution function

In addition to the spatial particle distribution, in order to better understand the effects of radiation reaction on their dynamics, we compared our results to those previously obtained by Tomczak & Pétri (2020), knowing that simulations with the radiation reaction enabled yield more realistic results. One way to find a conservative upper limit to the Lorentz factor reached by particles around neutron star has been given in the introduction (Sect. 1). Taking the potential drop estimate ΔΦ = ΩBR² = 10¹⁶ V and multiplying it by e/mc², we found the Lorentz factor γ_e = 10¹⁰ for electrons and γ_p = 10^6.7 for protons.

This estimate is however not accurate because it assumes a constant static electric field and no radiation reaction. This limit is therefore rather conservative. Our simulations are very different because the electric field varies in time and radiation reaction is taken into account. When looking for a correlation between the potential drop and the final Lorentz factor reached by the electrons in an orthogonal rotator, we obtained the plot shown in Fig. 19, where the Lorentz factor is shown against the potential along the particle trajectory in log–log scale. The populations of crashed, trapped, and ejected particles are shown with green, blue, and red symbols, respectively. The number of crashed particles, less than ten, shown in green symbols, is too limited to perform any significant statistical analysis. For the two other populations, we found no evidence for a correlation between the final Lorentz factor of the particle and the potential drop, as Pearson’s correlation coefficient is r = 0.055 for trapped electrons and r = −0.017 for ejected electrons. This demonstrates that the Lorentz factor does not significantly depend on the particle motion history but is rather controlled by the local conditions (i.e., parallel accelerating electric field and curvature radius).

Fig. 19. Final Lorentz factor of electrons as a function of the potential drop (normalised to e V/me c2) along their path. 2048 particles were simulated, in green crosses for crashed particles, in blue crosses for trapped particles and in red plus symbols for ejected particles. The magnetic field is for an orthogonal rotator χ = 90°.

Without radiation reaction, the correlation is very strong. We found an excellent agreement between the potential drop and the final Lorentz factor, as shown in Fig. 20, for all three kinds of particles: trapped, ejected, and crashed.

Fig. 20. Same as Fig. 19 but without radiation reaction. The line y = x in orange depicts the Lorentz factor if the full potential drop is used. The correlation is clearly visible for all particles.

Another way to estimate the true Lorentz factor consists of equating the power of the electric force and that of the energy losses due to the curvature radiation. Thus, we get:

$$\begin{array}{c}\hfill q\mathbf{E}·\mathit{v}=\frac{q^2}{6\pi {ϵ}_0}{\gamma }^4\frac{c}{{\rho }^2},\end{array}$$(13)

where $\rho =\sqrt{R{r}_{\mathrm{L}}}$ is a typical radius of the curvature of the trajectory of the particle (here, that of a field line close to the stellar surface); v is the speed of the particle (∼c); and ϵ₀ is the vacuum permittivity. By solving for the Lorentz factor, we found that

$$\begin{array}{c}\hfill \gamma ={\left(\frac{6\pi {ϵ}_{0}E{\rho }^2}{q}\right)}^{1/4}={\left(\frac{q{E}_{\Vert }}{mc\mathrm{\Omega }}\frac{r_{\mathrm{L}}}{r_{\mathrm{e}}}{\stackrel{\sim }{\rho }}^2\right)}^{1/4}.\end{array}$$(14)

This last expression uses quantities without dimensions, such as $\stackrel{\sim }{\rho }=\rho /r_{\mathrm{L}}$, and the electric strength parameter, with r_e being the electron classical radius and E_∥ = β ⋅ E as the accelerating electric field. In our case, we applied the expression to millisecond pulsars and got γ = 10^7.5, which is only a guess because the curvature radius can be very different from that on the stellar surface. The true curvature κ is found from the velocity vector derivative such that

$$\begin{array}{c}\hfill \kappa =\frac{1}{\rho }=\parallel \frac{\mathrm{d}\mathit{v}}{c^2\mathrm{d}t}\parallel .\end{array}$$(15)

This expression accurately captures the local curvature radius ρ along the trajectory. Therefore, by following the Lorentz factor from the Landau-Lifshitz approximation and comparing it to the radiation reaction limit estimate as given by Eq. (14), we show that the latter always finds higher Lorentz factors, see Fig. 21. To check that the results converged, several different time step integration parameters were used, two times as well as five times smaller without noticeable changes. Thus, our results have converged and are robust.

Fig. 21. Lorentz factor from the Landau-Lifshitz approximation (LLR), in red and orange colours for different time steps, compared to the radiation reaction limit (RRL) guess as given by Eq. (14), in green and blue colours for different time steps.

We note that the Lorentz factors of trapped particles span a large range from almost rest γ ∼ 10 to γ ≲ 10⁸. This does not necessarily mean that they always experience strong radiation damping. Indeed, the Lorentz factor variation is two-fold in this electromagnetic field environment. First, the radiation reaction decelerates the particles from a very high Lorentz factor of γ ∼ 10^{12 − 13} to γ ∼ 10^7 − 8, decreasing by several orders of magnitude their initial Lorentz factor. Second, at a low to moderate Lorentz factor, the electric field component E_∥ parallel to the magnetic field can accelerate but can also decelerate particles depending on the sign of E_∥. Thus, moderate energies do not necessarily mean strong radiation reaction but rather efficient electric field deceleration. Fluctuation was observed in the Lorentz factor on short timescales for trapped particles, as seen in Fig. 21. Depending on the final time, we picked out a γ factor value between the minimum and maximum of the possible interval γ ∼ [10, 10^7 − 8]. This fluctuation is a kind of stroboscopic effect, giving a sample of the Lorentz factor spreading at this interval. We checked this effect by looking at a sample of ten trapped particles and found that the Lorentz factor drastically fluctuates on very short timescales in this energy range, see Fig. 22.

Fig. 22. Fluctuation in the Lorentz factor for a sample of ten particles zoomed around the time Ω t = 1.

As shown in Fig. 23, protons reach γ ≃ 10¹¹ and electrons reach γ ≃ 10¹⁴ without radiation reaction, while with radiation reaction, they respectively reach γ ≃ 10¹⁰ and γ ≃ 10^10.5, which is quite different from the curvature radiation approximation. We note, however, that the particles reaching the highest energies follow the field lines with the largest curvature radii, while for our calculation, we took an averaged radius of the curvature.

Fig. 23. Comparison of the Lorentz factor distributions for simulations carried out with radiation reaction disabled for protons (a) and electrons (b), and with radiation reaction enabled for protons (c) and electrons (d). The number of particles is 14 336 in (a) and (b) and 229 376 in (c) and (d).

When comparing the influence of the particle species, we noticed that particles with a high mass are less affected than those with a low mass for a given charge. Indeed, as Fig. 24 shows, the loss of energy is higher for electrons than for protons since the highest energy electrons lost ∼4.8 orders of magnitude for their Lorentz factor, whereas the highest energy protons lost only approximately one order of magnitude. Since the energy of a particle is E = γmc², we found that without radiation reaction, the proton energy is about E_p = 10^10.5 × m_p c² = 4, 75 J and the electron energy E_e = 10^13.8 × m_e c² = 5, 17 J. With radiation reaction enabled, we got E_p = 10^9.5 × m_p c² = 0.48 J and E_e = 10¹⁰ × m_e c² = 0.000082 J. This means that even if the particles had approximately the same energy without radiation reaction, the fastest electrons have only 0.017% of the energy of the fastest protons after radiation reaction has been enabled. We nonetheless note that with radiation reaction, the Lorentz factor distribution for protons near a pulsar of inclination χ is similar in shape to that of electrons near a pulsar of inclination π − χ. Two modes, one at low energy with low statistics and another one at higher energy with up to N ∼ 10 000 at the peak for the cases, are shown in Fig. 24. The main difference is the positions of the extrema. For protons near a pulsar with an inclination of 60°, the low energy peak is at γ = 10^7.5 and the high energy one is at γ = 10^8.9, while for electrons near a pulsar with an inclination of 120°, the low energy and high energy peaks are reached at γ = 10^7.8 and γ = 10^8.8, respectively.

Fig. 24. Comparison of the Lorentz factor distributions for electrons near a pulsar with an inclination χ = 120°, without (a) and with (b) radiation reaction enabled, and for protons near a pulsar with an inclination χ = 60°, without (c) and with (d) radiation reaction enabled.

Crashed particles. Figure 25 presents a closer look at particles impacting the neutron stars. Since the number of protons impacting the surface is lower than ten when χ ≥ 90°, we could only interpret the Lorentz factor distribution for χ = 30° and χ = 60°. For both inclinations, the protons reach at most γ = 10^9.4, the peak of the distributions is γ = 10⁸, and the local maximum is at γ = 10^8.7. The main difference between these spectral distributions is at low energy. For χ = 60°, the distribution starts at γ = 10^6.5, while for χ = 30°, the distribution starts at γ = 10⁷. Even if the statistics are too low, protons reach at most γ = 10^7.8 and at least γ = 10^6.8 for χ = 90°. However, for this case as well as the χ = 150° case, it is certain that if given more particles, the overall shape of the distributions would drastically change, according to the initial positions.

Fig. 25. Lorentz factor distribution of protons impacting neutron stars with an inclination χ = 30° (a), χ = 60° (b), χ = 90° (c), and χ = 120° (d). Radiation reaction was enabled.

When taking a look at Fig. 26, we found that the overall shape of the spectra is slightly different for protons and electrons. In particular, the electrons do not have the high energy local maximum. The values of the Lorentz factors reached by electrons are lower than those reached by protons. For χ = 30° and χ = 60°, the peak of the Lorentz factor distribution is reached at γ = 10^7.8, while for χ = 150° and χ = 120°, the distributions of electrons peak at γ = 10^7.7 and γ = 10^7.5, respectively.

Fig. 26. Lorentz factor distribution of protons impacting neutron stars with an inclination χ = 150° (a), χ = 120° (b), χ = 90° (c), and χ = 60° (d). Radiation reaction was enabled.

Trapped particles. The distribution of the Lorentz factors of protons trapped around neutron stars is shown in Fig. 27. It always produces a mode ending at γ = 10⁵, but for χ = 120 ° , a few protons reached γ = 10⁸, and for χ = 150°, more particles were distributed between γ = 10⁵ and γ = 10⁸. For χ = 60°, the total number of particles is quite low and analysing the distribution function becomes problematic. With χ = 90°, the maximum of the distribution is a plateau between γ = 10^3.5 and γ = 10^4.5, with N = 300 particles per bin. The distribution begins with a power law with a slope of approximately one between γ = 10 and γ = 10^3.5. With χ = 120°, the distribution also starts at γ = 10, but the maximum is located at γ = 10^3.9, with N ∼ 300, and a local maximum is reached at γ = 10^4.7, with N = 10 particles per bin. Finally, with χ = 150°, most bins contain between N = 1 and N = 20 particles, meaning that the distribution is highly sensitive to noise. The number of trapped particles becomes very low in this configuration. However, protons with an energy of about γ = 10⁸ reach up to N = 10 particles per bin, meaning that this part of the distribution is not simply a random event and that it may become significant with a higher number of simulated particles.

Fig. 27. Lorentz factor distribution of protons trapped around neutron stars with an inclination of χ = 60° (a), χ = 90° (b), χ = 120° (c), and χ = 150° (d). Radiation reaction was enabled.

Again, looking at Fig. 28, proton and electron spectra are different for mirror inclinations (i.e., χ and π − χ). The electron distributions are always bimodal, but due to low statistics, the cases with χ = 120° and χ = 30° are hard to interpret. Nevertheless, the protons hardly exceed γ = 10⁸ in any inclination. When χ = 90°, the two maxima are located at γ = 10³ and γ = 10^6.5, and both are close to N ∼ 300 particles per bin. The minimum between the two modes is located at γ = 10^5.5, with N ∼ 50 particles per bin. For χ = 60, the two maxima are also at the same statistical level of N ∼ 130 particles per bin, but at γ = 10^3.5 and γ = 10^6.5, the minimum between the modes is found at γ = 10⁵ for N ∼ 30 particles per bin.

Fig. 28. Lorentz factor distribution of electrons trapped around neutron stars of inclination χ = 120° (a), χ = 90° (b), χ = 60° (c), and χ = 30° (d). Radiation reaction was enabled.

Ejected particles. Figure 29 shows the spectral distribution of the Lorentz factor of protons ejected away from neutron stars. The case χ = 30° shows a power law with a slope of approximately two between γ = 10^7.9 (N ∼ 1) and γ = 10^8.9 (N ∼ 100) followed by an abrupt cut-off. The inclination χ = 60°, however, starts at γ = 10^8.4 (apart from a few particles below this Lorentz factor) and grows fast until γ = 10^8.6 at N ∼ 30, where the growth is slower until γ = 10^9.1 at N ∼ 70, which is where the distribution ends. Taking a look at the orthogonal rotator, the distribution starts with some particles at γ = 10⁶, but then at γ = 10⁷, the distribution follows a power law until γ = 10^8.7 (N = 200), giving a slope of ∼1.35. Then the distribution forms a plateau before a cut-off, ending at γ = 10^9.3. Regarding the χ = 120° inclination, we noticed a rapid growth from γ = 10⁸ up to γ = 10⁹ followed by a short plateau that starts decreasing at γ = 10^9.2 and ending at γ = 10^9.5. Finally, for χ = 150°, the distribution starts as a power law of slope ∼2.3 until γ = 10^9.3 and ends in γ = 10^9.4 with a sharp cut-off.

Fig. 29. Lorentz factor distribution of protons ejected by neutron stars with an inclination of χ = 30° (a), χ = 60° (b), χ = 90° (c), χ = 120° (d), and χ = 150° (e). Radiation reaction was enabled.

Figure 30 shows the Lorentz factor of ejected electrons. If χ = 150°, the distribution grows between γ = 10^7.5 and γ = 10^8.9 at N ∼ 60, and it then drops to zero at γ = 10^9.2. The distribution of the Lorentz factors for χ = 120° starts at γ = 10^8.2 with low statistics and some empty bins and grows to N ∼ 60 at γ = 10⁹, and then the distribution drops, first slowly until γ = 10^9.1 and then faster to zero in γ = 10^9.2. In the case of the orthogonal rotator, three particles were found at γ < 10⁸, but most of the distribution starts at γ = 10^8.2. It then peaks at γ = 10^9.1 with N ∼ 300 and ends in γ = 10^9.9. If χ = 60°, the distribution starts with a growth between γ = 10⁸ (N ∼ 1) and γ = 10^9.5 (N ∼ 500). It then decreases following a power law of slope ∼1.9 until N ∼ 15 at γ = 10^10.4. After that, the distribution stops, except for three particles close to γ = 10^10.6. Finally, if χ = 30°, the distribution grows irregularly between γ = 10^7.3 and the peak at γ = 10^9.5, with N ∼ 1500, and it then decreases irregularly until γ = 10^10.3.

Fig. 30. Lorentz factor distribution of electrons ejected by neutron stars with an inclination of χ = 150° (a), χ = 120° (b), χ = 90° (c), χ = 60° (d), and χ = 30° (e). Radiation reaction was enabled.

4.4. Aligned and anti-aligned cases

Aligned case. As the aligned case, corresponding to χ = 0°, is a static field, it was treated separately from the other cases. In addition, even if in reality there are open field lines due to the plasma surrounding the neutron star, the Deutsch field does not have such magnetic field lines. Regarding the results of the simulations, we also found that in this case, all the protons crash onto the surface of the neutron star.

Inspecting Fig. 31, we found that particles impact the neutron star on its poles. When analysing the initial positions of the particles and linking them to their final Lorentz factors, particles injected closer to the neutron star were found to reach a lower Lorentz factor than those injected farther away. This phenomenon is due to the potential drop being greater for particles farther away from the neutron star than for those close to the surface. For a given starting radius, protons injected closer to the equator reach a lower energy than those injected close to the poles. The most probable explanation for this is the energy loss due to the curvature radiation. For instance, for a given radius, as we get closer to the equator (θ = 10°), the curvature radius gets smaller, meaning that a particle injected there will have more energy loss than other particles starting at the same radius.

Fig. 31. Proton dynamics for an aligned rotator (χ = 0°). Panel a: impact map on the neutron star. Panel b: Lorentz factor distribution of the protons impacting the surface. Panel c: initial colatitude and azimuth, and final Lorentz factor of the protons starting at r ∈ [0.8; 0.9]. Panel d: same as panel c but for protons starting at r ∈ [0.3; 0.4].

When looking at the Lorentz factor distribution, we noticed that protons reach up to γ = 10^9.8 and never go lower than γ = 10^6.5. Also, the peak of the distribution is located at γ = 10^8.1, with N ∼ 1300, but a local maximum was also found at γ = 10^8.6, with N ∼ 130.

Anti-aligned case. The anti-aligned case was also treated separately for reasons similar to the aligned case (i.e., a constant field and no open magnetic field line in the case of the Deutsch field). However, the same comment could be made that the plasma normally around the neutron star changes the fields and ‘opens’ the magnetic field lines beyond the light cylinder.

Looking at the starting positions of protons in Fig. 32, we observed that particles close to the poles are accelerated more efficiently than those close to the equator. We believe that, just like for the aligned case, the curvature radiation is more important for particles close to the equator than for particles close to the poles.

Fig. 32. Proton dynamics for the anti-aligned rotator (χ = 180°). Panel a: final position and Lorentz factor of protons in the x ∈ [ − 1; 1] slice. Panel b: the same in the z ∈ [ − 1; 1] slice. Panel c: initial colatitude and azimuth and final Lorentz factor of protons starting at r ∈ [0.8; 0.9]. Panel d: Lorentz factor distribution of the protons around the neutron star.

Regarding the final positions of the protons, it appeared that a thick disc forms around the equator because of radiation losses. Indeed, protons tend to oscillate between the north and south poles (as in a Van-Allen radiation belt), but as the radius of curvature of the field lines gets smaller, the radiation losses become more intense, and a particle that was previously oscillating between the poles finally becomes stuck in the equatorial plane. Inversely, farther away from the neutron star, the curvature radius of the magnetic field lines is larger, so protons continue to oscillate between the north and south poles and do not lose enough energy to become stuck in the equatorial plane. Particles starting close to the poles follow the field lines with curvature radii so large that by the end of the simulation they did not have time to reach the equatorial plane (forming a dome in each hemisphere), and since the energy loss is low because of the high curvature radius, these particles reach the highest energy. Additionally, Fig. 32 shows not a disc but rather a disc and rings. The rings are in fact particles that are still oscillating between the north and south poles sufficiently so that they may reach coordinates out of the range z ∈ [ − 1; 1] at other times, and the holes between the rings are there since other particles have similar behaviours but were not in the range z ∈ [ − 1; 1] by the end of the simulation.

Finally, regarding the spectral distribution of the Lorentz factors, we highlight three modes. First is a low energy mode with a maximum at γ = 10^3.5, with N ∼ 500, and ranging from γ ∼ 10 to γ = 10^5.3, with N = 30. The second mode is at γ ∈ [10^5.3; 10^7.2] and has the overall maximum of the distribution at γ = 10^6.9, at N ∼ 800. The last mode is at γ ∈ [10^7.2; 10^8.4], and it has its maximum at γ = 10^7.9, with N ∼ 200.

For comparison, we took a look at the aligned case for electrons, too. Figure 33 shows that electrons also form a disc by the end of the simulation, but due to their lower masses, radiation reaction makes them lose more energy than the protons, meaning that the disc formed by electrons is thinner than that of the protons. The Lorentz factor distribution is composed of a more populated low energy mode from γ = 1 to γ = 10⁵ with the maximum at γ = 10^2.5, while the high energy mode is between γ = 10^6.5 and γ = 10⁸. It appears that the high energy electrons are those injected close to the poles (in terms of colatitude) because these particles follow magnetic field lines going farther away from the neutron star, meaning that these electrons spend less time in a strong magnetic field and thus lose less energy than those that remained close to the neutron star and ultimately formed the disc.

Fig. 33. Electron dynamics for the anti-aligned rotator (χ = 180°). Panel a: final position and Lorentz factor of electrons in the x ∈ [ − 1; 1] slice. Panel b: the same as (a) but in the z ∈ [ − 1; 1] slice. Panel c: initial colatitude and azimuth and final Lorentz factor of protons starting at r ∈ [0.8; 0.9]. Panel d: Lorentz factor distribution of the protons around the neutron star.

5. Conclusions

In this paper, we studied the influence of radiation reaction on proton and electron dynamics near millisecond pulsars. We showed the drastic impact of radiative losses onto their trajectories and Lorentz factor. First of all, the evolution of the particles and their positions at the end of the simulations may still share some similarities when radiation reaction is enabled or disabled, but many particles have another behaviour, and almost all of them end the simulations with positions that are different from those in the simulations without radiation reaction. For example, the hot spots have different shapes and are located at different colatitudes regardless of whether radiation reaction is enabled. Regarding trapped particles, their radial position seems to be more impacted by radiation reaction than their colatitudes or azimuths. Moreover, the positions of particles are even more different when they are ejected away from the neutron star. Regarding the Lorentz factor distributions, it clearly appears that particles lose energy due to radiation reaction, as the Lorentz factors reached by the particles are more realistic than when radiation reaction is neglected. The interaction between particles is the next and last step to study the evolution of particles around a neutron star with realistic fields. We hope that the retroaction of the particles on the fields will help in obtaining results even closer to reality for the particle speeds and positions.

Acknowledgments

We are grateful to the referee for helpful comments and suggestions. This work has been supported by CEFIPRA grant IFC/F5904-B/2018 and ANR-20-CE31-0010.

References

All Tables

All Figures

Fig. 1. Particle trajectory in a constant magnetic field with radiation reaction. Panel a: trajectory of the simulated particle (in blue) compared to theoretical positions (in orange) for τ0 ωB = 10−5, γ0 = 104, initial speed along y. Panel b: relative error as a function of the proper time step dτ for the simulations in blue points, compared to the first order error expectations in red solid line.

In the text

	Fig. 2. Final state of the protons depending on their initial positions around the pulsar. The obliquity and the initial position are χ = 60°, r0 ∈ [0.3; 0.4] for (a) and r0 ∈ [0.8; 0.9] for (b), χ = 90°, r0 ∈ [0.3; 0.4] for (c) and r0 ∈ [0.8; 0.9] for (d), radiation reaction being enabled.
In the text

	Fig. 3. Final state of the electrons depending on their initial positions around the pulsar. The obliquity and the initial position are χ = 120°, r0 ∈ [0.3; 0.4] for (a) and r0 ∈ [0.8; 0.9] for (b), χ = 90°, r0 ∈ [0.3; 0.4] for (c) and r0 ∈ [0.8; 0.9] for (d), radiation reaction being enabled.
In the text

	Fig. 4. Map of the impact spots and Lorentz factor of protons on the surface for inclination χ = 30° (a), 60° (b), 90° (c), and 120° (d). Radiation reaction was enabled, and the magnetic axis is in the ϕ = 0 plane.
In the text

	Fig. 5. Map of the impact spots and Lorentz factor of electrons on the surface for inclination χ = 150° (a), 120° (b), 90° (c), and 60° (d). Radiation reaction was enabled, and the magnetic axis is in the ϕ = 0 plane.
In the text

	Fig. 6. Map of the final colatitude, azimuth, and radius (colour) of protons around neutron stars of inclination 60° (a), 90° (b), 120° (c), and 150° (d). Radiation reaction was enabled, and the magnetic axis is in the ϕ = 0 plane.
In the text

	Fig. 7. Map of the final colatitude, azimuth, and radius (colour) of electrons around neutron stars of inclination 60° (a), 90° (b), 120° (c), and 150° (d). Radiation reaction was enabled, and the magnetic axis is in the ϕ = 0 plane.
In the text

	Fig. 8. Map of the ejection colatitude and azimuth and Lorentz factor of protons on neutron stars of inclination 30° (a), 60° (b), 90° (c), and 120° (d). Radiation reaction was enabled.
In the text

	Fig. 9. Map of the ejection colatitude and azimuth and Lorentz factor of electrons on neutron stars of inclination 150° (a), 120° (b), 90° (c), and 60° (d). Radiation reaction was enabled.
In the text

	Fig. 10. Final Lorentz factor of the particles shown in Fig. 2.
In the text

	Fig. 11. Final Lorentz factor (a) and state (b) of the particles near a neutron star of inclination χ = 60°, with r0 ∈ [0.8; 0.9], and final Lorentz factor (c) and state (d) of the particles near a neutron star of inclination χ = 90°, with r0 ∈ [0.8; 0.9]. Radiation reaction was disabled in both cases.
In the text

Fig. 12. Two trajectories of the same particle injected at the same location but with radiation reaction in one case(‘rr’, in red solid lines) and without it in the other (‘no rr’, in dash-dotted blue lines). The upper panel corresponds to a crashed particle, the middle panel to an ejected particle, and the bottom panel to a trapped particle. The left column projects the trajectory on the (x, y) plane, the middle column projects it on the (x, z) plane, and the right column shows the evolution of the Lorentz factor in a log scale.

In the text

	Fig. 13. Same as Fig. 4 but with radiation reaction disabled.
In the text

	Fig. 14. Map of the final colatitude, azimuth, and radius (colour) of protons around neutron stars of inclination 60° (a), 90° (b), 120° (c), and 150° (d). Radiation reaction was disabled, and the magnetic axis is in the ϕ = 0 plane.
In the text

	Fig. 15. Map of the final colatitude, azimuth, and radius (colour) of electrons around neutron stars with an inclination of 90°. The simulations lasted half (a) and twice (b) the time of other simulations. There were 2048 particles, radiation reaction was enabled, and the magnetic axis is in the ϕ = 0 plane.
In the text

	Fig. 16. Map of the final positions of protons in the y, z plane for x ∈ [ − 1; 1] with χ = 120° (a) and χ = 150° (b), and in the x, y plane for z ∈ [ − 1; 1] with χ = 120° (c) and χ = 150° (d). Radiation reaction was enabled.
In the text

	Fig. 17. Comparison of final proton positions to the expected striped winds in the y, z plane (a) and in the x, y plane (b). Radiation reaction was enabled.
In the text

	Fig. 18. Map of the ejection colatitude and azimuth and Lorentz factor of protons on neutron stars of inclination 30° (a), 60° (b), 90° (c) and 120° (d), radiation reaction disabled.
In the text

	Fig. 19. Final Lorentz factor of electrons as a function of the potential drop (normalised to e V/me c2) along their path. 2048 particles were simulated, in green crosses for crashed particles, in blue crosses for trapped particles and in red plus symbols for ejected particles. The magnetic field is for an orthogonal rotator χ = 90°.
In the text

	Fig. 20. Same as Fig. 19 but without radiation reaction. The line y = x in orange depicts the Lorentz factor if the full potential drop is used. The correlation is clearly visible for all particles.
In the text

	Fig. 21. Lorentz factor from the Landau-Lifshitz approximation (LLR), in red and orange colours for different time steps, compared to the radiation reaction limit (RRL) guess as given by Eq. (14), in green and blue colours for different time steps.
In the text

	Fig. 22. Fluctuation in the Lorentz factor for a sample of ten particles zoomed around the time Ω t = 1.
In the text

	Fig. 23. Comparison of the Lorentz factor distributions for simulations carried out with radiation reaction disabled for protons (a) and electrons (b), and with radiation reaction enabled for protons (c) and electrons (d). The number of particles is 14 336 in (a) and (b) and 229 376 in (c) and (d).
In the text

	Fig. 24. Comparison of the Lorentz factor distributions for electrons near a pulsar with an inclination χ = 120°, without (a) and with (b) radiation reaction enabled, and for protons near a pulsar with an inclination χ = 60°, without (c) and with (d) radiation reaction enabled.
In the text

	Fig. 25. Lorentz factor distribution of protons impacting neutron stars with an inclination χ = 30° (a), χ = 60° (b), χ = 90° (c), and χ = 120° (d). Radiation reaction was enabled.
In the text

	Fig. 26. Lorentz factor distribution of protons impacting neutron stars with an inclination χ = 150° (a), χ = 120° (b), χ = 90° (c), and χ = 60° (d). Radiation reaction was enabled.
In the text

	Fig. 27. Lorentz factor distribution of protons trapped around neutron stars with an inclination of χ = 60° (a), χ = 90° (b), χ = 120° (c), and χ = 150° (d). Radiation reaction was enabled.
In the text

	Fig. 28. Lorentz factor distribution of electrons trapped around neutron stars of inclination χ = 120° (a), χ = 90° (b), χ = 60° (c), and χ = 30° (d). Radiation reaction was enabled.
In the text

	Fig. 29. Lorentz factor distribution of protons ejected by neutron stars with an inclination of χ = 30° (a), χ = 60° (b), χ = 90° (c), χ = 120° (d), and χ = 150° (e). Radiation reaction was enabled.
In the text

	Fig. 30. Lorentz factor distribution of electrons ejected by neutron stars with an inclination of χ = 150° (a), χ = 120° (b), χ = 90° (c), χ = 60° (d), and χ = 30° (e). Radiation reaction was enabled.
In the text

	Fig. 31. Proton dynamics for an aligned rotator (χ = 0°). Panel a: impact map on the neutron star. Panel b: Lorentz factor distribution of the protons impacting the surface. Panel c: initial colatitude and azimuth, and final Lorentz factor of the protons starting at r ∈ [0.8; 0.9]. Panel d: same as panel c but for protons starting at r ∈ [0.3; 0.4].
In the text

	Fig. 32. Proton dynamics for the anti-aligned rotator (χ = 180°). Panel a: final position and Lorentz factor of protons in the x ∈ [ − 1; 1] slice. Panel b: the same in the z ∈ [ − 1; 1] slice. Panel c: initial colatitude and azimuth and final Lorentz factor of protons starting at r ∈ [0.8; 0.9]. Panel d: Lorentz factor distribution of the protons around the neutron star.
In the text

	Fig. 33. Electron dynamics for the anti-aligned rotator (χ = 180°). Panel a: final position and Lorentz factor of electrons in the x ∈ [ − 1; 1] slice. Panel b: the same as (a) but in the z ∈ [ − 1; 1] slice. Panel c: initial colatitude and azimuth and final Lorentz factor of protons starting at r ∈ [0.8; 0.9]. Panel d: Lorentz factor distribution of the protons around the neutron star.
In the text