© The Authors 2023

1. Introduction

With the recent increase in computational power, performing full kinetic simulations of neutron star magnetospheres can now be envisaged. However, the difference in timescales between gyro frequency and stellar frequency prevents realistic values from being applied to these parameters. So far the only way to circumvent this scaling problem is to downsize these frequencies while still keeping them well separated by respecting the ordering of these frequencies. Although this is helpful for understanding the dynamics of charged particles in extreme environments, it does not permit estimations of the true efficiency of particle acceleration and radiation reaction to be made because the Lorentz factors reached are several orders of magnitude lower than the ones predicted from observations of very high energy photons (see the review by Philippov & Kramer 2022). Recently some attempts have been made to simulate realistic parameters by Tomczak & Pétri (2020) and Pétri (2022), but the computational time remains prohibitive, and only test particles have been investigated while neglecting their back reaction to the field.

It is highly desirable to overcome this limitation by employing an approximation known as the radiation reaction limit (RRL) regime, sometimes also called Aristotelian dynamics, for which the equation of motion with radiative friction is shortened by use of an algebraic expression for the particle velocity depending only on the local value of the electric and magnetic field. This idea was applied by Mestel et al. (1985) and Finkbeiner et al. (1989). Spectra and light curves in this regime were extensively studied by Pétri (2019) in a vacuum field. He found realistic Lorentz factors and photon energies in reasonable agreement with the spectra observed by Fermi/LAT (Abdo et al. 2013).

Recently, Chang et al. (2022) generalised the RRL velocity by including the Landau-Lifshitz term proportional to the velocity (Landau & Lifchitz 1989) and by computing the associated radiation spectra. They found a complicated formula that unfortunately does not apply to any electromagnetic field configurations. Moreover they introduced some hypotheses that are not well justified to derive an expression for the velocity. Following a different approach, Cai et al. (2022) studied the validity of the RRL equilibrium by describing the particle motion in a Frenet frame with a finite Lorentz factor. They introduced the principal null directions, which are the eigenvectors of the electromagnetic field tensors. The spatial part is equal to the Aristotelian spatial velocity, or stated differently, it is equal to the RRL velocity. Although their analysis is based on the Landau-Lifshitz equations, including the time evolution of the Lorentz factor, at the end of their derivation they had to resort to the computation of the curvature radius in order to estimate this aspect of the Lorentz factor. In this work, we attempt to estimate the Lorentz factor by evolving it in time from the initial conditions, but as we show, the curvature radiation interpretation leads to more accurate estimates of the Lorentz factor. Cai et al. (2022) applied their idea to a rather artificial magnetic field configuration. Our aim is to apply such techniques to realistic fields, such as a rotating magnetic dipole.

In this paper, we derive formulas for the velocity in an arbitrary electromagnetic field configuration starting from the reduced Landau-Lifshitz equation (LLR; i.e., neglecting the field time derivatives). In Sect. 2, we derive the algorithm for the new velocity field according to the LLRs and that we call improved radiation reaction (IRR). Some explicit expressions of this velocity are given in Sect. 3. In Sect. 4, we then quantify the improvement brought by the inclusion of the radiative friction term proportional to the velocity compared to the standard RRL. Conclusions and perspectives are touched in Sect. 5.

2. Strong radiation reaction regime

Radiation reaction can be thought of as a friction drag opposing some resistance to the Lorentz force. It acts as a brake and is appropriately depicted by a force opposite to the velocity vector. However, in the Landau-Lifshitz approximation, the radiation reaction force is opposite to the velocity only in the limit of ultra-relativistic particles. For an arbitrary particle speed, there are additional components along the electric field E, the magnetic field B, and the electric drift motion E ∧ B. We aim to quantify the effect of these additional forces in the particle trajectory by first deriving a new expression for the velocity.

2.1. Equation of motion

As an approximation of the Lorentz-Abraham-Dirac equation, we employ the Landau-Lifshitz expression according to Landau & Lifchitz (1989) such that

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

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

with the typical timescale related to the particle classical radius crossing time

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

with τ being the proper time, uⁱ = γ(c, v) as the 4-velocity, q being the particle charge, m as the mass, c as the speed of light, E and B as the electric and magnetic field, ε₀ as the vacuum permittivity, v as the particle velocity, and F^ik as the electromagnetic tensor.

To derive the velocity vector v and the Lorentz factor γ, it is judicious to switch to the 3+1 formalism by introducing the observer time dt = γ dτ. Therefore,

$$\begin{array}{cc}\hfill \frac{\mathrm{d}\mathit{p}}{\mathrm{d}t}& =q{\mathbf{F}}_{\mathrm{L}}+\gamma q{\tau }_{\mathrm{m}}[\frac{\mathrm{d}\mathbf{E}}{\mathrm{d}t}+\mathit{v}\wedge \frac{\mathrm{d}\mathbf{B}}{\mathrm{d}t}]\hfill \\ \hfill & +\frac{q^2{\tau }_{\mathrm{m}}}{m}[{\mathbf{F}}_{\mathrm{L}}\wedge \mathbf{B}+(\mathit{\beta }·\mathbf{E})\mathbf{E}/c]+\frac{q^2{\tau }_{\mathrm{m}}}{m{c}^2}{\gamma }^2[{(\mathit{\beta }·\mathbf{E})}^2-{\mathbf{F}}_{\mathrm{L}}^2]\mathit{v}\hfill \end{array}$$(3a)

