The dynamics of RE (runaway electrons) in fusion plasmas span a wide range of temporal scales, from the fast gyro-motion, 10 11 s, to the observational time scales, 10 2 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.

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 10 11 s to the observational time scales t 10 3 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 formulation21 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 τ c o l l 10 2 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.

Our starting point is the relativistic equations of motion

d x d t = v ,
(1)
d p d t = e [ E ( x ) + v × B ( x ) ] + F R ,
(2)

where x and v denote the position and velocity of an electron of charge −e and mass me, and p = m e γ v is the relativistic momentum with γ = ( 1 v 2 / c 2 ) 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.

The most complete, physically consistent model of F R is the Landau-Lifshitz representation21 of the Lorentz-Abraham-Dirac force

F R ( x , v ) = e 3 6 π ϵ 0 m e c 3 { γ [ D E Dt + v × D B Dt ] + e m e [ ( E · v ) c 2 E + ( E + v × B ) × B ] e γ 2 m e c 2 [ ( E + v × B ) 2 ( E · v ) 2 c 2 ] v } ,
(3)

where D / Dt = ( t + v · ) 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 R0 and circular cross section with radius redge. The coordinates ( r , ϑ , ζ ) are defined as x = ( R 0 + r cos ϑ ) sin ζ , y = ( R 0 + r cos ϑ ) cos ζ , and z = r sin ϑ , 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, ϕ = π / 2 ζ , of the standard cylindrical coordinate system.

The magnetic field model is

B ( r , ϑ ) = 1 1 + η cos ϑ [ B 0 e ̂ ζ + B ϑ ( r ) e ̂ ϑ ] ,
(4)

where η = r / R 0 is the aspect ratio, the constant B0 denotes the magnitude of the toroidal magnetic field, and B ϑ ( r ) = η B 0 / q ( r ) is the poloidal magnetic field with safety factor

q ( r ) = q 0 ( 1 + r 2 λ 2 ) .
(5)

The constant q0 is the safety factor at the magnetic axis, and the constant λ is determined from q0 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 = γ m e c 2 , is a constant of motion. If in addition the magnetic field is axisymmetric, the canonical angular momentum p ζ , is also a constant of motion. For the magnetic field model in Eq. (4)

p ζ = γ m e { R 0 2 ( 1 + η cos ϑ ) 2 ζ ̇ + λ 2 Ω e 2 q 0 log ( 1 + r 2 λ 2 ) } ,
(6)

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

E = E 0 1 + η cos ϑ e ̂ ζ ,
(7)

where the constant E0 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 method22 for advancing the electron's position x and momentum p 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 Δ t , there is a trade-off between accuracy and computational time. For the calculations of interest in this paper, it was observed that Δ t = τ e / 100 (where τ e = 2 π γ 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.

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 ζ , 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 Δ E ( t ) / E 0 = [ E ( t ) E 0 ] / E 0 , with initial energy E 0 = 20  MeV and four different initial pitch angles θ 0 = 10 ° , 30 ° , 50 ° and 70 ° , where the pitch angle is defined as θ = arctan ( v / v ) , with v and v the perpendicular and the parallel velocity component, respectively. Figure 1(b) shows Δ p ζ ( t ) / p ζ 0 = [ p ζ ( t ) p ζ 0 ] / p ζ 0 for the same simulation. The simulations followed runaway electrons in the magnetic field in Eq. (4) up to a time t = 1 × 10 3  s with Δ t = 6.5 × 10 12 s . Excellent energy conservation, up to ± 10 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 Δ p ζ / p ζ 0 ± 10 4 , for small pitch angles θ 0 30 ° , to Δ p ζ / p ζ 0 5 × 10 3 , for larger initial pitch angles θ 0 > 30 ° . 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 Δ p ζ / p ζ 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 Δ 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.

FIG. 1.

Time evolution of the change of relativistic energy, Δ E ( t ) / E 0 , and canonical angular momentum, Δ p ζ ( t ) / p ζ 0 , of passing (P) and trapped (T) runaways in typical KORC simulations. Panel (a): energy conservation for runaway electrons with initial energy E 0 = 20  MeV and four different initial pitch angles: θ 0 = 10 ° , 30 ° , 50 ° and 70 ° . In all cases the fluctuations of Δ E ( t ) / E 0 are observed to be between ± 10 12 . Panel (b): conservation of canonical angular momentum of runaway electrons in panel (a). The amplitude of the oscillations of Δ p ζ ( t ) / p ζ 0 for passing runaways varies between 10 4 and 5 × 10 3 , being larger for runaways with larger pitch angles. For trapped runaways, Δ p ζ ( t ) / p ζ 0 is two times larger than for passing particles.

