The dynamics of RE (runaway electrons) in fusion plasmas span a wide range of temporal scales, from the fast gyro-motion, $ \u223c 10 \u2212 11 $ s, to the observational time scales, $ \u223c 10 \u2212 2 \u2192 1 $ s. To cope with this scale separation, RE are usually studied within the bounce-average or the guiding center approximations. Although these approximations have yielded valuable insights, a study with predictive capabilities of RE in fusion plasmas calls for the incorporation of full orbit effects in configuration space in the presence of three-dimensional magnetic fields. We present numerical results on this problem using the Kinetic Orbit Runaway electrons Code that follows relativistic electrons in general electric and magnetic fields under the full Lorentz force, collisions, and radiation losses. At relativistic energies, the main energy loss is due to radiation damping, which we incorporate using the Landau-Lifshitz formulation of the Abraham-Lorentz-Dirac force. The main focus is on full orbit effects on synchrotron radiation. It is shown that even in the absence of magnetic field stochasticty, neglecting orbit dynamics can introduce significant errors in the computation of the total radiated power and the synchrotron spectra. The statistics of collisionless (i.e., full orbit induced) pitch angle dispersion, and its key role played on synchrotron radiation, are studied in detail. Numerical results are also presented on the pitch angle dependence of the spatial confinement of RE and on full orbit effects on the competition of electric field acceleration and radiation damping. Finally, full orbit calculations are used to explore the limitations of gyro-averaging in the relativistic regime. To explore the practical impact of the results, DIII-D and ITER-like parameters are used in the simulations.

## I. INTRODUCTION

The potential impact that runaway electrons (RE) can have to the safe operation of ITER is currently well-recognized. In particular, RE can severely damage the plasma facing components during a major disruption, see, for example, Refs. 1–4 and references therein. As a result, in recent years, the study of RE dynamics in toroidal plasmas has been an area of significant interest.

Numerical studies of RE can be roughly classified based on the numerical method used, the number of degrees of freedom considered, the type of radiation losses and collisional transport operators included, and the level of self-consistent coupling. The preferred numerical methods are either particle-based or continuum. Particle-based methods track the trajectories of large ensembles of RE and incorporate collisional transport using Monte-Carlo sampling techniques. Continuum methods, on the other hand, consider grid-based or spectral solutions of partial differential equations (e.g., the Fokker-Planck equation) for the RE distribution function. In this paper, we follow a particle-based approach that has the well-known advantages of flexible numerical implementation and parallelization.

Among the most physically and computationally impactful element in numerical studies of RE is the number of degrees of freedom incorporated. Physics considerations including symmetries, conservation laws, and gaps in the spatio-temporal scales (from the gyro-period $ t \u223c 10 \u2212 11 $ s to the observational time scales $ t \u223c 10 \u2212 3 \u2192 1 $ s) as well as computational limitations are the main motivations for considering reduced models that track only a fraction of the full 6 degrees of freedom of each electron. The ultimate level of approximation is provided by “0-D” test particle models that follow the time evolution of the statistical mean (first moment) of the parallel and perpendicular momentum of the RE particle distribution integrated into space, neglecting momentum space diffusion.^{5,6} The next level of description is based on bounce-average approximations in the large aspect-ratio limit that eliminate all the spatial degrees of freedom and follow the individual dynamics of electrons in a 2-D momentum space.^{7} Early Monte-Carlo implementations of this approach include Refs. 8 and 9, and more recent continuum numerical implementations and analytical work include Refs. 10–12. Models based on 3-D descriptions aim to augment the 2-D phase space approach by incorporating 1-D radially dependent information at some level, see, for example, Ref. 13. The next level of description is provided by models based on the guiding center approximation that eliminates the gyro-motion and follows the spatial degrees of freedom of the guiding center along with the parallel momentum, e.g., Refs. 14 and 15. The most detailed description, which is the one adopted in the present paper, follows RE orbits in fully 3-D electric and magnetic fields using the exact relativistic Lorentz force, that is, the individual electron dynamics are evolved in a 6-D phase space. Examples of previous work following this approach include Ref. 16 and more recently Refs. 17 and 18. Our goal in the present paper is to study the role of spatial (3-D configuration space) effects in the dynamics of relativistic RE electrons using full orbit information, avoiding potential limitations of reduced degrees-of-freedom descriptions.

As mentioned above, the type of radiation losses and collisional transport operators included as well as the level of self-consistent coupling are also key elements of the description of RE dynamics. Here, we neglect the role of self-consistent effects, i.e., we assume the electric and magnetic field as given. Although this approximation might be too restrictive to investigate some aspects of RE, we adopt it here because our main goal is to explore what is missing in reduced descriptions that neglect the spatial dynamics of RE, and these descriptions also typically neglect self-consistent effects. On the other hand, synchrotron radiation emission and the resulting damping due to the radiation reaction force is a critical element in the understanding of the dynamics of RE especially in the highly relativistic regime which is the regime of interest in the present paper. The modeling of the radiation reaction force is a subtle problem. In particular, the Lorentz-Abraham-Dirac (LAD) force which is usually advertised as the most accurate first-principles model, exhibits serious physical inconsistencies, see, for example, Refs. 19 and 20. This is why the incorporation of radiation damping in RE studies uses different levels of approximations of the LAD. A byproduct of these approximations is the elimination of physically troubling terms, for example, the dependence on the third derivative of the particle position. To address this problem, here we use the Landau-Lifshitz formulation^{21} of the radiation-reaction that has been shown to be the most accurate, physically consistent approximation of the LAD force.^{19,20} However, the main contribution of the present paper is not primarily the use of more sophisticated radiation force models but the incorporation of spatially dependent effects. For example, in radiation damping models commonly used in 2-D phase space Fokker-Planck descriptions, the magnetic field information enters as a single parameter (typically approximated as the value at the magnetic axis) devoid of all the spatial geometric information. The main question we address is what is the impact of the full orbit effects, in general, and the spatial magnetic field information, in particular, on synchrotron radiation of RE in the highly relativistic regime. Although it is straightforward to incorporate collisions (using Monte-Carlo methods) within the particle-based approach, we have opted not to include them because the time scales of interest are several orders of magnitude shorter than the relativistic electron collision time scale $ \tau c o l l \u223c 10 \u2212 2 \u2009 s $. That is, the space-dependent RE physics of interest is essentially collisionless.