$$\begin{array}{cc}\hfill \frac{\mathrm{d}\gamma }{\mathrm{d}t}& =\frac{q}{\mathit{mc}}[\mathit{\beta }·\mathbf{E}+{\tau }_{\mathrm{m}}\gamma \mathit{\beta }·\frac{\mathrm{d}\mathbf{E}}{\mathrm{d}t}+\frac{q{\tau }_{\mathrm{m}}}{mc}({\mathbf{F}}_{\mathrm{L}}·\mathbf{E}+{\gamma }^2[{(\mathit{\beta }·\mathbf{E})}^2-{\mathbf{F}}_{\mathrm{L}}^2])],\hfill \end{array}$$(3b)

where we define the vector field

$$\begin{array}{c}\hfill {\mathbf{F}}_{\mathrm{L}}=\mathbf{E}+\mathit{v}\wedge \mathbf{B},\end{array}$$(4)

the normalised velocity β = v/c, and the momentum by p = γ m v. In the constant field approximation, we drop the time derivatives and obtain the fundamental equation of motion for a particle as follows

$$\begin{array}{cc}\hfill \frac{\mathrm{d}\mathit{p}}{\mathrm{d}t}& =q{\mathbf{F}}_{\mathrm{L}}+\frac{q^2{\tau }_{\mathrm{m}}}{m}[{\mathbf{F}}_{\mathrm{L}}\wedge \mathbf{B}+(\mathit{\beta }·\mathbf{E})\mathbf{E}/c]\hfill \\ \hfill & +\frac{q^2{\tau }_{\mathrm{m}}}{m{c}^2}{\gamma }^2[{(\mathit{\beta }·\mathbf{E})}^2-{\mathbf{F}}_{\mathrm{L}}^2]\mathit{v}\hfill \end{array}$$(5a)

$$\begin{array}{cc}\hfill \frac{\mathrm{d}\gamma }{\mathrm{d}t}& =\frac{q}{\mathit{mc}}[\mathit{\beta }·\mathbf{E}+\frac{q{\tau }_{\mathrm{m}}}{mc}({\mathbf{F}}_{\mathrm{L}}·\mathbf{E}+{\gamma }^2[{(\mathit{\beta }·\mathbf{E})}^2-{\mathbf{F}}_{\mathrm{L}}^2])].\hfill \end{array}$$(5b)

2.2. Derivation of the velocity: First approach

The derivation of the particle velocity follows the procedure outlined by Mestel (1999). Nevertheless, instead of using a friction of the form −K v with K > 0, we use the three-dimensional version of the radiation reaction force, neglecting the space-time dependence of the electromagnetic field such that the radiative force reduces to the second and third term in the right-hand side of Eq. (5a). Writing the radiation reaction force as

$$\begin{array}{c}\hfill {\mathbf{F}}^{\mathrm{rad}}=K_2[{\mathbf{F}}_{\mathrm{L}}\wedge \mathbf{B}+(\mathit{\beta }·\mathbf{E})\mathbf{E}/c]-K_1\mathit{v}\end{array}$$(6)

and balancing the Lorentz force

$$\begin{array}{c}\hfill {\mathbf{F}}^{\mathrm{ext}}=q{\mathbf{F}}_{\mathrm{L}}\end{array}$$(7)

with this radiation reaction F^ext = F^rad, we arrive at

$$\begin{array}{c}\hfill q{\mathbf{F}}_{\mathrm{L}}=(K_1+K_2{B}^2)\mathit{v}-K_2(\mathbf{E}\wedge \mathbf{B}+(\mathbf{B}·\mathit{v})\mathbf{B}+(\mathit{\beta }·\mathbf{E})\mathbf{E}/c).\end{array}$$(8)

We note that there is no assumption about particles moving at the speed of light, their Lorentz factor is arbitrary, and v < c. This represents a novelty compared to all other radiation reaction expressions, which are always enforcing v = c. The coefficients K₁ and K₂ are deduced from Eq. (5a) and given by

$$\begin{array}{cc}& K_1=\frac{q^2{\tau }_{\mathrm{m}}}{m{c}^2}{\gamma }^2[{\mathbf{F}}_{\mathrm{L}}^2-{(\mathit{\beta }·\mathbf{E})}^2]\hfill \end{array}$$(9a)

$$\begin{array}{cc}& K_2=\frac{q^2{\tau }_{\mathrm{m}}}{m}.\hfill \end{array}$$(9b)

We note that these coefficients are algebraic, being positive whatever the sign of the charge q. However, K₁ depends on the Lorentz factor γ, and we leave it unconstrained. In order to solve Eq. (8), the velocity is advantageously decomposed into three components (σ, δ, η) such that

$$\begin{array}{c}\hfill \mathit{v}=\sigma \mathbf{E}+\delta \mathbf{B}+\eta \mathbf{E}\wedge \mathbf{B}.\end{array}$$(10)

These components must satisfy a linear system of three equations of unknowns (σ, δ, η) according to