FIG. 1.

Time evolution of the change of relativistic energy, Δ E ( t ) / E 0 , and canonical angular momentum, Δ p ζ ( t ) / p ζ 0 , of passing (P) and trapped (T) runaways in typical KORC simulations. Panel (a): energy conservation for runaway electrons with initial energy E 0 = 20  MeV and four different initial pitch angles: θ 0 = 10 ° , 30 ° , 50 ° and 70 ° . In all cases the fluctuations of Δ E ( t ) / E 0 are observed to be between ± 10 12 . Panel (b): conservation of canonical angular momentum of runaway electrons in panel (a). The amplitude of the oscillations of Δ p ζ ( t ) / p ζ 0 for passing runaways varies between 10 4 and 5 × 10 3 , being larger for runaways with larger pitch angles. For trapped runaways, Δ p ζ ( t ) / p ζ 0 is two times larger than for passing particles.

Close modal

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 redge = 0.5 m. Observations of RE disruptions in DIII-D25 suggest that for these plasmas the safety factor at the plasma edge varies between qedge = 2 and qedge = 3. We consider these two values for qedge and choose q0 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 , redge = 2.0 m, qedge = 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.

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 , p ) 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 ζ , when the magnetic field is toroidally symmetric. The conservation of p ζ 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 = q0. In this simple case, Eq. (6) provides the parametric representation ( R Δ ) 2 + Z 2 = R c 2 , of the RE orbit in the (R, Z) cylindrical coordinates plane, where r = ( R R 0 ) 2 + Z 2 , tan θ = Z / ( R R 0 ) , Δ = R o q 0 υ ζ / Ω e , υ ζ = R ζ ̇ , and R c 2 = 2 q 0 p ζ Ω e γ m e + Δ 2 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 Rc with the center shifted from the magnetic axis by Δ. However, care must be taken with this interpretation because υ ζ = R ζ ̇ is not a constant of the RE orbit. As a result, strictly speaking, Rc 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 θ 0 = 60 ° in DIII-D-like magnetic fields exhibit variations in ζ ̇ in the range ( 0.4 c , 0.8 c ) which leads to 35 % variations in Rc. 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.

FIG. 2.

Typical orbits, projected on a poloidal plane, of passing and trapped runaway electrons in the magnetic field of Eq. (4) with DIII-D and ITER like parameters, and E = 40  MeV. Panel (a) shows examples of confined (red trace) and unconfined (blue and grey traces) runaway electrons in a DIII-D-like magnetic field. The red cross marks the position of the magnetic axis, and the black crosses mark the initial position of the runaway electrons. Panel (b) shows examples of confined runaway electrons for an ITER-like magnetic field. In this case, most of the electrons remain confined. For reference, the figures also show the confinement regions for θ 0 = 50 ° and θ 0 = 60 ° (see Fig. 3).

FIG. 2.

Typical orbits, projected on a poloidal plane, of passing and trapped runaway electrons in the magnetic field of Eq. (4) with DIII-D and ITER like parameters, and E = 40  MeV. Panel (a) shows examples of confined (red trace) and unconfined (blue and grey traces) runaway electrons in a DIII-D-like magnetic field. The red cross marks the position of the magnetic axis, and the black crosses mark the initial position of the runaway electrons. Panel (b) shows examples of confined runaway electrons for an ITER-like magnetic field. In this case, most of the electrons remain confined. For reference, the figures also show the confinement regions for θ 0 = 50 ° and θ 0 = 60 ° (see Fig. 3).

Close modal

To study the dependence of confinement on the initial energy and pitch angle we consider ensembles of 4 × 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 studies28,29 indicate that the energy of runaway electrons in DIII-D tends to be in the 45 MeV–65 MeV range, with pitch angle θ 10 ° . 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 θ 0 = 0 ° to θ 0 = 80 ° . The higher energies and larger pitch angles are included for theoretical extrapolation purposes. The integration time t 10 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 10 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 qedge = 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.

FIG. 3.

Confinement regions for runaway electrons with E 0 = 40  MeV and different initial pitch angles. Panel (a): passing runaway electrons for DIII-D case with qedge = 2. Panel (b): passing runaway electrons for ITER case. Panel (c): same as panel (b) for trapped particles. Maps with similar shapes and smaller confining regions are obtained for runaway electrons with higher energies.