The numerical results presented here were obtained using the new code KORC (Kinetic Orbit Runaway electron Code). KORC is a modular Fortran 95, parallel code that we have recently developed to solve the dynamics of RE including radiation reaction forces and collisions (although, as mentioned above, the present results do not include collisions). The electric and magnetic fields accepted by KORC can be either defined analytically (like in the cases considered here) or pre-computed numerically (e.g., using the Variational Moments Equilibrium Code (VMEC), JFIT code, or any other equilibrium or dynamic MHD solver) on a grid. For the latter, KORC uses cubic splines interpolation. The parallelization of KORC was done using open multi-processing (MP) along with Message Passing Interface (MPI), resulting in excellent weak and strong scaling.

The organization of the rest of the paper is as follows. Section II describes the RE orbit equations of motion, the magnetic field model, and the numerical method. Section III focusses on full orbit effects on synchrotron radiation, which to the best of our knowledge has not been studied in detail in the fusion plasma physics literature before. We start with a study of spatial confinement of RE and its impact on synchrotron radiation. Following this, we discuss the statistics of collisionless pitch angle dispersion and its effect on synchrotron radiation. Of particular interest is the study of the role played by configuration-space orbit effects on the total radiation power and synchrotron spectra. This section also considers full orbit effects in the presence of radiation damping and electric field acceleration. Section IV discusses the magnetic field variability along RE orbits, questioning the validity of gyro-averaging in the relativistic regime. Section V presents the conclusions, and the Appendix presents details of the numerical method.

## II. FULL ORBIT MODEL AND NUMERICAL METHOD

Our starting point is the relativistic equations of motion

where ** x** and $ v $ denote the position and velocity of an electron of charge −

*e*and mass

*m*, and $ p = m e \gamma v $ is the relativistic momentum with $ \gamma = ( 1 \u2212 v 2 / c 2 ) \u2212 1 / 2 $. The first term on the right hand side of Eq. (2) is the Lorentz force, and $ F R $ is the radiation reaction force.

_{e}The most complete, physically consistent model of $ F R $ is the Landau-Lifshitz representation^{21} of the Lorentz-Abraham-Dirac force

where $ D / Dt = ( \u2202 \u2202 t + v \xb7 \u2207 ) $ is the advective derivative. The first term in curly brackets on the right hand side of Eq. (3) is negligible, and thus, it will not be included in the calculation. In particular, this term does not contribute to the radiation power, and its contribution to the reaction force is three to four orders of magnitude smaller than the dominant third term in curly brackets.

We assume a toroidal domain with major radius *R*_{0} and circular cross section with radius *r _{edge}*. The coordinates $ ( r , \u03d1 , \zeta ) $ are defined as $ x = ( R 0 + r \u2009 cos \u2009 \u03d1 ) \u2009 sin \u2009 \zeta , \u2009 y = ( R 0 + r \u2009 cos \u2009 \u03d1 ) \u2009 cos \u2009 \zeta $, and $ z = r \u2009 sin \u2009 \u03d1 $, where (

*x*,

*y*,

*z*) are the Cartesian coordinates. In these coordinates,

*r*denotes the minor radius, ϑ the poloidal angle, and

*ζ*the toroidal angle. Note that in this right-handed toroidal coordinate system, the toroidal angle

*ζ*rotates clockwise, that is, it is anti-parallel to the azimuthal angle, $ \varphi = \pi / 2 \u2212 \zeta $, of the standard cylindrical coordinate system.

The magnetic field model is

where $ \eta = r / R 0 $ is the aspect ratio, the constant *B*_{0} denotes the magnitude of the toroidal magnetic field, and $ B \u03d1 ( r ) = \eta B 0 / q ( r ) $ is the poloidal magnetic field with safety factor

The constant *q*_{0} is the safety factor at the magnetic axis, and the constant *λ* is determined from *q*_{0} and the value of *q*(*r*) at the plasma edge $ r = r e d g e $.

As it is well-known, when $ E = F R = 0 $, the relativistic energy, $ E = \gamma m e c 2 $, is a constant of motion. If in addition the magnetic field is axisymmetric, the canonical angular momentum $ p \zeta $, is also a constant of motion. For the magnetic field model in Eq. (4)

where $ \Omega e = e B 0 / \gamma m e $ is the electron relativistic cyclotron frequency at the magnetic axis. For the electric field, we adopt the following model

where the constant *E*_{0} denotes the electric field at the magnetic axis. As mentioned before, the fields accepted by the KORC can be either defined analytically (like in the case considered here) or pre-computed numerically (e.g., using VMEC, JFIT or any other equilibrium or dynamic MHD solver) on a grid. For the latter, the KORC uses cubic splines interpolation.

The full orbit dynamics of RE are computed using the code KORC that solves Eqs. (1)–(3) using a modified leap-frog method^{22} for advancing the electron's position ** x** and momentum

**under the Lorentz force, and the scheme in Ref. 23 for the radiation reaction force. A detailed description of the numerical method can be found in the Appendix. As expected, when choosing the time step $ \Delta t $, there is a trade-off between accuracy and computational time. For the calculations of interest in this paper, it was observed that $ \Delta t = \tau e / 100 $ (where $ \tau e = 2 \pi \gamma m e / e B 0 $ is the relativistic electron gyro-period) exhibits very good conservation properties without compromising the runtime. Most importantly, smaller time steps did not change the physics results reported in any significant way.**