$$\begin{array}{cc}& q(1-\eta B^2)=[K_1+K_2{B}^2]\sigma -K_2[\sigma E^2+\delta (\mathbf{E}·\mathbf{B})]/c^2\hfill \end{array}$$(11a)

$$\begin{array}{cc}& q\eta (\mathbf{E}·\mathbf{B})=K_1\delta -K_2(\mathbf{E}·\mathbf{B})\sigma \hfill \end{array}$$(11b)

$$\begin{array}{cc}& q\sigma =[K_1+K_2{B}^2]\eta -K_2.\hfill \end{array}$$(11c)

We recall that K₁ is unconstrained. Therefore, in order to fully solve the system, an additional condition is required for K₁. To this end, we could enforce v = c, but as a generalisation, we impose an user-defined Lorentz factor γ = (1 − v²/c²)^−1/2. Equations (10) and (11), supplemented with the condition on the speed v, fully determine the velocity vector v. Solving for δ, we get

$$\begin{array}{c}\hfill \delta =\frac{\mathbf{E}·\mathbf{B}}{K_1}[q\eta +K_2\sigma ],\end{array}$$(12)

reducing the system to a 2 × 2 size. Indeed, the smaller linear system to be solved reads

$$\begin{array}{c}\hfill \left(\begin{array}{cc}{\mathrm{K}}_1+{\mathrm{K}}_2({\mathrm{B}}^2-\frac{{\mathrm{E}}^2}{{\mathrm{c}}^2})-\frac{{\mathrm{K}}_2^2}{{\mathrm{K}}_1}{\left(\frac{\mathrm{E}·\mathrm{B}}{\mathrm{c}}\right)}^2& q[{\mathrm{B}}^2-\frac{{\mathrm{K}}_2}{{\mathrm{K}}_1}{\left(\frac{\mathrm{E}·\mathrm{B}}{\mathrm{c}}\right)}^2]\\ -q& {\mathrm{K}}_1+{\mathrm{K}}_2{\mathrm{B}}^2\end{array}\right)\left(\begin{array}{c}\sigma \\ \eta \end{array}\right)=\left(\begin{array}{c}q\\ K_2\end{array}\right).\end{array}$$(13)

Deviation from the standard RRL arises because of the terms containing K₂, which are usually neglected in the ultra-relativistic regime. In order to deduce the number of relevant free parameters in the problem, it is preferable to employ quantities without dimensions, as explained in the next section.

2.3. Dimensionless system

As generally required for numerical simulations, we introduce several useful quantities without dimensions relevant for the computation of the velocity. Following our previous work in Pétri (2022), the 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 time and length scale as well as electromagnetic field strengths such that the length scale L = c/ω; the time scale T = 1/ω; the magnetic field strength B_n = m ω/|q|; the electric field strength E_n = c B_n; and the typical electromagnetic force strength F_n = |q| E_n. The two important parameters defining the family of solutions are the field strength parameters a_B and a_E and the radiation reaction efficiency k₂ = ω τ_m according to the following definitions

$$\begin{array}{cc}& a_B=\frac{B}{B_n}=\frac{{\omega }_{\mathrm{B}}}{\omega }\hfill \end{array}$$(14a)

$$\begin{array}{cc}& a_E=\frac{E}{E_n}=\frac{{\omega }_{\mathrm{E}}}{\omega }.\hfill \end{array}$$(14b)

The external force becomes

$$\begin{array}{c}\hfill \frac{{\mathbf{F}}^{\mathrm{ext}}}{F_n}=\text{sign}(q)(\mathit{e}+\mathit{\beta }\wedge \mathit{b})\end{array}$$(15)

with sign(q) = q/|q| and the radiative force

$$\begin{array}{c}\hfill \frac{{\mathbf{F}}^{\mathrm{rad}}}{F_n}=k_2[\mathit{e}\wedge \mathit{b}+\mathit{b}\wedge (\mathit{b}\wedge \mathit{\beta })+(\mathit{\beta }·\mathit{e})\mathit{e}]-k_1\mathit{\beta }\end{array}$$(16)

with the normalised fields e = E/E_n, b = B/B_n and

$$\begin{array}{c}\hfill k_1=\frac{K_1}{|q|B_n};k_2=\frac{K_2{B}_n}{|q|}=\omega {\tau }_{\mathrm{m}}.\end{array}$$(17)

The velocity expansion coefficients are also normalised according to

$$\begin{array}{c}\hfill \stackrel{\sim }{\sigma }=\sigma B_n;\stackrel{\sim }{\eta }=\eta B_n^2;\stackrel{\sim }{\delta }=\delta B_n/c.\end{array}$$(18)

The normalised system to be solved then reads with ζ = sign(q)

$$\begin{array}{c}\hfill \left(\begin{array}{cc}{k}_1+k_2(b^2-e^2)-\frac{k_2^2}{k_1}{(\mathit{e}·\mathit{b})}^2& \zeta [b^2-\frac{k_2}{k_1}{(\mathit{e}·\mathit{b})}^2]\\ -\zeta & k_1+k_2{b}^2\end{array}\right)\left(\begin{array}{c}\stackrel{\sim }{\sigma }\\ \stackrel{\sim }{\eta }\end{array}\right)=\left(\begin{array}{c}\zeta \\ k_2\end{array}\right).\end{array}$$(19)