FIG. 3.

Confinement regions for runaway electrons with E 0 = 40  MeV and different initial pitch angles. Panel (a): passing runaway electrons for DIII-D case with qedge = 2. Panel (b): passing runaway electrons for ITER case. Panel (c): same as panel (b) for trapped particles. Maps with similar shapes and smaller confining regions are obtained for runaway electrons with higher energies.

Close modal

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 qedge = 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 qedge (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 ° θ 40 ° , for larger pitch angles the confinement starts to decrease significantly.

FIG. 4.

Percentage of confined runaway electrons (top row) and total synchrotron radiated power due to confined electrons (bottom row), for DIII-D (with qedge = 2 and qedge = 3) and ITER-like magnetic fields.

FIG. 4.

Percentage of confined runaway electrons (top row) and total synchrotron radiated power due to confined electrons (bottom row), for DIII-D (with qedge = 2 and qedge = 3) and ITER-like magnetic fields.

Close modal

As shown in Fig. 4, RE with large energies ( E 0 50  MeV) and large pitch angles ( θ 0 > 40 ° ) 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 = j N P R j , where

P R j = e 4 6 π ϵ 0 m e 2 c 3 { E 2 + ( v × B ) · E γ 2 [ ( E + v × B ) 2 ( E · v ) 2 c 2 ] } ,
(8)

is the radiated power by the j-th confined electron of the ensemble, with B and E 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 qedge = 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 qedge = 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 θ 50 ° . It is important to note that for the inferred pitch angles in DIII-D plasmas25  θ 10 ° , P R monotonically increases with E 0 for qedge = 2. However, for qedge = 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 qedge. For the ITER case, the maximum of P R is at θ 50 ° , 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.

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 θ ( 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)

f ( θ ) = d r d p δ ( θ θ ) f ( r , p ) ,
(9)

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 θ ( θ , t = 0 ) = δ ( θ θ 0 ) / π . 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 θ 0 50 ° in this simulation, the spreading of f ( θ ) depends strongly on θ0. Figure 6 shows the standard deviation of the pitch angle PDF σ θ = [ θ θ ] 2 1 / 2 , where a = a f ( θ ) d θ , 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 ( θ 0 30 ° in the figure) σ θ increases with the energy but for large initial pitch angles ( θ 0 40 ° in the figure) σ θ decreases with the energy.

FIG. 5.

Steady state pitch angle probability distribution function for different values of θ0 and E 0 = 40  MeV. In all cases, the initial probability distribution function was a delta function, f ( θ , t = 0 ) = δ ( θ θ 0 ) / π , and the observed spreading in the (time asymptotic) steady state results from collisionless pitch angle dispersion.

FIG. 5.

Steady state pitch angle probability distribution function for different values of θ0 and E 0 = 40  MeV. In all cases, the initial probability distribution function was a delta function, f ( θ , t = 0 ) = δ ( θ θ 0 ) / π , and the observed spreading in the (time asymptotic) steady state results from collisionless pitch angle dispersion.

Close modal
FIG. 6.

Standard deviation of steady state pitch angle probability distribution as function of energy for different values of the initial pitch angle.

FIG. 6.

Standard deviation of steady state pitch angle probability distribution as function of energy for different values of the initial pitch angle.

Close modal

As a first step to characterize f ( θ ) , 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 ( θ ) , the data is plotted using rescaled variables: the horizontal axis is shifted by θ0 and rescaled by 1 / σ θ and the vertical axis is rescaled by σ θ . The different colored lines correspond to the PDFs for different values of the initial energy = 10 , 20 , 60 MeV and the initial pitch angle θ 0 = 10 ° , , 50 ° 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

f ( θ ) = σ θ 1 F [ θ θ 0 σ θ ] ,
(10)

where σ θ 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

F ( x ) = e α cos x 2 π I o ( α ) ,
(11)

where α = 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 θ 0 = 0 ° . The reason being that, as shown in Fig. 5, initial condition with θ 0 = 0 ° are special because full orbit effect preclude electrons with vanishing pitch angle, and as a result, the PDF rapidly develops a non-generic structure.

FIG. 7.

Rescaled probability distribution functions of pitch angle for initial energies E 0 = 10 , 20 , 60 MeV and initial pitch angles θ 0 = 10 ° , , 50 ° resulting in thirty different ensembles of RE. The collapse of the data into a coherent trend lends support to the ansatz in Eq. (10) with F given in Eq. (11) and plotted with the solid black line.

FIG. 7.

Rescaled probability distribution functions of pitch angle for initial energies E 0 = 10 , 20 , 60 MeV and initial pitch angles θ 0 = 10 ° , , 50 ° resulting in thirty different ensembles of RE. The collapse of the data into a coherent trend lends support to the ansatz in Eq. (10) with F given in Eq. (11) and plotted with the solid black line.

Close modal

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 θ ( t ) is

P R j = e 4 B 2 v 2 6 π ϵ 0 m e 2 c 3 γ 2 sin 2 θ ,
(12)

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

κ = | r ̇ × r ¨ | | r ̇ | 3 ,
(13)

denote the instantaneous curvature of the RE orbit with velocity r ̇ and acceleration r ¨ . Neglecting the electric field, and using the Lorentz force

κ = e γ m e v 3 | v × ( v × B ) | = e B γ m e v sin θ .
(14)

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

P R j = e 2 6 π ϵ 0 c 3 γ 4 v 4 κ 2 .
(15)

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, θ ( 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: κ B 0 and κ B 0 θ 0 . The computation of κ B 0 takes into account the pitch angle dependence θ = θ ( 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 B0 is the amplitude of the magnetic field at the magnetic axis. The more extreme approximation κ B 0 θ 0 neglects the collisionless pitch angle dispersion, i.e., it assumes θ = θ 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 κ B 0 and κ B 0 θ 0 . Both κ and κ B 0 exhibit a fast oscillation of the order of the gyro-period and slower oscillation of the order of the poloidal transit time.

FIG. 8.

Comparison of different models for the instantaneous curvature of an electron with E 0 = 40  MeV and θ 0 = 10 ° . The approximated curvature κ B 0 (green trace) underestimates the actual curvature κ (black trace) of the electron orbit. Shown in red is the constant curvature κ B 0 θ 0 which neglects collisionless pitch angle dispersion and the spatial variations of the magnetic field.

FIG. 8.

Comparison of different models for the instantaneous curvature of an electron with E 0 = 40  MeV and θ 0 = 10 ° . The approximated curvature κ B 0 (green trace) underestimates the actual curvature κ (black trace) of the electron orbit. Shown in red is the constant curvature κ B 0 θ 0 which neglects collisionless pitch angle dispersion and the spatial variations of the magnetic field.

Close modal

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 × 10 4 RE electrons initially distributed uniformly on a poloidal plane with θ 0 = 10 ° 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), κ B 0 (green dots) and κ B 0 θ 0 (red dot). As expected, when using κ B 0 θ 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 κ B 0 θ 0 which, as shown in Fig. 8, is a poor approximation of κ that leads to the underestimation of the radiated power observed Figure 9.

FIG. 9.

Full orbit effects on synchrotron radiation power in ensembles of runaway electrons with θ 0 = 10 ° and E 0 = 40  MeV (top panel) and θ 0 = 10 ° and E 0 = 80  MeV (bottom panel). The black dots show the radiated power per electron, PR in Eq. (15), calculated with the curvature formula in Eq. (14) including collisionless pitch-angle dispersion and spatial variations of the magnetic field using the full orbit simulation data. The green dots show PR computed using κ B 0 , i.e., neglecting the spatial variations of the magnetic field. The red point denotes the result using κ B 0 θ 0 , i.e., neglecting the spatial variations of the magnetic field and colissionless pitch angle dispersion. The bars in the horizontal (vertical) axis denote the histogram of pitch angles (radiated power).

FIG. 9.

Full orbit effects on synchrotron radiation power in ensembles of runaway electrons with θ 0 = 10 ° and E 0 = 40  MeV (top panel) and θ 0 = 10 ° and E 0 = 80  MeV (bottom panel). The black dots show the radiated power per electron, PR in Eq. (15), calculated with the curvature formula in Eq. (14) including collisionless pitch-angle dispersion and spatial variations of the magnetic field using the full orbit simulation data. The green dots show PR computed using κ B 0 , i.e., neglecting the spatial variations of the magnetic field. The red point denotes the result using κ B 0 θ 0 , i.e., neglecting the spatial variations of the magnetic field and colissionless pitch angle dispersion. The bars in the horizontal (vertical) axis denote the histogram of pitch angles (radiated power).

Close modal

Figure 10 compares the ensemble average of PR computed using the full orbit information with the approximate power radiated Papp, neglecting all the spatial information, i.e., using κ B 0 θ 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, Papp seems to be a good approximation of P R only for runaway electrons with small energy E 0 10  MeV. On the other hand, for the ITER case, Papp is observed to be a fair approximation of P R for most energies and initial pitch angles, specially for the pitch angles θ 0 = 10 ° and 70 ° .

FIG. 10.

Comparison of numerically computed, full orbit, radiation power, P R , and approximate radiation power Papp. Top row shows the the ratio, P R / P a p p and bottom row the difference P R P a p p for runaway electrons with different initial energies and initial pitch angles. Left column DIII-D-like case with q e d g e = 2 , right column ITER-like case.

FIG. 10.

Comparison of numerically computed, full orbit, radiation power, P R , and approximate radiation power Papp. Top row shows the the ratio, P R / P a p p and bottom row the difference P R P a p p for runaway electrons with different initial energies and initial pitch angles. Left column DIII-D-like case with q e d g e = 2 , right column ITER-like case.

Close modal

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 by27 

P R = 1 3 c e 2 ε 0 γ 2 λ 3 λ c / λ K 5 / 3 ( η ) d η ,
(16)

where λ c = ( 2 λ 0 / 3 ) ( m c 2 / E ) 3 is the critical wavelength, λ 0 = 2 π / κ 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 κ B 0 θ 0 , i.e., assuming B = B0 and θ = θ 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 and p ) can give misleading results.

FIG. 11.

Comparison of the orbit averaged (dashed line) and the full orbit (solid line) computation of the synchrotron power spectrum for 40 MeV runaway electron in the case of DIII-D-like (left) and ITER-like (right) magnetic fields. In panels (a) and (d) θ 0 = 20 ° ; panels (b) and (e) θ 0 = 30 ° ; and panels (c) and (f) θ 0 = 50 ° .

FIG. 11.

Comparison of the orbit averaged (dashed line) and the full orbit (solid line) computation of the synchrotron power spectrum for 40 MeV runaway electron in the case of DIII-D-like (left) and ITER-like (right) magnetic fields. In panels (a) and (d) θ 0 = 20 ° ; panels (b) and (e) θ 0 = 30 ° ; and panels (c) and (f) θ 0 = 50 ° .

Close modal

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 = 4 V/m, a typical value observed in DIII-D experiments. The sign of E0 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 × 10 3 runaway electrons with an initial energy of E 0 = 40  MeV, and initial pitch angles of θ 0 = 10 ° , 30 ° and 50 ° in a DIII-D type magnetic field with qedge = 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 θ 0 = 10 ° , 30 ° and 50 ° , respectively. The total integration time was t = 10 2 s, which is, the characteristic time scale of flat-top RE in DIII-D plasmas.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 , and perpendicular p , 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 ( θ ) in Eq. (9), and the marginal probability density function of energy

f ( E ) = d r d p δ ( E E ) f ( r , p ) ,
(17)

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 ( θ ) , 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 θ 0 = 50 ° , while for electrons with a small pitch angle θ 0 = 10 ° 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 ( θ ) towards smaller pitch angles and a reduction of its spreading, see Fig. 12(e). The average energy gain is approximately 25% for electrons with θ 0 = 50 ° and 30% for electrons with θ 0 = 10 ° . 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 θ 0 30 ° keep being energized, while for electrons with θ 0 = 50 ° the energy losses are larger than the acceleration of the electric field. For these simulations, the pitch angle distribution function f ( θ ) 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.

FIG. 12.

Summary of results of simulations of runaway electrons including synchrotron energy losses and a toroidal electric field. In these simulations, three ensembles of runaway electrons are advanced for t = 10 ms in a magnetic field with characteristic parameters of DIII-D with qedge = 2. All the electrons are initialized with an energy of E 0 = 40  MeV, and initial pitch angles of θ 0 = 10 ° , 30 ° and 50 ° . Panel (a): parallel p and perpendicular p components of the linear momentum of each electron with respect to the local magnetic field at the final simulation time for the simulation including only radiation losses. Panel (b): pitch angle distribution functions f ( θ ) at the final simulation time for the simulation of panel (a). Panel (c): energy distribution functions f ( E ) at the final simulation time. Panels (d)–(f) show the corresponding results for the simulation including only the toroidal electric field. Panels (g)–(i) show the results of the simulation including both the synchrotron radiation losses and the toroidal electric field.

FIG. 12.

Summary of results of simulations of runaway electrons including synchrotron energy losses and a toroidal electric field. In these simulations, three ensembles of runaway electrons are advanced for t = 10 ms in a magnetic field with characteristic parameters of DIII-D with qedge = 2. All the electrons are initialized with an energy of E 0 = 40  MeV, and initial pitch angles of θ 0 = 10 ° , 30 ° and 50 ° . Panel (a): parallel p and perpendicular p components of the linear momentum of each electron with respect to the local magnetic field at the final simulation time for the simulation including only radiation losses. Panel (b): pitch angle distribution functions f ( θ ) at the final simulation time for the simulation of panel (a). Panel (c): energy distribution functions f ( E ) at the final simulation time. Panels (d)–(f) show the corresponding results for the simulation including only the toroidal electric field. Panels (g)–(i) show the results of the simulation including both the synchrotron radiation losses and the toroidal electric field.

Close modal

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.

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 τ e = 2 π γ 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 MeV E 0 100 MeV , and the initial pitch angles in the range 0 ° θ 0 80 ° . For each set of E 0 and θ0 we uniformly distribute an ensemble of 4 × 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

δ B = 100 × max ( | B ( x ( t 0 + τ e ) ) B ( x ( t 0 ) ) | B ( x ( t 0 ) ) ) ,
(18)

and magnetic field direction

Δ B = 100 × max ( | B ( x ( t 0 + τ e ) ) B ( x ( t 0 ) ) | | B ( x ( t 0 ) ) | ) ,
(19)

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

Figure 13 shows the spatial distribution of δB and Δ B for DIII-D parameters with qedge = 2 for the two sets of initial conditions E 0 = 60  MeV and θ 0 = 0 ° , and E 0 = 60  MeV and θ 0 = 40 ° . For small values of θ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 δB, and at the plasma boundary for Δ B . The regions of large δB at the top and bottom of the plasma occur due to RE streaming from regions close to R = R0 to regions of lower (higher) magnetic field at R > R 0 ( R < R 0 ), while the region of large Δ 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 γ m e v 2 / ( e B 0 R 0 ) . Here v is the parallel component of the RE velocity. As the initial pitch angle θ0 increases, δB becomes larger and Δ B becomes smaller and the regions of large δB and Δ 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).