*p*The numerical integration was tested by monitoring the conservation of energy, $E$, in the case $ E = F R = 0 $, and the conservation of canonical momentum, $ p \zeta $, in the case of toroidally symmetric magnetic fields. It was observed that the conservation properties of the integration scheme are comparable to previous studies resolving the full orbit dynamics of runaways.^{17,18,24} As a typical example, Fig. 1(a) shows $ \Delta E ( t ) / E 0 = [ E ( t ) \u2212 E 0 ] / E 0 $, with initial energy $ E 0 = 20 $ MeV and four different initial pitch angles $ \theta 0 = 10 \xb0 , 30 \xb0 , \u2009 50 \xb0 $ and $ 70 \xb0 $, where the pitch angle is defined as $ \theta = arctan ( v \u22a5 / v \u2225 ) $, with $ v \u22a5 $ and $ v \u2225 $ the perpendicular and the parallel velocity component, respectively. Figure 1(b) shows $ \Delta p \zeta ( t ) / p \zeta 0 = [ p \zeta ( t ) \u2212 p \zeta 0 ] / p \zeta 0 $ for the same simulation. The simulations followed runaway electrons in the magnetic field in Eq. (4) up to a time $ t = 1 \xd7 10 \u2212 3 $ s with $ \Delta t = 6.5 \xd7 10 \u2212 12 \u2009 s $. Excellent energy conservation, up to $ \xb1 10 \u2212 12 $, is observed in all cases, independently of the initial energy and position, pitch angle, and whether particles are passing or trapped. The momentum conservation for passing particles exhibits small oscillations varying from $ \Delta p \zeta / p \zeta 0 \u223c \xb1 10 \u2212 4 $, for small pitch angles $ \theta 0 \u2264 30 \xb0 $, to $ \Delta p \zeta / p \zeta 0 \u223c 5 \xd7 10 \u2212 3 $, for larger initial pitch angles $ \theta 0 > 30 \xb0 $. In the case of trapped particles, these oscillations are twice as big as for the case of passing particles. While energy conservation is not sensitive to the chosen time step, the amplitude of the oscillations of $ \Delta p \zeta / p \zeta 0 $ decreases linearly with the time step. In all cases studied, the behavior of the energy and canonical angular momentum did not exhibit a qualitative change with $ \Delta t $. In particular, the numerical errors did not show any secular growth, guaranteeing the fidelity and stability of the computations for large times. Simulations (not reported here) with non-axisymmetric (integrable and chaotic) magnetic fields show similar energy conservation properties. Of course, in these cases, the canonical angular momentum is not a constant of motion due to the lack of toroidal symmetry.

## III. FULL ORBIT EFFECTS ON SYNCHROTRON RADIATION

In this section, we study space dependent, full orbit effects on the synchrotron radiation of relativistic RE in toroidal geometry. In all the calculations we use the magnetic field model in Eq. (4) with DIII-D and ITER-like parameters. For DIII-D we set $ B 0 = 2.19 $ T, $ R 0 = 1.5 $ m, and minor radius *r _{edge}* = 0.5 m. Observations of RE disruptions in DIII-D

^{25}suggest that for these plasmas the safety factor at the plasma edge varies between

*q*= 2 and

_{edge}*q*= 3. We consider these two values for

_{edge}*q*and choose

_{edge}*q*

_{0}so that $ q e d g e / q 0 = 1.25 $ in both cases. In the case of ITER we set $ B 0 = 5.3 $ T, $ R 0 = 6.2 $,

*r*= 2.0 m,

_{edge}*q*= 3, and $ q e d g e / q 0 = 1.25 $. The use of different $ q e d g e / q 0 $ ratios do not seem to modify the results significantly.

_{edge}### A. Radial confinement and synchrotron radiation power

Although the spatial confinement of RE has been addressed before in the literature, e.g., Refs. 14, 15, and 26, we study it here because an accurate computation of synchrotron emission should only consider confined RE. Although this might seem an obvious statement, reduced models that do not track the spatial location of RE (e.g., bounced-averaged 2-D, $ ( p \u2225 , p \u22a5 ) $ phase space descriptions) face the risk of including RE that, because of their pitch angle and energy, might not be actually confined. In all our numerical computations of synchrotron radiation only confined RE in DIII-D and ITER-like plasmas are considered.

As mentioned in Sec. II, in the absence of an electric field, radiation damping, and collisions, the full orbit dynamics exactly conserves the toroidal canonical angular momentum $ p \zeta $, when the magnetic field is toroidally symmetric. The conservation of $ p \zeta $ can be used to gain valuable insights on the radial confinement of RE.^{26} To illustrate the basic idea consider the magnetic field model in Eq. (4) with a constant safety factor *q* = *q*_{0}. In this simple case, Eq. (6) provides the parametric representation $ ( R \u2212 \Delta ) 2 + Z 2 = R c 2 $, of the RE orbit in the (*R*, *Z*) cylindrical coordinates plane, where $ r = ( R \u2212 R 0 ) 2 + Z 2 , \u2009 \u2009 tan \u2009 \theta = Z / ( R \u2212 R 0 ) , \u2009 \Delta = R o \u2212 q 0 \upsilon \zeta / \Omega e , \u2009 \upsilon \zeta = R \zeta \u0307 $, and $ R c 2 = 2 q 0 p \zeta \Omega e \gamma m e + \Delta 2 \u2212 R 0 2 $. As customarily done in the literature, for the case of passing particles, it is tempting to interpret this parametric representation as a circle or radius *R _{c}* with the center shifted from the magnetic axis by Δ. However, care must be taken with this interpretation because $ \upsilon \zeta = R \zeta \u0307 $ is

*not*a constant of the RE orbit. As a result, strictly speaking,

*R*and Δ are not constant, and RE orbits are not simple shifted circles even in the case of passing particles. For example, passing RE with $ E 0 = 40 $ MeV and $ \theta 0 = 60 \xb0 $ in DIII-D-like magnetic fields exhibit variations in $ \zeta \u0307 $ in the range $ ( 0.4 c , 0.8 c ) $ which leads to $ \u223c 35 % $ variations in

_{c}*R*. To illustrate the complexity of the full orbit dynamics, Fig. 2 shows some typical passing and trapped RE orbits in DIII-D and ITER-like plasmas.

_{c}To study the dependence of confinement on the initial energy and pitch angle we consider ensembles of $ 4 \xd7 10 4 $ mono-energetic and mono-pitch angle RE initially distributed uniformly on a poloidal plane. In these simulations, the *q* profile depends on *r* according to Eq. (5). Recent studies^{28,29} indicate that the energy of runaway electrons in DIII-D tends to be in the 45 MeV–65 MeV range, with pitch angle $ \theta \u223c 10 \xb0 $. Motivated by this, the initial energy in the computations is in the range $ E 0 = 10 $ MeV–100 MeV, and the initial pitch angle varies from $ \theta 0 = 0 \xb0 $ to $ \theta 0 = 80 \xb0 $. The higher energies and larger pitch angles are included for theoretical extrapolation purposes. The integration time $ t \u223c 10 \u2212 6 $ s, corresponds approximately to ten toroidal turns, guaranteeing that the electrons undergo at least three poloidal transits. For the purpose of studying radial confinement in the non-stochastic magnetic field under consideration, there is no need to extend the total integration time. In this time scale, collisional transport is negligible, and thus, as discussed in the introduction, will not be taken into consideration. Although the electric field acceleration and the radiation damping can be easily incorporated in the calculation of radial confinement, we will assume $ E = F R = 0 $ because we are interested in radial confinement at fixed energy. At any rate, in the time scale of interest $ t \u223c 10 \u2212 6 $ s, the acceleration and damping are negligible. In the simulations, we stop following an electron once it crosses the plasma edge at $ r = r e d g e $, and keep advancing only those that remain confined.