In the above system, k₂ is fixed by the nature of the charged particle (q, m) and the typical frequency ω. There is no freedom to choose it arbitrarily. However k₁ is undetermined and needs to be fixed by an additional constraint on the velocity. Choosing the Lorentz factor γ, the coefficient k₁ is found from the condition v²/c² = 1 − γ⁻².

Actually k₁ is related to the velocity v by Eq. (9) in the LLR approximation. But solving for the coefficients σ, δ, η, and v is only a function of k₁. Thus Eq. (9) is a non-linear equation for k₁ solely, which connects back to the fact that Eq. (8) is a non-linear equation for v involving the Lorentz factor. This first approach has the drawback of implicitly including the Lorentz factor in the linear system via the parameter K₁. In the next sub-section, we develop a second approach that is quadratic in the velocity and does not contain the Lorentz factor.

2.4. Derivation of the velocity: Second approach

In a second approach, called velocity radiation reaction (VRR), instead of cancelling the relativistic momentum time derivative $\frac{\mathrm{d}\mathit{p}}{\mathrm{d}t}$, we decided to cancel the velocity time derivative $\frac{\mathrm{d}\mathit{v}}{\mathrm{d}t}$ given in the Landau-Lifshitz approximation by

$$\begin{array}{cc}\hfill \gamma m\frac{\mathrm{d}\mathit{v}}{\mathrm{d}t}& =q[\mathbf{E}+\mathit{v}\wedge \mathbf{B}-(\mathit{\beta }·\mathbf{E})\mathit{\beta }]\hfill \\ \hfill & +\frac{q^2{\tau }_{\mathrm{m}}}{m}[\mathbf{E}\wedge \mathbf{B}+(\mathit{v}\wedge \mathbf{B})\wedge \mathbf{B}+(\mathit{\beta }\wedge \mathbf{E})\wedge \mathbf{E}/c\hfill \end{array}$$(20)

$$\begin{array}{cc}\hfill & +\mathit{\beta }·(\mathbf{E}\wedge \mathbf{B})\mathit{\beta }].\hfill \end{array}$$(21)

The advantage of this approach is that it sticks closer to the Aristotelian regime. Indeed, if the term involving τ_m is removed, we retrieve the RRL and the associated Aristotelian velocity expression that exactly satisfies

$$\begin{array}{c}\hfill \mathbf{E}+\mathit{v}\wedge \mathbf{B}-(\mathit{\beta }·\mathbf{E})\mathit{\beta }=0.\end{array}$$(22)

Translated into normalised units, we get

$$\begin{array}{cc}& \zeta k_2[\mathit{e}\wedge \mathit{b}+(\mathit{\beta }\wedge \mathit{b})\wedge \mathit{b}+(\mathit{\beta }\wedge \mathit{e})\wedge \mathit{e}+\mathit{\beta }·(\mathit{e}\wedge \mathit{b})\mathit{\beta }]\hfill \\ \hfill & +\mathit{e}+\mathit{\beta }\wedge \mathit{b}-(\mathit{\beta }·\mathit{e})\mathit{\beta }=0.\hfill \end{array}$$(23)

This expression is quadratic in β, and unlike the previous approach, it does not involve the Lorentz factor. It can thus be solved by standard root finding techniques for a fixed Lorentz factor (or equivalently a fixed velocity norm). In a simple prescription, we set the velocity norm to ∥v∥=c again, but any Lorentz factor can be imposed. Departure from the RRL arises due to the term involving ζ k₂. In the next section, we discuss the different approximations to the particle Lorentz factor.

3. Approximations of the particle Lorentz factor

In this section, we explore several approximations to estimate the particle velocity without resorting to a full time integration of the equation of motion. We first remind the standard expression for the velocity, how it compares to our new expression and then discuss the Lorentz factor estimation.

3.1. Friction opposite to velocity

Starting from the radiation reaction description of Mestel (1999) where the radiative friction is opposite to the particle velocity vector v, we write

$$\begin{array}{c}\hfill q(\mathbf{E}+\mathit{v}\wedge \mathbf{B})=K\mathit{v},\end{array}$$(24)

where K is a positive parameter related to the power radiated by the particle. Solving for the velocity, we find

$$\begin{array}{c}\hfill (B^2+\frac{K^2}{q^2})\mathit{v}=\frac{K}{q}\mathbf{E}+\mathbf{E}\wedge \mathbf{B}+\frac{q}{K}(\mathbf{E}·\mathbf{B})\mathbf{B}.\end{array}$$(25)

Moreover K is the only positive solution of the bi-quadratic equation

$$\begin{array}{c}\hfill K^4{v}^2-q^2(E^2-v^2{B}^2)K^2-q^4{(\mathbf{E}·\mathbf{B})}^2=0.\end{array}$$(26)

Thus, it satisfies

$$\begin{array}{c}\hfill \frac{K}{|q|}=\sqrt{\frac{E^2-v^2{B}^2+\sqrt{{(E^2-v^2{B}^2)}^2+4v^2{(\mathbf{E}·\mathbf{B})}^2}}{2v^2}}.\end{array}$$(27)

So far there are no constraints on the particle speed v < c. In the limit of v = c, we retrieve the velocity expression used in the literature, namely,

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