FIG. 13.

Spatial distribution of the maximum change of magnetic field amplitude δB, and direction, Δ B , along full orbit over one gyro-period for an ensemble of RE in a DIII-D-like magnetic field with qedge = 2. Panels (a) and (b) show δB and Δ B for E 0 = 60  MeV and θ 0 = 0 ° . Localized regions of large δB at the top and bottom of the plasma are observed due to electrons streaming from regions close to the magnetic axis to regions of lower (higher) magnetic field at R > R 0 ( R < R 0 ). On the other hand, large values of Δ B at the plasma boundary occur due to radial transport driven by the curvature drift. Panels (c) and (d) show δB and Δ B with E 0 = 60  MeV and θ 0 = 40 ° . The regions of large δB and Δ B become less localized as the pitch angle increases because the electrons explore more regions of the space due to their larger parallel displacement over on gyro-period.

FIG. 13.

Spatial distribution of the maximum change of magnetic field amplitude δB, and direction, Δ B , along full orbit over one gyro-period for an ensemble of RE in a DIII-D-like magnetic field with qedge = 2. Panels (a) and (b) show δB and Δ B for E 0 = 60  MeV and θ 0 = 0 ° . Localized regions of large δB at the top and bottom of the plasma are observed due to electrons streaming from regions close to the magnetic axis to regions of lower (higher) magnetic field at R > R 0 ( R < R 0 ). On the other hand, large values of Δ B at the plasma boundary occur due to radial transport driven by the curvature drift. Panels (c) and (d) show δB and Δ B with E 0 = 60  MeV and θ 0 = 40 ° . The regions of large δB and Δ B become less localized as the pitch angle increases because the electrons explore more regions of the space due to their larger parallel displacement over on gyro-period.

