# An new radiation reaction approximation for particle dynamics in the strong field regime

Jérôme Pétri

► To cite this version:

Jérôme Pétri. An new radiation reaction approximation for particle dynamics in the strong field regime. Astronomy and Astrophysics - A&A, 2023, 677, pp.A72. 10.1051/0004-6361/202346515. hal-04174884

HAL Id: hal-04174884

https://hal.science/hal-04174884

Submitted on 9 Nov 2023

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# A new radiation reaction approximation for particle dynamics in the strong field regime

J. Pétri

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

Received 28 March 2023 / Accepted 14 July 2023

## ABSTRACT

**Context.** Following particle trajectories in the intense electromagnetic field of a neutron star is prohibited by the large ratio between the cyclotron frequency $\omega_B$ and the stellar rotation frequency $\Omega$. No fully kinetic simulations on a macroscopic scale and with realistic field strengths have been performed so far due to the huge computational cost implied by this enormous scale of separation.

**Aims.** In this paper, we derive new expressions for the particle velocity subject to strong radiation reaction that are intended to be more accurate than the current state-of-the-art expression in the radiation reaction limit regime, the so-called Aristotelian regime.

**Methods.** We shortened the timescale hierarchy by solving the particle equation of motion in the radiation reaction regime, where the Lorentz force is always and immediately balanced by the radiative drag, and including a friction not necessarily opposite to the velocity vector, as derived in the Landau-Lifshitz approximation.

**Results.** Starting from the reduced Landau-Lifshitz equation (i.e., neglecting the field time derivatives), we found expressions for the velocity depending only on the local electromagnetic field configuration and on a new parameter related to the field strength that controls the strength of the radiative damping. As an example, we imposed a constant Lorentz factor $\gamma$ during the particle motion. We found that for ultra-relativistic velocities satisfying $\gamma \gtrsim 10$, the difference between strong radiation reaction and the radiation reaction limit becomes negligible.

**Conclusions.** The new velocity expressions produce results similar in accuracy to the radiation reaction limit approximation. We therefore do not expect this new method to improve the accuracy of neutron star magnetosphere simulations. The radiation reaction limit is a simple but accurate, robust, and efficient way to follow ultra-relativistic particles in a strong electromagnetic field.

**Key words.** magnetic fields – stars: neutron – relativistic processes – methods: numerical – acceleration of particles – radiation: dynamics

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

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

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

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

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