The numerical results are reported by plotting the initial conditions (color coded by the initial value of the pitch angle, *θ*_{0}) in the poloidal plane corresponding to confined RE for a given initial energy $ E 0 $. Figure 3 compares a DIII-D-like case with *q _{edge}* = 2 with an ITER-like case, both with $ E = 40 $ MeV. For this energy, there are no confined trapped RE in DIII-D and thus only the confinement of passing electrons is plotted in panel (a). The situation is different in the ITER case that exhibits both passing RE, shown in panel (b), and trapped RE, shown in panel (c). Note that the regions of confinement are nested and shrink as the pitch angle increases (decreases) in the case of passing (trapped) RE.

A summary of the numerical results is reported in Fig. 4 that shows the percentage of RE that remain confined in the DIII-D (with *q _{edge}* = 2 and 3) and ITER like cases, for a wide range of initial energies and pitch angles. In the DIII-D case it is observed that increasing

*q*(i.e., decreasing the poloidal magnetic field) degrades the confinement. As expected, in both DIII-D and ITER like cases less confinement is observed when the energy increases. However, compared to DIII-D, the ITER case shows a very robust confinement with more than 50% of the runaways being confined for all the range of energies considered in the simulations. In all cases, the confinement properties remain nearly constant for a given energy and pitch angles in the range $ 0 \xb0 \u2264 \theta \u2264 40 \xb0 $, for larger pitch angles the confinement starts to decrease significantly.

_{edge}As shown in Fig. 4, RE with large energies ( $ E 0 \u2265 50 $ MeV) and large pitch angles ( $ \theta 0 > 40 \xb0 $) are poorly confined. However, because of their relatively large energies and pitch angles, these electrons might have a significant contribution to the total synchrotron radiation emitted, and thus it might be misleading to ignore them when accounting for the total measured power synchrotron emission. For each ensemble in the calculations reported in Fig. 4, we compute the total synchrotron radiated power $ P R = \u2211 j N P R j $, where

is the radiated power by the *j*-th confined electron of the ensemble, with ** B** and

**given by Eqs. (4) and (7), respectively. Figure 4 shows $ P R $ for ensembles of RE with different energies and pitch angles. It is observed that, for the DIII-D case with**

*E**q*= 2, the larger contribution to the total synchrotron radiation comes from confined runaway electrons with energies in the range of 50 MeV–80 MeV. For

_{edge}*q*= 3, the larger contribution to $ P R $ comes from confined electrons with energies in the range of 40 MeV–70 MeV. In both cases, the maximum of $ P R $ is observed at pitch angles of $ \theta \u223c 50 \xb0 $. It is important to note that for the inferred pitch angles in DIII-D plasmas

_{edge}^{25}$ \theta \u2248 10 \xb0 , \u2009 P R $ monotonically increases with $ E 0 $ for

*q*= 2. However, for

_{edge}*q*= 3 the total synchrotron radiated power first increases with the energy of the runaway electrons up to $ E 0 = 70 $ MeV, and then decreases for higher energies because confinement of runaway electrons degrades for high energies and large

_{edge}*q*. For the ITER case, the maximum of $ P R $ is at $ \theta \u223c 50 \xb0 $, and $ P R $ monotonically increases with $ E 0 $, this is mainly because the number of confined electrons does not vary greatly for the range of energies considered in our simulations.

_{edge}### B. Statistics of collisionless pitch angle dispersion

Full orbit effects imply that even in the absence of collisions, charge particles in toroidal geometry exhibit a dispersion in their pitch angle resulting from the spatial variations of the magnetic field. This effect, which is neglected in 2-D phase space Fokker-Planck models and treated approximately in guiding center models, has an impact on the accurate evaluation of synchrotron radiation emission. The goal of this subsection is to explore this problem in the context of full orbit computations. Collisionless pitch angle dispersion in RE has been recently discussed in Refs. 17 and 18. Here we go beyond this work by presenting a statistical study, and most importantly, by discussing the impact of collisionless pitch angle dispersion on synchrotron radiation.

The time evolution of the pitch angle $ \theta ( t ) $ depends not only on its initial value, *θ*_{0}, and the parameters of the magnetic field but also on the initial energy $ E 0 $ and initial spatial location of the electron. Thus, the precise computation of collisionless pitch angle dispersion requires solving the equations of motion in Eqs. (1) and (2) for an ensemble of confined runaway electrons with given initial values of $ E 0 $ and *θ*_{0}. An observable of particular interest is the steady state pitch angle marginal probability density function (PDF)

where $ f ( r , p ) $ denotes the 6-D PDF which is estimated using the full orbit computation of the RE ensemble. Figure 5 shows an example computed by following a large ensemble of initially mono-energetic and mono-pitch angle RE. By construction, $ f \theta ( \theta , t = 0 ) = \delta ( \theta \u2212 \theta 0 ) / \pi $. However, in these cases, even in the absence of electric field acceleration, radiation damping, and collisions, the asymptotic steady state PDF shows a significant spreading. Note that despite the similarity of the confinement regions of runaways with $ \theta 0 \u2264 50 \xb0 $ in this simulation, the spreading of $ f ( \theta ) $ depends strongly on *θ*_{0}. Figure 6 shows the standard deviation of the pitch angle PDF $ \sigma \theta = \u27e8 [ \theta \u2212 \u27e8 \theta \u27e9 ] 2 \u27e9 1 / 2 $, where $ \u27e8 a \u27e9 = \u222b a f ( \theta ) d \theta $, as a function of the RE energy for different values of the initial pitch angle. In all cases reported in the figure, it is observed that for a given energy, the spreading increases with the value of the initial pitch angle *θ*_{0}. On the other hand, for small initial pitch angles ( $ \theta 0 \u2264 30 \xb0 $ in the figure) $ \sigma \theta $ increases with the energy but for large initial pitch angles ( $ \theta 0 \u2265 40 \xb0 $ in the figure) $ \sigma \theta $ decreases with the energy.