Close modal

A statistical summary of all the simulations performed for the qedge = 3 DIII-D type magnetic field is presented in Fig. 14 that shows the ensemble averages of the magnitude δ B , and the direction, Δ B 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 δ B , 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 Δ B increases with E 0 but decreases with increasing θ0. The values of δ B and Δ B in the DIII-D case are observed to be independent of the value of qedge, being as large as 15% and 60%, respectively. For ITER-type magnetic fields, similar trends are observed; however, the values of δ B and Δ B tend to be an order of magnitude smaller than for the DIII-D case.

FIG. 14.

Ensemble average of the maximum change of the amplitude δ B (left panel) and direction Δ B (right panel) of the magnetic field observed by RE orbits along one gyro-period τe. The values of δ B are observed to increase with both E 0 and θ0, while the values of Δ B increase with increasing E 0 and decrease with increasing θ0; in both cases, they are independent of the chosen qedge. For the ITER case, the trends are the same as for DIII-D but the values are nearly ten times smaller.

FIG. 14.

Ensemble average of the maximum change of the amplitude δ B (left panel) and direction Δ B (right panel) of the magnetic field observed by RE orbits along one gyro-period τe. The values of δ B are observed to increase with both E 0 and θ0, while the values of Δ B increase with increasing E 0 and decrease with increasing θ0; in both cases, they are independent of the chosen qedge. For the ITER case, the trends are the same as for DIII-D but the values are nearly ten times smaller.