assuming particles moving at the speed of light ∥v_±∥=c. In the equation, v₊ represents positively charged particles, whereas v₋ represents negatively charged particles. The electromagnetic field strength E₀ and B₀ are deduced from the electromagnetic invariants ${\mathbf{E}}^2-c^2{\mathbf{B}}^2=E_0^2-c^2{B}_0^2$ and E ⋅ B = E₀ B₀ with the constraint E₀ ≥ 0. Therefore, the radiated power is 𝒫_R = q F_L ⋅ v = |q| c E₀ and K = |q| E₀/c ≥ 0.

Except for this ultra-relativistic limit for which we assume v = c, another relation is required to set the particle Lorentz factor γ. To this end, we equate the radiated power according to the local curvature radius ρ_c of the particle trajectory as

$$\begin{array}{c}\hfill {\mathcal{P}}_{\mathrm{R}}=\frac{q^2}{6\pi {\epsilon }_0}{\gamma }^4\frac{c}{{\rho }_{\mathrm{c}}^2}={\gamma }^4\frac{{\tau }_{\mathrm{m}}m{c}^4}{{\rho }_{\mathrm{c}}^2}=|q|c{E}_0\end{array}$$(29)

from which the Lorentz factor becomes

$$\begin{array}{c}\hfill \gamma ={\left(\frac{|q|E_0}{{\tau }_{\mathrm{m}}m{c}^3}{\rho }_{\mathrm{c}}^2\right)}^{1/4}={\left(\frac{\mathit{\beta }·\mathit{e}}{k_2}\frac{{\rho }_{\mathrm{c}}^2}{r_{\mathrm{L}}^2}\right)}^{1/4}.\end{array}$$(30)

Moreover, the curvature is found from the acceleration by

$$\begin{array}{c}\hfill {\kappa }_{\mathrm{c}}=\frac{1}{{\rho }_{\mathrm{c}}}\approx \left||\frac{\mathrm{d}\mathit{\beta }}{c\mathrm{d}t}\right||.\end{array}$$(31)

As long as γ ≫ 1, the RRL Eq. (28) remains a very good approximation.

In Eq. (24), the strength of the damping K is undetermined but usually set by the particle velocity v or equivalently by its Lorentz factor γ. Looking at LLR, we observed that the radiation reaction force term proportional to γ² is also opposite to the velocity v. We could therefore identify K₁ with K to get

$$\begin{array}{c}\hfill K{\tau }_{\mathrm{m}}{v}^2=m{c}^2.\end{array}$$(32)

Hence, K/|q|≥m/|q| τ_m = 9 × 10¹¹ T for electrons and positrons. This value is much too high. As a consequence, there are two equations, Eqs. (27) and (32), for the two unknowns K and v.

In the ultra-relativistic limit K ≈ |q| E₀/c and

$$\begin{array}{c}\hfill {\beta }^2\approx \frac{mc}{|q|E_0{\tau }_{\mathrm{m}}}=\frac{1}{{\omega }_{{\mathrm{E}}_0}{\tau }_{\mathrm{m}}}.\end{array}$$(33)

Keeping the velocity less than the speed of light leads to E₀ ≥ 2.7 × 10²⁰ V m⁻¹, which is even larger than the critical value of E_crit ≈ 1.3 × 10¹⁸ V m⁻¹. Therefore, this idea fails and gives the same value as before for the magnetic equivalent of c E₀ = 9 × 10¹¹ T. We must conclude that the only reasonable way to compute the Lorentz factor is via the curvature radiation power Eq. (29). Before switching back to the LLR equation, we check how the new velocity approximation compares to the simple prescription presented in this paragraph.

3.2. Comparison to the radiation reaction limit

In the literature about approximated radiation reaction formulas, only the force opposite to the velocity is considered. Translated into our more general approach dealing with the full set of terms in the LLR equation, we enforce k₂ = 0. The system (19) then simplifies into

$$\begin{array}{c}\hfill \left(\begin{array}{cc}{k}_1& \zeta b^2\\ -\zeta & k_1\end{array}\right)\left(\begin{array}{c}\stackrel{\sim }{\sigma }\\ \stackrel{\sim }{\eta }\end{array}\right)=\left(\begin{array}{c}\zeta \\ 0\end{array}\right).\end{array}$$(34)

The solution is readily found with

$$\begin{array}{cc}& \stackrel{\sim }{\eta }=\frac{1}{k_1^2+b^2}\hfill \end{array}$$(35a)

$$\begin{array}{cc}& \stackrel{\sim }{\sigma }=\zeta k_1\stackrel{\sim }{\eta }\hfill \end{array}$$(35b)

$$\begin{array}{cc}& \stackrel{\sim }{\delta }=\zeta \frac{\mathit{e}·\mathit{b}}{k_1}\stackrel{\sim }{\eta }.\hfill \end{array}$$(35c)

This solution is exactly the same as the one presented in Eq. (25) except that all quantities are now normalised. The speed is not explicitly imposed to be equal to the speed of light. Therefore, k₁ needs to be deduced, for instance, from the Lorentz factor, as in the previous discussion about curvature radiation power.