with $\tau$ being the proper time, $u^i = \gamma(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, $\varepsilon_0$ 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 $\gamma$, it is judicious to switch to the 3+1 formalism by introducing the observer time $dt = \gamma dr$. Therefore,

$$ \frac{d\mathbf{p}}{dt} = q\mathbf{F}_L + \gamma q \tau_m \left[ \frac{d\mathbf{E}}{dt} + v \wedge \frac{d\mathbf{B}}{dt} \right] + \frac{q^2 \tau_m}{m} \left[ \mathbf{F}_L \wedge \mathbf{B} + (\beta \cdot \mathbf{E}) \mathbf{E}/c \right] + \frac{q^2 \tau_m}{m c^2} \gamma^2 [(\beta \cdot \mathbf{E})^2 - \mathbf{F}_L^2] v \quad (3a) $$

$$ \frac{d\gamma}{dt} = \frac{q}{mc} \left[ \beta \cdot \mathbf{E} + \tau_m \gamma \beta \cdot \frac{d\mathbf{E}}{dt} + \frac{q \tau_m}{mc} (\mathbf{F}_L \cdot \mathbf{E} + \gamma^2 [(\beta \cdot \mathbf{E})^2 - \mathbf{F}_L^2]) \right], \quad (3b) $$

where we define the vector field

$$ \mathbf{F}_L = \mathbf{E} + v \wedge \mathbf{B}, \quad (4) $$

the normalised velocity $\beta = v/c$, and the momentum by $\mathbf{p} = \gamma m v$. In the constant field approximation, we drop the time derivatives and obtain the fundamental equation of motion for a particle as follows

$$ \frac{d\mathbf{p}}{dt} = q \mathbf{F}_L + \frac{q^2 \tau_m}{m} [\mathbf{F}_L \wedge \mathbf{B} + (\beta \cdot \mathbf{E}) \mathbf{E}/c] + \frac{q^2 \tau_m}{m c^2} \gamma^2 [(\beta \cdot \mathbf{E})^2 - \mathbf{F}_L^2] v \quad (5a) $$

$$ \frac{d\gamma}{dt} = \frac{q}{mc} \left[ \beta \cdot \mathbf{E} + \frac{q \tau_m}{mc} \left( \mathbf{F}_L \cdot \mathbf{E} + \gamma^2 [(\beta \cdot \mathbf{E})^2 - \mathbf{F}_L^2] \right) \right]. \quad (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

$$ \mathbf{F}^{\text{rad}} = K_2 [\mathbf{F}_L \wedge \mathbf{B} + (\beta \cdot \mathbf{E}) \mathbf{E}/c] - K_1 v \quad (6) $$

and balancing the Lorentz force

$$ \mathbf{F}^{\text{ext}} = q \mathbf{F}_L \quad (7) $$

with this radiation reaction $\mathbf{F}^{\text{ext}} = \mathbf{F}^{\text{rad}}$, we arrive at

$$ q \mathbf{F}_L = (K_1 + K_2 B^2) v - K_2 (\mathbf{E} \wedge \mathbf{B} + (\mathbf{B} \cdot v) \mathbf{B} + (\beta \cdot \mathbf{E}) \mathbf{E}/c). \quad (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_1$ and $K_2$ are deduced from Eq. (5a) and given by

$$ K_1 = \frac{q^2 \tau_m}{m c^2} \gamma^2 [\mathbf{F}_L^2 - (\beta \cdot \mathbf{E})^2] \quad (9a) $$

$$ K_2 = \frac{q^2 \tau_m}{m}. \quad (9b) $$

We note that these coefficients are algebraic, being positive whatever the sign of the charge $q$. However, $K_1$ depends on the Lorentz factor $\gamma$, and we leave it unconstrained. In order to solve Eq. (8), the velocity is advantageously decomposed into three components $(\sigma, \delta, \eta)$ such that

$$ v = \sigma \mathbf{E} + \delta \mathbf{B} + \eta \mathbf{E} \wedge \mathbf{B}. \quad (10) $$

These components must satisfy a linear system of three equations of unknowns $(\sigma, \delta, \eta)$ according to

$$ q(1 - \eta B^2) = [K_1 + K_2 B^2] \sigma - K_2 [\sigma E^2 + \delta (\mathbf{E} \cdot \mathbf{B})] / c^2 \quad (11a) $$

$$ q \eta (\mathbf{E} \cdot \mathbf{B}) = K_1 \delta - K_2 (\mathbf{E} \cdot \mathbf{B}) \sigma \quad (11b) $$

$$ q \sigma = [K_1 + K_2 B^2] \eta - K_2. \quad (11c) $$

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

$$ \delta = \frac{\mathbf{E} \cdot \mathbf{B}}{K_1} [q \eta + K_2 \sigma], \quad (12) $$
reducing the system to a 2 × 2 size. Indeed, the smaller linear system to be solved reads

$$ \left( K_1 + K_2 \left( B^2 - \frac{E^2}{c^2} \right) - \frac{K_2^2}{K_1} \left( \frac{E \mathbf{B}}{c} \right)^2 - q \right) \left[ B^2 - \frac{K_2}{K_1} \left( \frac{E \mathbf{B}}{c} \right)^2 \right] \left( \frac{\sigma}{\eta} \right) = \left( \frac{q}{K_2} \right)^2. \quad (13) $$

Deviation from the standard RRL arises because of the terms containing $K_2$, 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 $\omega$ 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/\omega$; the time scale $T = 1/\omega$; the magnetic field strength $B_n = m\omega/q$; the electric field strength $E_n = cB_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_2 = \omega\tau_m$ according to the following definitions

$$ a_B = \frac{B}{B_n} = \frac{\omega_B}{\omega} \quad (14a) $$

$$ a_E = \frac{E}{E_n} = \frac{\omega_E}{\omega}. \quad (14b) $$

The external force becomes

$$ \frac{\mathbf{F}_{\text{ext}}}{F_n} = \text{sign}(q) (\mathbf{e} + \boldsymbol{\beta} \wedge \mathbf{b}) \quad (15) $$

with $\text{sign}(q) = q/|q|$ and the radiative force

$$ \frac{\mathbf{F}_{\text{rad}}}{F_n} = k_2 [\mathbf{e} \wedge \mathbf{b} + \mathbf{b} \wedge (\mathbf{b} \wedge \boldsymbol{\beta}) + (\boldsymbol{\beta} \cdot \mathbf{e}) \mathbf{e}] - k_1 \boldsymbol{\beta} \quad (16) $$

with the normalised fields $\mathbf{e} = \mathbf{E}/E_n$, $\mathbf{b} = \mathbf{B}/B_n$ and

$$ k_1 = \frac{K_1}{|q| B_n}; \quad k_2 = \frac{K_2 B_n}{|q|} = \omega \tau_m. \quad (17) $$

The velocity expansion coefficients are also normalised according to

$$ \tilde{\sigma} = \sigma B_n; \quad \tilde{\eta} = \eta B_n^2; \quad \tilde{\delta} = \delta B_n/c. \quad (18) $$

The normalised system to be solved then reads with $\zeta = \text{sign}(q)$

$$ \left( k_1 + k_2 \left( b^2 - e^2 \right) - \frac{k_2}{k_1} (\mathbf{e} \cdot \mathbf{b})^2 - \zeta \left[ b^2 - \frac{k_2}{k_1} (\mathbf{e} \cdot \mathbf{b})^2 \right] \right) \left( \tilde{\sigma} \right) = \left( \zeta \frac{k_1 + k_2 b^2}{k_2} \right)^2. \quad (19) $$

In the above system, $k_2$ is fixed by the nature of the charged particle $(q, m)$ and the typical frequency $\omega$. There is no freedom to choose it arbitrarily. However $k_1$ is undetermined and needs to be fixed by an additional constraint on the velocity. Choosing the Lorentz factor $\gamma$, the coefficient $k_1$ is found from the condition $v^2/c^2 = 1 - \gamma^{-2}$.

Actually $k_1$ is related to the velocity $v$ by Eq. (9) in the LLR approximation. But solving for the coefficients $\sigma, \delta, \eta$, and $v$ is only a function of $k_1$. Thus Eq. (9) is a non-linear equation for $k_1$ 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_1$. 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 $d\mathbf{p}/dt$, we decided to cancel the velocity time derivative $d\mathbf{v}/dt$ given in the Landau-Lifshitz approximation by

$$ \gamma m \frac{d\mathbf{v}}{dt} = q[\mathbf{E} + \mathbf{v} \wedge \mathbf{B} - (\boldsymbol{\beta} \cdot \mathbf{E})\boldsymbol{\beta}] + \frac{q^2 \tau_m}{m} [\mathbf{E} \wedge \mathbf{B} + (\mathbf{v} \wedge \mathbf{B}) \wedge \mathbf{B} + (\boldsymbol{\beta} \wedge \mathbf{E}) \wedge \mathbf{E}/c + \boldsymbol{\beta} \cdot (\mathbf{E} \wedge \mathbf{B})\boldsymbol{\beta}] \quad (20) $$

$$ + \boldsymbol{\beta} \cdot (\mathbf{E} \wedge \mathbf{B})\boldsymbol{\beta}. \quad (21) $$

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

$$ \mathbf{E} + \mathbf{v} \wedge \mathbf{B} - (\boldsymbol{\beta} \cdot \mathbf{E})\boldsymbol{\beta} = 0. \quad (22) $$

Translated into normalised units, we get

$$ \zeta k_2 [\mathbf{e} \wedge \mathbf{b} + (\boldsymbol{\beta} \wedge \mathbf{b}) \wedge \mathbf{b} + (\boldsymbol{\beta} \wedge \mathbf{e}) \wedge \mathbf{e} + \boldsymbol{\beta} \cdot (\mathbf{e} \wedge \mathbf{b})\boldsymbol{\beta}] + \mathbf{e} + \boldsymbol{\beta} \wedge \mathbf{b} - (\boldsymbol{\beta} \cdot \mathbf{e})\boldsymbol{\beta} = 0. \quad (23) $$

This expression is quadratic in $\boldsymbol{\beta}$, 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 $||\mathbf{v}|| = c$ again, but any Lorentz factor can be imposed. Departure from the RRL arises due to the term involving $\zeta k_2$. 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 $\mathbf{v}$, we write

$$ q (\mathbf{E} + \mathbf{v} \wedge \mathbf{B}) = K \mathbf{v}, \quad (24) $$

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

$$ \left( B^2 + \frac{K^2}{q^2} \right) \mathbf{v} = \frac{K}{q} \mathbf{E} + \mathbf{E} \wedge \mathbf{B} + \frac{q}{K} (\mathbf{E} \cdot \mathbf{B}) \mathbf{B}. \quad (25) $$
Moreover *K* is the only positive solution of the bi-quadratic
equation

$$
K^4 v^2 - q^2 (E^2 - v^2 B^2) K^2 - q^4 (\mathbf{E} \cdot \mathbf{B})^2 = 0. \quad (26)
$$

Thus, it satisfies

$$
\frac{K}{|q|} = \sqrt{\frac{E^2 - v^2 B^2 + \sqrt{(E^2 - v^2 B^2)^2 + 4 v^2 (\mathbf{E} \cdot \mathbf{B})^2}}{2 v^2}} . \quad (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,

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

assuming particles moving at the speed of light $||v_{\pm}|| = c$. In the equation, $v_+$ represents positively charged particles, whereas $v_-$ represents negatively charged particles. The electromagnetic field strength $E_0$ and $B_0$ are deduced from the electromagnetic invariants $\mathbf{E}^2 - c^2 \mathbf{B}^2 = E_0^2 - c^2 B_0^2$ and $\mathbf{E} \cdot \mathbf{B} = E_0 B_0$ with the constraint $E_0 \ge 0$. Therefore, the radiated power is $\mathcal{P}_R = q \mathbf{F}_L \cdot \mathbf{v} = |q| c E_0$ and $K = |q| E_0/c \ge 0$.

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

$$
\mathcal{P}_R = \frac{q^2}{6\pi\epsilon_0} \gamma^4 \frac{c}{\rho_c^2} = \gamma^4 \frac{\tau_m m c^4}{\rho_c^2} = |q| c E_0 \quad (29)
$$

from which the Lorentz factor becomes

$$
\gamma = \left( \frac{|q| E_0}{\tau_m m c^3} \rho_c^2 \right)^{1/4} = \left( \frac{\boldsymbol{\beta} \cdot \mathbf{e}}{k_2} \frac{\rho_c^2}{r_L^2} \right)^{1/4} . \quad (30)
$$

Moreover, the curvature is found from the acceleration by

$$
\kappa_c = \frac{1}{\rho_c} \approx \left\| \frac{d\beta}{c dt} \right\|. \qquad (31)
$$

As long as $\gamma \gg 1$, the RRL Eq. (28) remains a very good approx-
imation.

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 radia-
tion reaction force term proportional to *γ*² is also opposite to the
velocity *v*. We could therefore identify *K*₁ with *K* to get

$$
K \tau_m v^2 = m c^2. \tag{32}
$$

Hence, $K/|q| \ge m/|q| \tau_m = 9 \times 10^{11}$ 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 \approx |q| E_0/c$ and

$$
\beta^2 \approx \frac{mc}{|q| E_0 \tau_m} = \frac{1}{\omega_{E_0} \tau_m}. \quad (33)
$$

Keeping the velocity less than the speed of light leads to $E_0 \ge 2.7 \times 10^{20} \text{ V m}^{-1}$, which is even larger than the critical value of $E_{\text{crit}} \approx 1.3 \times 10^{18} \text{ V m}^{-1}$. Therefore, this idea fails and gives the same value as before for the magnetic equivalent of $c E_0 = 9 \times 10^{11} \text{ 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 pre- scription 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_2 = 0$. The system (19) then
simplifies into

$$
\begin{pmatrix} k_1 & \zeta b^2 \\ -\zeta & k_1 \end{pmatrix} \begin{pmatrix} \tilde{\sigma} \\ \tilde{\eta} \end{pmatrix} = \begin{pmatrix} \zeta \\ 0 \end{pmatrix}. \qquad (34)
$$

The solution is readily found with

$$
\tilde{\eta} = \frac{1}{k_1^2 + b^2} \qquad (35a)
$$

$$
\tilde{\sigma} = \zeta k_1 \tilde{\eta} \qquad (35b)
$$

$$
\tilde{\delta} = \zeta \frac{\mathbf{e} \cdot \mathbf{b}}{k_1} \tilde{\eta}. \qquad (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_1$ needs to be deduced, for instance, from the Lorentz factor, as
in the previous discussion about curvature radiation power.

When $k_2 \neq 0$, $k_1$ 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_1$ in the system (19) is

$$
k_1 \approx \frac{E_0}{E_n} = \frac{|q| E_0}{m c \omega} = a_{E_0}. \qquad (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 inde-
pendent of $\gamma$.

3.3. Lorentz factor from LLR

The system of equations (19) solves the particle velocity vec-
tor by assuming the coefficient $k_1$ 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 $\gamma$.
We need to solve for $\gamma$ such that Eq. (9) is verified. However, as
we show in the next section, the IRR algorithm finds very simi-
lar velocities compared to the ‘standard’ algorithm. In this case
Eq. (24) also holds, approximately. Then it can be shown that

$$
\gamma^2 q^2 [(\boldsymbol{\beta} \cdot \mathbf{E})^2 - \mathbf{F}_L^2] \approx -K_1^2 v^2, \quad (37)
$$

which gives values of $K_1$ very different from the expectation in
Eq. (9).

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

$$
\beta \cdot \mathbf{E} + \frac{q \tau_m}{m c} \left( \mathbf{F}_L \cdot \mathbf{E} + \gamma^2 \left[ (\beta \cdot \mathbf{E})^2 - \mathbf{F}_L^2 \right] \right) = 0. \quad (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) depend-
ing only on $\gamma$. However, here also, as the solution is close to the
expression (28), we found instead that

$$
\mathbf{F}_L \cdot \mathbf{E} + \gamma^2 [(\boldsymbol{\beta} \cdot \mathbf{E})^2 - \mathbf{F}_L^2] \approx 0, \quad (39)
$$

which is not compatible with Eq. (38).

Actually, the second term in Eq. (38) corresponds to the
opposite of the curvature radiation power $\mathcal{P}_R$. It is related to the
curvature $\kappa_c$ in the case of an ultra-relativistic particle such that

$$
\kappa_c = \frac{|q|}{\gamma^2 m c^2} \sqrt{\gamma^2 [\mathbf{F}_L^2 - (\boldsymbol{\beta} \cdot \mathbf{E})^2] - \mathbf{F}_L \cdot \mathbf{E}}. \quad (40)
$$

If the term $\mathbf{F}_L \cdot \mathbf{E}$ is negligible, we retrieve the result (Kelner et al. 2015)

$$
\kappa_c \approx \frac{|q|}{\gamma m c^2} \sqrt{\mathbf{F}_L^2 - (\boldsymbol{\beta} \cdot \mathbf{E})^2}. \quad (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 $\mathcal{P}_R$, again
by replacing $E_0$ by $\boldsymbol{\beta} \cdot \mathbf{E}$ in Eqs. (29) and (30).

A final trial consisted of integrating the Lorentz factor dif-
ferential 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 acceler-
ate 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

$$
\frac{d\mathbf{E}}{dt} = \frac{\partial \mathbf{E}}{\partial t} + \mathbf{v} \cdot \nabla \mathbf{E} \quad (42)
$$

and with a similar expression for **B**. These expressions are, how-
ever, unwieldy to implement because they require the comput-
ing of partial time and space derivatives $\partial_t$ and $\partial_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 $\Omega$ and set $\omega = \Omega$. The numerical setup, elec-
tromagnetic field configuration, and initial conditions for parti-
cle 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_2 b^2 / k_1$ 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

<figure><img src="image_2.png" /><figcaption>Fig. 1. Ratio of the coefficient $k_1$ and $k_2$ expressed as $k_2 b^2/k_1$ in log scale for a sample of eight trajectories starting at several distances from the surface given by $r_0/r_L \approx \{1, 0.37, 0.14, 0.05\}$.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 2. Time evolution of the ratio $c K_1 / |q| E_0$ in log scale for a sample of eight representative trajectories. Both numbers are identical to eight digits of precision.</figcaption></figure>

the **E**, **B**, and **E** ∧ **B** directions. To check if this situation hap-
pens for electrons in the rotating magnetic dipole, we plotted
this ratio in a log scale (see Fig. 1) for a sample of eight tra-
jectories starting at different locations $r_0$ within the light cylin-
der. The radial distances are given in the legend of the figure
and were normalised to the light cylinder radius $r_L$ such that
$r_0/r_L \approx \{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.

Figure 2 shows the deviation of $K_1$ from $|q|E_0/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_\parallel = \beta \cdot \mathbf{E}$ to the value $E_0$ corresponding to the parallel electric field in the strict radiation reaction regime. Because $E_\parallel < 0$ for negatively charged particles, we plotted $\zeta E_\parallel/E_0$ to keep positive numbers. Both values of the parallel electric field are identical to more than 12 digits of precision.
<figure><img src="image_4.png" /><figcaption>Fig. 3. Time evolution of the true parallel electric field $E_\parallel$ compared to the estimated parallel electric field $E_0$ for a sample of eight trajectories. Because $E_\parallel < 0$ is negative for electrons, we plot $\zeta E_\parallel/E_0$.</figcaption></figure>

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

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 $\kappa_c$ in Eq. (31) was estimated with a finite difference approximation

$$ \kappa_c = \left\| \frac{v^{n+1/2} - v^{n-1/2}}{c^2 dt} \right\|. \qquad (43) $$

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.

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_0/r_L \approx \{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,

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

corresponding to the first row with $r_0/r_L \approx 1$ and to the second row with $r_0/r_L \approx 0.37$, all trajectories agree and overlap. Only some trajectories starting from the stellar surface, corresponding to the fourth row with $r_0/r_L \approx 0.05$, do not match, though the behaviour remains the same. For particles starting at $r_0/r_L \approx 0.14$, third row, the agreement is also excellent.

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 $d\ln E/dt$ and $d\ln B/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 $\gamma$. A simple criterion for dropping these terms is $\gamma k_2 \ll 1$. Thus, we can confidently ignore these time derivatives even outside the light cylinder.
<figure><img src="image_7.png" /><figcaption>**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 $r_0/r_L \approx \{1, 0.37, 0.14, 0.05\}$, from top to bottom.</figcaption></figure>

**5. Conclusions**

Tracking a charged particle motion in an ultra-strong elec-
tromagnetic 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
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Time evolution of the electric field derivative $d\ln E/dt$ (blue) and magnetic field derivative $d\ln B/dt$ (red), in normalised units, along a sample of three particle trajectories depicted by solid lines, dashed lines, and dotted lines.</figcaption></figure>

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 $\gamma^2$. It almost always remains negligible compared to the $\gamma^2$ 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_0 \gtrsim 10^{22} \text{ W cm}^{-2}$ (Gonoskov et al. 2022), corresponding to magnetic field strengths on the order of $B \gtrsim 10^7 \text{ 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.

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

Abdo, A. A., Ajello, M., Allafort, A., et al. 2013, ApJS, 208, 17

Cai, Y., Gralla, S. E., & Paschalidis, V. 2022, arXiv e-prints [arXiv:2209.07469]

Chang, S., Zhang, L., Jiang, Z., & Li, X. 2022, MNRAS, 513, 925

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

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

Gonoskov, A., Blackburn, T., Marklund, M., & Bulanov, S. 2022, Rev. Mod. Phys., 94, 045001

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

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

Mestel, L. 1999, Stellar Magnetism (Oxford: Clarendon)

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

Philippov, A., & Kramer, M. 2022, ARA&A, 60, 495

Pétri, J. 2019, MNRAS, 484, 5669

Pétri, J. 2022, A&A, 666, A5

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