Close modal

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 ρ L = v / Ω e , is of the order ρ L 0.02 m in (a)–(b), and of the order ρ L 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 1 m. However, the accurate evaluation of the variations of the magnetic field in one gyro-period requires considering, in addition to ρL, the parallel Larmor radius length scale ρ = v τ e = 2 π v / Ω e . For the results reported in Fig. 13, independently of the pitch angle ρ 0.6 m, which is significantly larger than ρL and comparable to LB. 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.

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 . 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, 10 11 s, and the time scales of interest, 10 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 ( θ , t = 0 ) = δ ( θ θ 0 ) / π . In the absence of collisions, radiation damping, and electric fields, reduced models trivially imply the invariance of the pitch angle distribution, i.e., f ( θ , t > 0 ) = δ ( θ θ 0 ) / π . 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.

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

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

The variable x is computed at the half time steps i 1 / 2 , i + 1 / 2 , , while p is computed at the integer time steps i , i + 1 , . In what follows, we will distinguish the momentum calculated at the time level 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

x i + 1 / 2 x i 1 / 2 Δ t = v i ,
(A1)
p L i + 1 p i Δ t = q ( E i + 1 / 2 + v i + v L i + 1 2 × B i + 1 / 2 ) ,
(A2)

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

p L i + 1 = s [ p + ( p · t ) t + p × t ] ,
(A3)
γ i + 1 = σ + σ 2 + 4 ( τ 2 + p * 2 ) 2 ,
(A4)