As a first step to characterize $ f ( \theta ) $, Fig. 7 shows the numerically computed PDF (constructed from the histograms of the pitch angle data) for all the initial energies and all the initial pitch angle ensembles. To characterize the general functional form of $ f ( \theta ) $, the data is plotted using rescaled variables: the horizontal axis is shifted by *θ*_{0} and rescaled by $ 1 / \sigma \theta $ and the vertical axis is rescaled by $ \sigma \theta $. The different colored lines correspond to the PDFs for different values of the initial energy $ = 10 , 20 , \u2026 60 \u2009 MeV $ and the initial pitch angle $ \theta 0 = 10 \xb0 , \u2026 , 50 \xb0 $ resulting in thirty different cases. The observation that, to a good approximation, the data of all the cases collapse into a coherent trend indicates that

where $ \sigma \theta $ depends on $ E 0 $ and *θ*_{0} as indicated in Fig. 6, and *F* is a “universal” (i.e., independent of $ E 0 $, *θ*_{0}) function. As shown in Fig. 7, *F* is well approximated by a Von-Misses PDF of the form

where $ \alpha = 1.75 $ and $ I o $ denotes the modified Bessel function of zeroth order, included for normalization purposes. The parameter *α* is expected to depend on the details of the magnetic field and the aspect ratio. Note that in this calculation we have excluded data corresponding to $ \theta 0 = 0 \xb0 $. The reason being that, as shown in Fig. 5, initial condition with $ \theta 0 = 0 \xb0 $ are special because full orbit effect preclude electrons with vanishing pitch angle, and as a result, the PDF rapidly develops a non-generic structure.

### C. Effects of collisionless pitch angle dispersion on synchrotron radiation

A problem of particular interest is to explore the role of the observed collisionless pitch angle dispersion on synchrotron radiation. In the absence of electric field, Eq. (8) implies that the instantaneous total power radiated by an electron that at time *t* is at $ x ( t ) $ with velocity $ v ( t ) $ and pitch angle $ \theta ( t ) $ is

where $ B 2 = | B ( x ) | 2 $ and $ v 2 = | v ( t ) | 2 $. Let

denote the instantaneous curvature of the RE orbit with velocity $ r \u0307 $ and acceleration $ r \xa8 $. Neglecting the electric field, and using the Lorentz force

Using this expression, Eq. (12) can be equivalently written as

The geometry of the orbit enters in the evaluation of the curvature, and thus in the evaluation of the instantaneous radiated power, throughout the pitch angle dependence, $ \theta ( t ) $ and the magnetic field dependence on the electron's position $ B = | B ( x ( t ) ) | $. To explore these two dependences, Fig. 8 compares the curvature in Eq. (14) using the full orbit information with two approximate expressions: $ \kappa B 0 $ and $ \kappa B 0 \u200a \theta 0 $. The computation of $ \kappa B 0 $ takes into account the pitch angle dependence $ \theta = \theta ( t ) $, (which is taken from the full orbit integration) but neglects the spatial dependence of the magnetic field, i.e., it assumes $ B 2 = B 0 2 $ where *B*_{0} is the amplitude of the magnetic field at the magnetic axis. The more extreme approximation $ \kappa B 0 \u200a \theta 0 $ neglects the collisionless pitch angle dispersion, i.e., it assumes $ \theta = \theta 0 = $ constant, in addition to the assumption $ B 2 = B 0 2 $. As shown in Fig. 8, the temporal variability of *κ* is significantly reduced in the approximations $ \kappa B 0 $ and $ \kappa B 0 \u200a \theta 0 $. Both *κ* and $ \kappa B 0 $ exhibit a fast oscillation of the order of the gyro-period and slower oscillation of the order of the poloidal transit time.

To study in further detail the impact of full orbit effects on the computation of the curvature and the synchrotron radiation, we computed the total radiated power for an ensemble of $ N = 4 \xd7 10 4 $ RE electrons initially distributed uniformly on a poloidal plane with $ \theta 0 = 10 \xb0 $ and two values of the initial energy, $ E 0 = 40 $ MeV and $ E 0 = 80 $ MeV. Figure 9 shows the results of the computation according to Eq. (15) using the exact curvature, *κ*, based on the full orbit numerical integration (black dots), $ \kappa B 0 $ (green dots) and $ \kappa B 0 \u200a \theta 0 $ (red dot). As expected, when using $ \kappa B 0 \u200a \theta 0 $ the power radiated is the same for all the RE in the ensemble a result quite different to the exact calculation (black dots) that exhibits a significant variation resulting from the collisionless pitch angle dispersion shown in the histogram below the horizontal axes in the plots of Fig. 9. The corresponding histograms of the radiated power are shown in the vertical axes. It is important to reiterate that in these calculations there is no electric field acceleration, no radiation damping, and no collisions. Therefore, the observed changes in the curvature are solely due to the spatial effects of the RE orbit. Under these conditions, in the context Fokker-Planck 2-D phase space models, *κ* will remain constant and equal to $ \kappa B 0 \u200a \theta 0 $ which, as shown in Fig. 8, is a poor approximation of *κ* that leads to the underestimation of the radiated power observed Figure 9.

Figure 10 compares the ensemble average of *P _{R}* computed using the full orbit information with the approximate power radiated

*P*, neglecting all the spatial information, i.e., using $ \kappa B 0 \u200a \theta 0 $ for various ensembles of runaway electrons with different energies and initial pitch angles for DIII-D and ITER type magnetic fields. For the DIII-D case,

_{app}*P*seems to be a good approximation of $ \u27e8 P R \u27e9 $ only for runaway electrons with small energy $ E 0 \u2248 10 $ MeV. On the other hand, for the ITER case,

_{app}*P*is observed to be a fair approximation of $ \u27e8 P R \u27e9 $ for most energies and initial pitch angles, specially for the pitch angles $ \theta 0 = 10 \xb0 $ and $ 70 \xb0 $.

_{app}In addition to the total radiated power, the collisionless pitch angle dispersion modifies the radiation spectrum. The synchrotron power spectral density for a relativistic electron is given by^{27}

where $ \lambda c = ( 2 \lambda 0 / 3 ) ( m c 2 / E ) 3 $ is the critical wavelength, $ \lambda 0 = 2 \pi / \kappa $ and $ E = m 2 c 4 + c 2 p 2 $ are the total relativistic energies. As in the case of the total power, the spectrum depends on the geometry of the RE orbit because, as shown in Eq. (14), the curvature *κ* depends on the local value of the magnetic field $ B ( r ) $, at the instantaneous location of the RE, and on the instantaneous value of the pitch angle, *θ*. To explore the full orbit effects on the spectrum, Fig. 11 compares the evaluation of Eq. (16) using the exact RE curvature (i.e., computed using the full orbit information) with the approximate formula $ \kappa B 0 \theta 0 $, i.e., assuming *B* = *B*_{0} and $ \theta = \theta 0 $. In addition to the previously discussed underestimation of the total radiated power (i.e., the integrated area below the shown curves) it is observed that the full orbit effects introduce a shift of the spectrum towards low wavelengths. On the other hand, the two spectra seem to converge for high enough wavelengths. It is interesting to compare these results with those in Ref. 28 that showed that computing synchrotron radiation spectra assuming a mono-energetic RE ensemble with a prescribed pitch angle (rather than using a probability distribution function of $ p \u2225 $ and $ p \u22a5 $) can give misleading results.