When k₂ ≠ 0, k₁ is the root of a polynomial of high degree with no analytical expression. We compute this coefficient by applying a root finding algorithm via Newton-Raphson method. A good initial guess for k₁ in the system (19) is

$$\begin{array}{c}\hfill k_1\approx \frac{E_0}{E_n}=\frac{|q|E_0}{mc\omega }=a_{E_0}.\end{array}$$(36)

We checked that very few iterations are required to converge to a highly accurate solution with several digits of precision. Some simulations are shown in the next section. Finally, in the last approach, we use the full terms in the original LLR equation and solve for the velocity while taking into account the term independent of γ.

3.3. Lorentz factor from LLR

The system of equations (19) solves the particle velocity vector by assuming the coefficient k₁ is freely adjustable. Actually, from the LLR equation, it is not tuneable and must be determined self-consistently with the expression (9) in which there are no free parameters once the velocity v is fixed. This approach would lead to a first algorithm for finding the particle Lorentz factor γ. We need to solve for γ such that Eq. (9) is verified. However, as we show in the next section, the IRR algorithm finds very similar velocities compared to the ‘standard’ algorithm. In this case Eq. (24) also holds, approximately. Then it can be shown that

$$\begin{array}{c}\hfill {\gamma }^2{q}^2[{(\mathit{\beta }·\mathbf{E})}^2-{\mathbf{F}}_{\mathrm{L}}^2]\approx -K_1^2{v}^2,\end{array}$$(37)

which gives values of K₁ very different from the expectation in Eq. (9).

In a second alternative algorithm using the LLR equation, we can try to set dγ/dt = 0 and solve for the value of K₁ such that it satisfies

$$\begin{array}{c}\hfill \mathit{\beta }·\mathbf{E}+\frac{q{\tau }_{\mathrm{m}}}{mc}({\mathbf{F}}_{\mathrm{L}}·\mathbf{E}+{\gamma }^2[{(\mathit{\beta }·\mathbf{E})}^2-{\mathbf{F}}_{\mathrm{L}}^2])=0.\end{array}$$(38)

Here again there are no free parameters once the velocity v is fixed. This equation constrains the Lorentz factor because it depends on v, which is fully solved once K₁ is fixed. Therefore, the procedure consists of finding the root of Eq. (38) depending only on γ. However, here also, as the solution is close to the expression (28), we found instead that

$$\begin{array}{c}\hfill {\mathbf{F}}_{\mathrm{L}}·\mathbf{E}+{\gamma }^2[{(\mathit{\beta }·\mathbf{E})}^2-{\mathbf{F}}_{\mathrm{L}}^2]\approx 0,\end{array}$$(39)

which is not compatible with Eq. (38).

Actually, the second term in Eq. (38) corresponds to the opposite of the curvature radiation power 𝒫_R. It is related to the curvature κ_c in the case of an ultra-relativistic particle such that

$$\begin{array}{c}\hfill {\kappa }_{\mathrm{c}}=\frac{|q|}{{\gamma }^{2}m{c}^2}\sqrt{{\gamma }^2[{\mathbf{F}}_{\mathrm{L}}^2-{(\mathit{\beta }·\mathbf{E})}^2]-{\mathbf{F}}_{\mathrm{L}}·\mathbf{E}}.\end{array}$$(40)

If the term F_L ⋅ E is negligible, we retrieve the result (Kelner et al. 2015)

$$\begin{array}{c}\hfill {\kappa }_{\mathrm{c}}\approx \frac{|q|}{\gamma m{c}^2}\sqrt{{\mathbf{F}}_{\mathrm{L}}^{\mathbf{2}}-{(\mathit{\beta }·\mathbf{E})}^{\mathbf{2}}}.\end{array}$$(41)

The curvature would vanish in the limiting case investigated in this section because, by construction, dv/dt = 0. Consequently, as in the previous section, the best procedure to compute the Lorentz factor is through the curvature radiation power 𝒫_R, again by replacing E₀ by β ⋅ E in Eqs. (29) and (30).

A final trial consisted of integrating the Lorentz factor differential Eq. (5b) in time from the initial conditions. Because the regime is close to the RRL, the second term expressing the power radiated almost always vanishes. Contrary to the magnetic field, only the electric field produces work and is able to accelerate particles. The results are less good compared to the curvature radius approach. Actually, the curvature radius represents only an auxiliary variable to compute the Lorentz factor. It could be derived straightforwardly from the definition of Eq. (31) but at the expense of computing the Lagrangian time derivatives of the electric and magnetic fields as

$$\begin{array}{c}\hfill \frac{\mathrm{d}\mathbf{E}}{\mathrm{d}t}=\frac{\partial \mathbf{E}}{\partial t}+\mathit{v}·\mathbf{\nabla }\mathbf{E}\end{array}$$(42)

and with a similar expression for B. These expressions are, however, unwieldy to implement because they require the computing of partial time and space derivatives ∂_t and ∂_r. We prefer to compute the curvature from a finite difference approximation of Eq. (31), which is an equivalent description but much simpler to implement numerically. In the next section, we explore the efficiency and accuracy of the above mentioned methods for a rotating magnetic dipole with an electric quadrupole component.

4. Simulations around a rotating dipole