where we have defined

p = p i + q Δ t ( E i + 1 / 2 + v i 2 × B i + 1 / 2 ) ,
(A5)

τ = ( q Δ t / 2 ) B i + 1 / 2 , t = τ / γ i + 1 , p * = p · τ / m e c , σ = γ 2 τ 2 , γ = 1 + p 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

p R i + 1 p i Δ t = F R ( x i + 1 / 2 , p i + 1 / 2 ) ,
(A6)

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

p i + 1 = p L i + 1 + p R i + 1 p i .
(A7)
1.
M.
Lehnen
,
K.
Aleynikova
,
P. B.
Aleynikov
,
D. J.
Campbell
,
P.
Drewelow
,
N. W.
Eidietis
,
Yu.
Gasparyan
,
R. S.
Granetz
,
Y.
Gribov
,
N.
Hartmann
 et al,
J. Nucl. Mater.
463
,
39
(
2015
).
2.
E. M.
Hollmann
,
P. B.
Aleynikov
,
T.
Fülöp
,
D. A.
Humphreys
,
V. A.
Izzo
,
M.
Lehnen
,
V. E.
Lukash
,
G.
Papp
,
G.
Pautasso
,
F.
Saint-Laurent
, and
J. A.
Snipes
,
Phys. Plasmas
22
,
021802
(
2015
).
3.
T. C.
Hender
,
J. C.
Wesley
,
J.
Bialek
,
A.
Bondeson
,
A. H.
Boozer
,
R. J.
Buttery
,
A.
Garofalo
,
T. P.
Goodman
,
R. S.
Granetz
,
Y.
Gribov
 et al,