### D. Full orbit effects in the presence of electric field acceleration and radiation damping

To conclude we consider the long term dynamics of runaway electrons including both synchrotron radiation losses and acceleration due to the toroidal electric field in Eq. (7) with $ E 0 = \u2212 4 $ V/m, a typical value observed in DIII-D experiments. The sign of *E*_{0} is chosen to guarantee that ** E** is anti-parallel to the toroidal magnetic field in order to accelerate runaway electrons streaming in the same direction as the magnetic field. We consider three ensembles of $ 3 \xd7 10 3 $ runaway electrons with an initial energy of $ E 0 = 40 $ MeV, and initial pitch angles of $ \theta 0 = 10 \xb0 , \u2009 30 \xb0 $ and $ 50 \xb0 $ in a DIII-D type magnetic field with

*q*= 2. For these parameters the balance between electric field acceleration and synchrotron radiation damping is approximately reached at the limit energies 146 MeV, 47 MeV and 25 MeV for runaways with $ \theta 0 = 10 \xb0 , \u2009 30 \xb0 $ and $ 50 \xb0 $, respectively. The total integration time was $ t = 10 \u2212 2 $ s, which is, the characteristic time scale of flat-top RE in DIII-D plasmas.

_{edge}^{25}To contrast the role played by the acceleration and damping, we performed simulations with radiation damping only, with electric field acceleration only, and with both radiation damping and electric field acceleration. To ensure that the electrons in these simulations are initially confined we chose initial conditions inside the dashed black circle in Fig. 3(a).

The results are summarized in Fig. 12 that shows the parallel $ p \u2225 $, and perpendicular $ p \u22a5 $, components of the momentum of each electron in the ensemble superimposed to energy and pitch angle iso-contours. We also show the marginal probability distribution functions of pitch angle, $ f ( \theta ) $ in Eq. (9), and the marginal probability density function of energy

where $ f ( r , p ) $ denotes the 6-D probability density function which is estimated using the full orbit computation of the ensemble of RE. The net effect of the synchrotron radiation damping is to decelerate the electrons without changing significantly the pitch angle distribution function $ f ( \theta ) $, shown in Fig. 12(b), with respect to the case where no radiation losses are included (cf. Fig. 5). We observe a large energy loss of up to 35% for electrons with $ \theta 0 = 50 \xb0 $, while for electrons with a small pitch angle $ \theta 0 = 10 \xb0 $ the average energy loss is very small, only 5%. In Fig. 12(c) we show the corresponding energy distribution functions $ f ( E ) $ at the final simulation time. Figures 12(d)–12(f) show the corresponding results for simulations with the electric field but no radiation damping. We observe a preferential acceleration of the runaway electrons along the parallel direction that generates a shift of $ f ( \theta ) $ towards smaller pitch angles and a reduction of its spreading, see Fig. 12(e). The average energy gain is approximately 25% for electrons with $ \theta 0 = 50 \xb0 $ and 30% for electrons with $ \theta 0 = 10 \xb0 $. Notice the very different shape of $ f ( E ) $ for this simulation with respect to Fig. 12(c). Finally, Figs. 12(g)–12(i) show results including both synchrotron radiation losses and electric field acceleration. In this case, the radiation losses and the electric field are so that electrons with $ \theta 0 \u2264 30 \xb0 $ keep being energized, while for electrons with $ \theta 0 = 50 \xb0 $ the energy losses are larger than the acceleration of the electric field. For these simulations, the pitch angle distribution function $ f ( \theta ) $ is similar to the second simulation where only the electric field is included, and the shape of the energy distribution functions $ f ( E ) $ are similar to those of the first simulation where only the radiation losses are included.

These results highlight the role of collisionless pitch angle dispersion in the otherwise straightforward (in the collisionless regime) tendency of the electric field to align the particle velocity to the magnetic field. In particular, as panels (b), (e) and (h) of Fig. 12 show, the orbit effects tend to frustrate the alignment. This also seems to indicate that the simple idea of an attractor (resulting from the balance of synchrotron emission, collisions, and acceleration) might require reconsideration to incorporate the role of the full orbit induced pitch angle dispersion.

## IV. MAGNETIC FIELD VARIATION

As discussed in Sec. I the use of orbit-averaged reduced models rests on the assumption that the magnetic field variations along the RE electron orbits is negligible. In particular, the guiding center approximation requires these variations to be small in the scale of the gyro-motion.^{17,30,31} Bounced averaged 2-D Fokker-Planck phase space descriptions are even more restrictive since they neglect the configuration space information resulting from the magnetic field and RE orbit geometry. The goal of this section is to explore numerically the validity these assumptions for DIII-D and ITER like magnetic fields. Going beyond the recent work in Ref. 18, we present a statistical study uncovering the full spatial variation of the magnitude and direction of the magnetic field for large ensembles of RE as the function of the energy and the pitch angle.

We start by considering the variability of the magnetic field that RE with different initial energies $ E 0 $ and pitch angles *θ*_{0} experience along their orbits. The simulations follow RE for one local gyro-period $ \tau e = 2 \pi \gamma m e / e B ( x ( t 0 ) ) $, where $ x ( t 0 ) $ denotes the initial condition, in the same magnetic fields used in the simulations reported in Sec. III. As before, we set the electric field and the radiation reaction force to zero because, for this time scale, they are negliglible. The initial energies of the RE are in the range $ 10 \u2009 MeV \u2264 E 0 \u2264 100 \u2009 MeV $, and the initial pitch angles in the range $ 0 \xb0 \u2264 \theta 0 \u2264 80 \xb0 $. For each set of $ E 0 $ and *θ*_{0} we uniformly distribute an ensemble of $ 4 \xd7 10 4 $ RE on a poloidal plane and follow them for a time *τ _{e}*. For each electron, we calculate the maximum change of the magnetic field amplitude