As a typical macroscopic frequency, we used the neutron star rotation frequency Ω and set ω = Ω. The numerical setup, electromagnetic field configuration, and initial conditions for particle position and velocity are exactly the same as in Pétri (2022). We simulated a sample of test particles evolving in the Deutsch (1955) electromagnetic field.

4.1. Radiation reaction limit accuracy

The improved version of the radiation reaction regime in the first approach differs significantly from the standard version only whenever the ratio k₂ b²/k₁ becomes comparable or greater than one. This means that the braking force no longer aligns with the particle velocity vector and that it also involves friction in the E, B, and E ∧ B directions. To check if this situation happens for electrons in the rotating magnetic dipole, we plotted this ratio in a log scale (see Fig. 1) for a sample of eight trajectories starting at different locations r₀ within the light cylinder. The radial distances are given in the legend of the figure and were normalised to the light cylinder radius r_L such that r₀/r_L ≈ {1, 0.37, 0.14, 0.05}. As can be seen in this plot, this ratio is mostly much lower than one, meaning that the improved version of radiation reaction does not significantly differ from the straightforward RRL, except for sparse events of very few trajectories. Moreover, as we later show, even in these cases, the trajectories are not drastically affected by the corrections brought through the IRR expression.

Fig. 1. Ratio of the coefficient k1 and k2 expressed as k2 b2/k1 in log scale for a sample of eight trajectories starting at several distances from the surface given by r0/rL ≈ {1, 0.37, 0.14, 0.05}.

Figure 2 shows the deviation of K₁ from |q| E₀/c in the IRR version for the same sample shown in Fig. 1. The ratio equals one to very high accuracy. Both parameters are identical up to eight digits of precision. This supports the fact that the improvement is marginal. Finally, in Fig. 3, we compare the parallel electric field E_∥ = β ⋅ E to the value E₀ corresponding to the parallel electric field in the strict radiation reaction regime. Because E_∥ < 0 for negatively charged particles, we plotted ζ E_∥/E₀ to keep positive numbers. Both values of the parallel electric field are identical to more that 12 digits of precision.

Fig. 2. Time evolution of the ratio c K1/|q| E0 in log scale for a sample of eight representative trajectories. Both numbers are identical to eight digits of precision.

Fig. 3. Time evolution of the true parallel electric field E∥ compared to the estimated parallel electric field E0 for a sample of eight trajectories. Because E∥ < 0 is negative for electrons, we plot ζ E∥/E0.

Based on all the above observations, we did not expect to observe a drastic change in the particle dynamics between the RRL and IRR regime. For a more quantitative analysis, we plotted the Lorentz factor evolution in time for the same sample of electrons. No difference in the Lorentz factor evaluation was observed between both approximations. We therefore concluded that there is no advantage in including the correction brought by the Landau-Lifshitz equation with a friction term not anti-aligned with the velocity vector.

4.2. Comparison between IRR, VRR, and LLR

We have checked that the IRR does not bring significant improvements compared to radiation reaction. To complete the analysis of accuracy and efficiency of the IRR approximation, we compared it to the more reliable LLR equation of motion.

Figure 4 shows the evolution of the Lorentz factor for the three descriptions of the particle motion, RRL, IRR, and LLR. The curves only differ by their initial condition. The RRL and its improved version show Lorentz factor estimates agreeing with the LLR computations to reasonably good accuracy. We note that the time evolution of the Lorentz factor is reproduced with the associated fluctuations for one of the trajectories. For the RRL and IRR case, the Lorentz factors were computed according to expression (30). The curvature κ_c in Eq. (31) was estimated with a finite difference approximation

$$\begin{array}{c}\hfill {\kappa }_{\mathrm{c}}=\left||\frac{{\mathit{v}}^{n+1/2}-{\mathit{v}}^{n-1/2}}{c^2\mathrm{d}t}\right||.\end{array}$$(43)

Fig. 4. Evolution of the Lorentz factor for the three approximations of the equation of motion of an electron. Colours are as follows: RRL in dashed blue, IRR in dashed green, and LLR in solid red lines.

It only involved the value of the electromagnetic field at two neighbouring times: t^n + 1/2 and t^{n − 1/2}.

If we integrate the time evolution of the Lorentz factor instead, we get less accurate estimates of the Lorentz factor, as shown in Fig. 5 in green dashed lines for the VRR approach and indicated as VRR2. The blue dashed lines show the Lorentz factor computed via the curvature and give similar results to IRR in Fig. 4, indicated as VRR1.

Fig. 5. Evolution of the Lorentz factor for the VRR approximation in dashed green and dashed blue lines compared to LLR in solid red lines. VRR2 stands for evaluation by integration in time of the Lorentz factor, whereas VRR1 stands for evaluation of the Lorentz factor by the curvature radius.

Representing another check of the efficiency of the IRR and radiation reaction approximation, Fig. 6 overlaps the LLR trajectories shown in black solid lines onto the IRR trajectories shown in coloured, thick solid lines for a sample of particles starting at r₀/r_L ≈ {1, 0.37, 0.14, 0.05} from the top to the bottom row, respectively see the legend in the right-column panels). For the trajectories starting well above the stellar surface, corresponding to the first row with r₀/r_L ≈ 1 and to the second row with r₀/r_L ≈ 0.37, all trajectories agree and overlap. Only some trajectories starting from the stellar surface, corresponding to the fourth row with r₀/r_L ≈ 0.05, do not match, though the behaviour remains the same. For particles starting at r₀/r_L ≈ 0.14, third row, the agreement is also excellent.