Nucl. Fusion
47
,
S128
S202
(
2007
).
4.
A. H.
Boozer
,
Phys. Plasmas
22
,
032504
(
2015
).
5.
G.
Fussmann
,
Nucl. Fusion
19
,
327
(
1979
).
6.
J. R.
Martín-Solís
,
J. D.
Alvarez
,
R.
Sánchez
, and
B.
Esposito
,
Phys. Plasmas
5
,
2370
(
1998
).
7.
M. N.
Rosenbluth
and
S. V.
Putvinski
,
Nucl. Fusion
37
,
1355
(
1997
).
8.
J. A.
Heikkinen
and
S. K.
Sipila
,
Comput. Phys. Commun.
76
,
215
(
1993
).
9.
L. G.
Eriksson
and
P.
Helander
,
Comput. Phys. Commun.
154
,
175
(
2003
).
10.
M.
Landreman
,
A.
Stahl
, and
T.
Fülöp
,
Comput. Phys. Commun.
185
,
847
(
2014
).
11.
A.
Stahl
,
E.
Hirvijoki
,
J.
Decker
,
O.
Embréus
, and
T.
Fülöp
,
Phys. Rev. Lett.
114
,
155002
(
2015
).
12.
P.
Aleynikov
and
B. N.
Breizman
,
Phys. Rev. Lett.
114
,
155001
(
2015
).
13.
E.
Nilsson
,
J.
Decker
,
Y.
Peysson
,
R. S.
Granetz
,
F.
Aint-Laurent
, and
M.
Vlainic
,
Plasma Phys. Controlled Fusion
57
,
095006
(
2015
).
14.
G.
Papp
,
M.
Drevlak
,
T.
Fülöp
, and
P.
Helander
,
Nucl. Fusion
51
,
043004
(
2011
).
15.
X.
Guan
,
H.
Qin
, and
N. J.
Fisch
,
Phys. Plasmas
17
,
092502
(
2010
).
16.
A. J.
Russo
,
Nucl. Fusion
31
,
117
(
1991
).
17.
L.
Jian
,
W.
Yulei
, and
Q.
Hong
,
Nucl. Fusion
56
,
064002
(
2016
).
18.
W.
Yulei
,
Q.
Hong
, and
L.
Jian
,
Phys. Plasmas
23
,
062505
(
2016
).
19.
H.
Spohn
,
Europhys. Lett.
50
,
287
(
2000
).
20.
F.
Rohrlich
,
Phys. Lett. A
283
,
276
(
2001
).
21.
L.
Landau
and
E.
Lifshitz
,
The Classical Theory of Fields
, 3rd revised English ed., translated by M. Hamermesh (
Pergamon Press
,
New York
,
1971
).
22.
J.-L.
Vay
,
Phys. Plasmas
15
,
056701
(
2008
).
23.
M.
Tamburini
,
F.
Pegoraro
,
A.
Di Piazza
,
C. H.
Keitel
, and
A.
Macchi
,
New J. Phys.
12
,
123005
(
2010
).
24.
R.
Zhang
,
J.
Liu
,
H.
Qin
,
Y.
Wang
,
Y.
He
, and
Y.
Sun
,
Phys. Plasmas
22
,
044501
(
2015
).
25.
E. M.
Hollmann
,
M. E.
Austin
,
J. A.
Boedo
,
N. H.
Brooks
,
N.
Commaux
,
N. W.
Eidietis
,
D. A.
Humphreys
,
V. A.
Izzo
,
A. N.
James
,
T. C.
Jernigan
 et al,
Nucl. Fusion
53
,
083004
(
2013
).
26.
H.
Knoepfel
and
D. A.
Spong
,
Nucl. Fusion
19
,
785
(
1979
).
27.
J.
Schwinger
,
Phys. Rev.
75
,
1912
(
1949
).
28.
A.
Stahl
,
M.
Landreman
,
G.
Papp
,
E.
Hollmann
, and
T.
Fülöp
,
Phys. Plasmas
20
,
093302
(
2013
).
29.
J. H.
Yu
,
E. M.
Hollmann
,
N.
Commaux
,
N. W.
Eidietis
,
D. A.
Humphreys
,
A. N.
James
,
T. C.
Jernigan
, and
R. A.
Moyer
,
Phys. Plasmas
20
,
042113
(
2013
).
30.
T. G.
Northrop
,
The Adiabatic Motion of Charged Particles
(
Interscience
,
New York
,
1963
).
31.
R. G.
Littlejohn
,
Phys. Fluids
24
,
1730
(
1981
).