and magnetic field direction

where $ x ( t ) $ is the electron's position at time *t*, and $ max ( f ) $ denotes the maximum value of *f* between the times *t*_{0} and $ t 0 + \tau e $. Only confined (passing and trapped) electrons are taken into account in this calculation.

Figure 13 shows the spatial distribution of *δ _{B}* and $ \Delta B $ for DIII-D parameters with

*q*= 2 for the two sets of initial conditions $ E 0 = 60 $ MeV and $ \theta 0 = 0 \xb0 $, and $ E 0 = 60 $ MeV and $ \theta 0 = 40 \xb0 $. For small values of

_{edge}*θ*

_{0}there are well defined regions in space where the magnitude and direction of the magnetic field exhibit large changes, especially at the top and bottom of the plasma for

*δ*, and at the plasma boundary for $ \Delta B $. The regions of large

_{B}*δ*at the top and bottom of the plasma occur due to RE streaming from regions close to

_{B}*R*=

*R*

_{0}to regions of lower (higher) magnetic field at $ R > R 0 $ ( $ R < R 0 $), while the region of large $ \Delta B $ at the plasma boundary occurs mainly due to the radial transport of the electrons driven by the upward curvature drift $ v c u r v \u2248 \gamma m e v \u2225 2 / ( e B 0 R 0 ) $. Here $ v \u2225 $ is the parallel component of the RE velocity. As the initial pitch angle

*θ*

_{0}increases,

*δ*becomes larger and $ \Delta B $ becomes smaller and the regions of large

_{B}*δ*and $ \Delta B $ become less localized because the electrons can explore more regions of the space due to their larger gyro-radius (see panels (c) and (d) of Fig. 13).

_{B}A statistical summary of all the simulations performed for the *q _{edge}* = 3 DIII-D type magnetic field is presented in Fig. 14 that shows the ensemble averages of the magnitude $ \u27e8 \delta B \u27e9 $, and the direction, $ \u27e8 \Delta B \u27e9 $ of the magnetic field variation for a wide range of $ E 0 $ and

*θ*

_{0}. It is observed that the maximum change of the magnetic field amplitude along a RE orbit $ \u27e8 \delta B \u27e9 $, monotonically increases with both the initial energy $ E 0 $ and initial pitch angle

*θ*

_{0}. On the other hand, the maximum change in the direction of the magnetic field $ \u27e8 \Delta B \u27e9 $ increases with $ E 0 $ but decreases with increasing

*θ*

_{0}. The values of $ \u27e8 \delta B \u27e9 $ and $ \u27e8 \Delta B \u27e9 $ in the DIII-D case are observed to be independent of the value of

*q*, being as large as 15% and 60%, respectively. For ITER-type magnetic fields, similar trends are observed; however, the values of $ \u27e8 \delta B \u27e9 $ and $ \u27e8 \Delta B \u27e9 $ tend to be an order of magnitude smaller than for the DIII-D case.

_{edge}It is interesting to discuss the relation of these results with the Larmor radius and magnetic field length scales. In particular, for the results reported in Fig. 13, the perpendicular Larmor radius $ \rho L = v \u22a5 / \Omega e $, is of the order $ \rho L \u223c 0.02 $ m in (a)–(b), and of the order $ \rho L \u223c 0.10 $ m in (c)–(d). In both cases, *ρ _{L}* is certainly smaller than the length scale of variation of the magnetic field $ L B \u223c 1 $ m. However, the accurate evaluation of the variations of the magnetic field in one gyro-period requires considering, in addition to

*ρ*, the parallel Larmor radius length scale $ \rho \u2225 = v \u2225 \tau e = 2 \pi v \u2225 / \Omega e $. For the results reported in Fig. 13, independently of the pitch angle $ \rho \u2225 \u223c 0.6 $ m, which is significantly larger than

_{L}*ρ*and comparable to

_{L}*L*. That is, the large parallel displacements experienced by electrons in the relativistic regime might signal the breakdown of reduced models based on gyro-averaging despite the fact that the perpendicular Larmor radius might be small.

_{B}## V. SUMMARY AND CONCLUSIONS

The majority of previous works on RE have been based on bounce-averaged Fokker-Planck models determining the evolution of the RE distribution function in a reduced 2-D phase space, e.g., $ p | | $ and $ p \u22a5 $. Although these reduced descriptions have provided very valuable insights, there is the pressing need to incorporate configuration space dynamical information. Motivated by this, in this paper we have presented numerical results on full orbit effects on the dynamics or RE in magnetically confined fusion plasmas in toroidal geometry. Our study is based on the direct numerical integration of large ensembles of RE fully resolving the gyro-motion and using the exact expression of the Lorentz force and the Landau-Lifshitz consistent formulation of the Abraham-Lorentz-Dirac radiation reaction force. To cope with the wide gap between the electron gyro-period time scale, $ \u223c 10 \u2212 11 $ s, and the time scales of interest, $ \u223c 10 \u2212 2 $ s we have used an efficient method that exhibits excellent conservation properties and stability.

Of particular interest has been the impact of full orbit effects on synchrotron radiation. In the relativistic regime of interest to RE in fusion plasmas, synchrotron radiation is largely determined by the curvature, *κ*, of the RE orbit. In general, *κ* depends on the instantaneous value of the pitch angle and the local magnetic field at the spatial location of the RE. Reduced descriptions that ignore the spatial (configuration-space) degrees of freedom (e.g., 2-D phase space Fokker-Planck models) cannot account for spatial variations of the magnetic field along the RE orbit and are thus forced to assume a constant magnetic field (usually taken to be the value at the magnetic axis). This restrictive assumption does not allow the exploration of the potentially highly nontrivial effects that the 3-D magnetic field geometry can have on synchrotron radiation. Here we have shown that even in the relatively simple case of an integrable magnetic field (i.e., even in the absence of magnetic field stochasticity) the constant magnetic field assumption can lead to a significant underestimation of synchrotron radiation.