Fig. 6. Comparison of IRR and LLR trajectories of a relevant sample of electrons. The solid black lines depict the LLR trajectories, and the coloured symbols show the IRR approximation. In each horizontal panel, particles start at a fixed spherical radius given by r0/rL ≈ {1, 0.37, 0.14, 0.05}, from top to bottom.

Finally, we stress that all the above results rely on the assumption that dE/dt = dB/dt = 0 are valid. This is correct as long as the terms involving dE/dt and dB/dt remain small compared to the other terms in Eq. (3). We checked this a posteriori by computing the time derivatives dlnE/dt and dlnB/dt, in normalised units, during a full simulation span. The time evolution of these derivatives is shown in Fig. 7. They remain of order unity, with a maximal value of about ten. If multiplied by the correct factors given in the Landau-Lifshitz equation, we noticed that these terms in the radiation reaction force indeed stay at a negligible level, even when multiplied by a factor γ. A simple criterion for dropping these terms is γ k₂ ≪ 1. Thus, we can confidently ignore these time derivatives even outside the light cylinder.

Fig. 7. Time evolution of the electric field derivative dlnE/dt (blue) and magnetic field derivative dlnB/dt (red), in normalised units, along a sample of three particle trajectories depicted by solid lines, dashed lines, and dotted lines.

5. Conclusions

Tracking a charged particle motion in an ultra-strong electromagnetic field is computationally a very demanding task. However, finding accurate approximations able to follow these ultra-relativistic trajectories with radiative friction is a central problem in modelling realistic neutron star magnetospheres. In this paper, we extended the velocity vector expression in the RRL by including a radiative force linear in velocity as derived from the Landau-Lifshitz equation. We showed that integrating the particle trajectories with this new expression gives very similar results to the ‘standard’ radiation reaction expression of the Aristotelian dynamics. The Lorentz factors are identical in both cases. A new parameter was introduced to control the strength of this force linear in velocity compared to the ultra-relativistic term proportional to γ². It almost always remains negligible compared to the γ² term anti-aligned with the velocity vector. Including such a refinement in the radiation reaction regime to obtain more accurate solutions is therefore not recommended because it also requires more computational time for no benefit.

Nevertheless, we observed some discrepancy between the Landau-Lifshitz solution and the IRR solution for some particles starting from regions close to the surface where the field strength is maximal. In such cases, the Landau-Lifshitz integration scheme is recommended if accuracy becomes an issue in obtaining reliable results. An alternative approach therefore would be to evolve the velocity vector in time and the Lorentz factor using the ultra-relativistic equation of motion approximation for a charged particle while assuming that the speed is and remains very close to the speed of light.

Another possible application beyond neutron stars but not explored in this work is using lasers in the extreme light regime to investigate high energy physics in ultra-strong electromagnetic fields in the laboratory. Indeed current technology pushes the laser nominal intensity above I₀ ≳ 10²² W cm⁻² (Gonoskov et al. 2022), corresponding to magnetic field strengths on the order of B ≳ 10⁷ T, which are similar to field strengths met around compact objects in high energy astrophysics. At such laser intensities, the field strength is expected to reach the radiation dominated regime and even the strong field quantum electrodynamics domain where electron-positron pair cascades are triggered.

Acknowledgments

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

References

All Figures

	Fig. 1. Ratio of the coefficient k1 and k2 expressed as k2 b2/k1 in log scale for a sample of eight trajectories starting at several distances from the surface given by r0/rL ≈ {1, 0.37, 0.14, 0.05}.
In the text

	Fig. 2. Time evolution of the ratio c K1/|q| E0 in log scale for a sample of eight representative trajectories. Both numbers are identical to eight digits of precision.
In the text

	Fig. 3. Time evolution of the true parallel electric field E∥ compared to the estimated parallel electric field E0 for a sample of eight trajectories. Because E∥ < 0 is negative for electrons, we plot ζ E∥/E0.
In the text

	Fig. 4. Evolution of the Lorentz factor for the three approximations of the equation of motion of an electron. Colours are as follows: RRL in dashed blue, IRR in dashed green, and LLR in solid red lines.
In the text

	Fig. 5. Evolution of the Lorentz factor for the VRR approximation in dashed green and dashed blue lines compared to LLR in solid red lines. VRR2 stands for evaluation by integration in time of the Lorentz factor, whereas VRR1 stands for evaluation of the Lorentz factor by the curvature radius.
In the text

	Fig. 6. Comparison of IRR and LLR trajectories of a relevant sample of electrons. The solid black lines depict the LLR trajectories, and the coloured symbols show the IRR approximation. In each horizontal panel, particles start at a fixed spherical radius given by r0/rL ≈ {1, 0.37, 0.14, 0.05}, from top to bottom.
In the text

	Fig. 7. Time evolution of the electric field derivative dlnE/dt (blue) and magnetic field derivative dlnB/dt (red), in normalised units, along a sample of three particle trajectories depicted by solid lines, dashed lines, and dotted lines.
In the text