In addition to the local magnetic field, synchrotron radiation has a key dependence on the pitch angle *θ*. The role of collisions and electric field acceleration on *θ* are well captured in 2-D phase space descriptions. Reduced models can also account for the (much smaller) changes in *θ* due to radiation damping. However, the role of the less known collisionless pitch angle dispersion requires full orbit information and is outside the scope of reduced descriptions. Collisionless pitch angle dispersion results from the fact that *θ* is not a constant of motion and therefore it changes as a function of time even in the absence of collisions, radiation damping, and electric field acceleration. To illustrate this effect, we considered ensembles of initially mono-pitch angle RE, i.e., such that $ f ( \theta , t = 0 ) = \delta ( \theta \u2212 \theta 0 ) / \pi $. In the absence of collisions, radiation damping, and electric fields, reduced models trivially imply the invariance of the pitch angle distribution, i.e., $ f ( \theta , t > 0 ) = \delta ( \theta \u2212 \theta 0 ) / \pi $. However, full orbit numerical results show that the toroidal geometry of the 3-D magnetic field gives rise to the spreading of *θ* that leads to a stationary pitch angle probability distribution function with finite variance. Most importantly, the numerical results show that the collisionless dispersion of the pitch angle plays a nontrivial role in both the total radiation power and the radiation spectrum.

A topic of interest in the recent literature has been the study of the energy limit of RE resulting from the time-asymptotic equilibration of electric field acceleration, radiation damping, and Coulomb drag. Although the existence of phase space “attractors” describing such asymptotic states is a plausible assumption in the context of reduced models, care must be taken when spatial degrees of freedom are included in the description. In particular, it is not obvious that the projection of the full 6-D dynamics on a reduced 2-D phase space plane will result on a simple attractor. Our numerical simulations in the presence of electric field acceleration and radiation damping indicate that full orbit effects, and in particular collisionless pitch angle dispersion, have a nontrivial effects on the evolution of the pitch angle and energy distributions. This motivates the need for a detailed study of the asymptotic dynamics in the energy limit regime and the existence of attractors based on full orbit information. We plan to address these problems in a future publication.

The use of a simplified electric field model in the calculations presented here is motivated by the fact that for the time scales of the calculations presented (e.g., instantaneous synchrotron radiation and colissionless pitch angle dispersion) the details of the electric field temporal evolution are not relevant. However, for problems tracking the electric field acceleration over long times (e.g., acceleration from thermal energies) an accurate modeling of the time dependence of the loop voltage is required. This is a problem that we plan to address in a future publication.

Our main goal has been to study the role played by configuration-space orbit effects on the dynamics of relativistic runaway electrons, with particular emphasis on synchrotron radiation. These effects can be studied using the exact equations of motion or by using approximate descriptions of the orbits like those provided by the guiding center model. Faced with the pressing need to develop codes with predictive capabilities, we decided to avoid the approximate guiding center model, and perform full orbit calculations. It is important to realize that with the current hardware computing capabilities and advanced numerical methods full orbit calculations like those discussed in this paper are quite accessible in terms of runtime and accuracy. In the relativistic regime of interest, synchrotron radiation is dominated by the gyro-motion which is not resolved in the guiding center description. Another problem of interest demanding the full resolution of gyro-motion is the interaction of the magnetic ripple with the RE beam. At a more fundamental level, we have shown that the variations of the magnetic field in one gyro-period can be significant, thus questioning the validity of gyro-averaging (the cornerstone of the guiding center model) in the study of relativistic runaway electrons. Due to the machine size and magnitude of the magnetic field, these variations are more significant in DIII-D-like plasmas that in ITER-like plasmas. Having said this, we want to indicate that a detailed comparison between the guiding center and the full-orbit physics predictions and computational runtime for the problems discussed here has not been done and it is outside the scope of this paper.

## ACKNOWLEDGMENTS

We thank Jeremy Lore for the valuable comments on the manuscript. Research sponsored by the Laboratory Directed Research and Development Program of Oak Ridge National Laboratory, managed by UT-Battelle, LLC, for the U.S. Department of Energy. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.

This manuscript has been authored by UT-Battelle, LLC under Contract No. DE-AC05-00OR22725 with the U.S. Department of Energy. The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a non-exclusive, paid-up, irrevocable, worldwide license to publish or reproduce the published form of this manuscript, or allow others to do so, for United States Government purposes. The Department of Energy will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan (http://energy.gov/downloads/doe-public-access-plan).

### APPENDIX: NUMERICAL METHOD

KORC uses a numerical scheme based on a combination of two numerical methods that results in an accurate, stable, and easily parallelizable algorithm for solving the dynamics of runaway electrons. The first method is the modified relativistic leap-frog algorithm of Ref. 22 which is a conservative and stable technique for advancing the electron's position ** x** and momentum

**in Eqs. (1)–(2) using only the relativistic Lorentz force. This method does not include numerical spurious forces that detriment the conservation properties of other leap-frog methods such as the relativistic Boris and corrected relativistic Boris methods. The second method is the one proposed in Ref. 23. We use this method to include the radiation reaction force in Eq. (3) keeping intact the conservation and stability properties of the modified relativistic leap-frog method.**

*p*The variable ** x** is computed at the half time steps $ i \u2212 1 / 2 , i + 1 / 2 , \u2026 $, while

**is computed at the integer time steps $ i , i + 1 , \u2026 $. In what follows, we will distinguish the momentum calculated at the time level**

*p**i*+ 1 including both the Lorentz force and the radiation reaction force as $ p i + 1 $, while $ p L i + 1 $ and $ p R i + 1 $ will denote the momentum calculated independently using the Lorentz force alone and the radiation reaction force alone, respectively.

The discretized equations of motion to advance the change in the momentum due to the Lorentz force are

where $ \Delta t $ is the time step, *q* denotes the charge, $ p j = m e \gamma j v j $, and $ \gamma j = 1 / 1 + v j 2 / c 2 $. The evolution of the relativistic *γ* factor is given by $ \gamma i + 1 = 1 + ( p L i + 1 / m e c ) 2 = 1 + p L i + 1 \xb7 p \u2032 / m e 2 c 2 $, which together with Eq. (A2), implies

where we have defined

$ \tau = ( q \Delta t / 2 ) B i + 1 / 2 , \u2009 t = \tau / \gamma i + 1 , \u2009 p * = p \u2032 \xb7 \tau / m e c , \u2009 \sigma = \gamma \u2032 2 \u2212 \tau 2 , \u2009 \gamma \u2032 = 1 + p \u2032 2 / m e 2 c 2 $, and $ s = 1 / ( 1 + t 2 ) $.

The discretized equation of motion to advance the change in the momentum due to the radiation reaction force force is

where $ p i + 1 / 2 = ( p L i + 1 + p i ) / 2 $. Finally, using $ p L i + 1 $ from Eq. (A3) and $ p R i + 1 $ from Eq. (A6) the momentum at time level *i* + 1 is given by