Calculation of collisionless pitch-angle scattering of runaway electrons with synchrotron radiation via high-order guiding-centre equation

Abstract

Recently, the collisionless pitch-angle scattering for relativistic runaway electrons (REs) in toroidal geometries such as tokamaks was discovered through a full orbit simulation approach (Liu et al., Nucl. Fusion, vol. 56, 2016, p. 064002), and it was then theoretically investigated that a new expression for the magnetic moment, including the second-order corrections, could essentially reproduce the so-called collisionless pitch-angle scattering process (Liu et al., Nucl. Fusion, vol. 58, 2018, p. 106018). In this paper, with synchrotron radiation, extensive numerical verification of the validity of the high-order guiding-centre theory is given for simulations involving REs by incorporating such an expression for the magnetic moment into our particle tracing code. A high-order guiding-centre simulation approach with synchrotron radiation (HGSA) is applied. Synchrotron radiation plays an essential role in the life cycle of REs. The energy of REs first increases and then becomes saturated until the electric field acceleration is balanced by the radiation dissipation. Unfortunately, the process cannot be simulated accurately with the standard guiding-centre model, i.e. the first-order guiding-centre model. Remarkably, it is found that the HGSA can effectively produce the fundamental process of REs. Since the time scale of the energy saturation of REs is close to seconds, the computational cost becomes significant. In order to save costs, it is necessary to estimate the time of energy saturation. An analytical estimate is derived for the time it takes for synchrotron drag to balance an accelerating electric field and the provided formula has been numerically verified. Test calculations reveal that HGSA is favourable for exploiting the dynamics of REs in tokamak plasmas.

1 Introduction

Plasma disruption (Lehnen, Aleynikova & Aleynikov 2015; Tang, Zeng & Chen 2021) is almost inevitable in the operation of the tokamak. One of the most dangerous hazards of plasma disruption is the generation of large amounts of runaway electrons (REs) (Paz-Soldan, Cooper & Aleynikov 2017; Breizman et al. 2019; Tinguely, Izzo & Garnier 2021). The underlying reason is that the friction force decreases with increasing speed for fast electrons (Dreicer 1959). The electrons can then be accelerated by an electric field without limitation (Aleynikov & Breizman 2015; Stahl et al. 2015), which is now known as the runaway phenomenon. Thus, it is possible that a small population of REs becomes a large current carrier (Chen, Huang & Luo 2016; Cai et al. 2021). In recent years, both experiments and simulations have illustrated that REs have an impact on plasma instability (Cai & Fu 2015; Liu, Hirvijoki & Fu 2018a; Spong, Heidbrink & Paz-Soldan 2018; Liu, Zhao & Jardin 2021). The energy of the RE is close to several tens to hundreds of MeV. Attention should be paid to the following points: the relativistic effect is significant; the radiation reaction (Finken, Watkins & Rusbüldt 2011; del Castillo-Negrete et al. 2018; Hoppe et al. 2018; Zhang, Zhou & Hu 2018,  2021) is important; when the inductive electric field becomes strong enough, the runaway avalanche (Nilsson et al. 2015; Svensson et al. 2021) develops after the thermal quench. Without proper control, these energetic REs will eventually hit the first wall of the device in a localized manner, posing a serious threat to the safe operation of the device.

In a recent investigation (Liu, Wang & Qin 2016), collisionless pitch-angle scattering was found to be much stronger than collisional scattering and invalidates the first-order guiding-centre model for REs. Taking synchrotron radiation, the loop electric field and the magnetic field into account, the phase space of REs has been studied by using the full orbit model. The energy increases until the synchrotron drag balances the electric field. When the parallel momentum becomes dozens of times $m_0c$ ($m_0$ is the rest mass of an electron, and c is the speed of light), the absolute magnitude of the perpendicular momentum changes a lot during one gyro-period, i.e. fast oscillations. The term ‘fast’ relates to the large change of the perpendicular momentum in one gyro-period of approximately $10^{-12}\sim 10^{-10}\,\mathrm {s}$, compared with the bounce period of REs, which is close to $10^{-7}\sim 10^{-6}\,\mathrm {s}$. With fast oscillations, the perpendicular momentum grows with increasing parallel momentum. This phenomenon is called ‘collisionless pitch-angle scattering’. It is a way in which the parallel momentum is converted to perpendicular momentum in the absence of collisions. However, for the long-term simulation of REs, the computing costs of the full orbit model are nearly unacceptable.

According to the results in Liu, Qin & Hirvijoki (2018b), a new expression for the magnetic moment is introduced, where the magnetic moment $\mu$ can still be considered as a good adiabatic invariant. The magnetic moment $\mu$ is decided not only by the perpendicular momentum but also by the parallel momentum and the magnetic-field-line curvature vector $\boldsymbol {\kappa }$, thus providing the aforementioned collisionless momentum transfer. Based on the adiabatic invariant, the Hamiltonian high-order equations of motion were derived. However, for more realistic simulation of REs, the radiation term needs to be considered.

In this paper, focusing on synchrotron radiation, we investigate the dynamics of relativistic REs by utilizing different models in the Chinese Fusion Engineering Testing Reactor (CFETR) steady-state scenario via a particle tracing code (PTC) (Wang et al. 2021), which was developed recently. The code includes both the full orbit and guiding-centre orbit models to track charged particles in a tokamak with a three-dimensional field. In PTC, the computational domain is divided into two regions: the core region and the scrape-off-layer region. The meshes are generated by a constrained Delaunay triangulation method. On each flux surface, the vertices are placed at equidistant arc lengths in the poloidal section. With the unstructured triangle mesh, the geometry of the tokamak poloidal section can be described precisely, and the flux coordinate can be interpolated conveniently. From the g-eqdsk file, bicubic interpolation on a uniform grid is used to obtain the first derivative of the magnetic field of the vertices of the triangular meshes. Bilinear finite element methods are used to calculate the second derivative of the magnetic field at any point. It now also includes a relativistic (Wimmel 1983) and high-order modification. Extensive numerical verifications of the validity of the high-order guiding-centre theory are given for simulations involving REs. The simulation results prove that the high-order guiding-centre approach (HGSA) is more appropriate than the first-order guiding-centre approach for REs. HGSA reproduces the ‘collisionless pitch-angle scattering’, but the first-order approach cannot. The perpendicular momentum is the dominant factor of the synchrotron radiation. Synchrotron radiation plays an essential role in the life cycle of REs. The energy of REs first increases and then becomes saturated until the electric field acceleration is balanced by the radiation dissipation. HGSA has reproduced the basic physical process of energy saturation. Since the time scale of the energy saturation of REs is close to seconds, the computational cost becomes significant. Avoiding unnecessary waste of calculation costs, an analytical formula is derived to calculate the balance time. In certain test cases, the provided formula has been numerically verified. Test calculations reveal that HGSA is admittedly favourable for exploiting the dynamics of REs in tokamak plasmas.

In § 2.1, we use the relativistic full orbit model to advance the test particle in the six-dimensional phase space. We give some physical explanations for the details of the development of the perpendicular momentum. In § 2.2, we carry out the work by using the high-order guiding-centre model, calculate the radiation effect and derive the analytical formula to estimate the energy balance time. In § 2.3, we analyse the relationship between the conservation of the magnetic moment and the change of magnetic field in one gyro-period and then discuss the applicability of the high-order guiding-centre equations. Finally, in § 3, we conclude our work.

2 Methods and results

2.1 Relativistic full orbit model

To analyse the phase-space characteristics of energetic particles, the test particle method is an effective numerical method. Based on this method, many important physical results have been produced (Sommariva et al. 2017; Hao, Chen & Cai 2020; Zhao, Wang & Wang 2020). The method advances the particle motion through time-dependent orbital equations, while the electromagnetic fields are not solved self-consistently. Two kinds of particle motion equation are used currently, one is the standard guiding-centre equations, i.e. the first-order guiding-centre equations, and the other is to use the Lorentz orbital equations, i.e. the full orbit equations. Considering computational efficiency, the guiding-centre model could be made more applicable by averaging out the fast gyro-motion. However, the first-order guiding-centre model may not be suitable for the simulation of REs (Liu et al. 2016).

For the relativistic full orbit calculation, the equations of motion can be written as

\begin{equation} \left. \begin{gathered} \frac{\rm{d} \boldsymbol{x}}{\rm{d} t}=\boldsymbol{v}, \\ \frac{\rm{d} \boldsymbol{p}}{\rm{d} t}=q\left(\boldsymbol{E}+ \boldsymbol{v}\times\boldsymbol{B} \right), \end{gathered} \right\} \end{equation}

where $\boldsymbol {x},\ \boldsymbol {v},\ \boldsymbol {p}=\gamma m_0 \boldsymbol {v}$ represent the physical location, velocity and mechanical momentum of the particle respectively, $\boldsymbol {E}$ and $\boldsymbol {B}$ denote the electric and magnetic fields and the Lorentz factor $\gamma$ is given as

(2.2)\begin{equation} \gamma=\sqrt{1+\frac{p^2}{m_0^2c^2}}=\frac{1}{\sqrt{1-(v/c)^2}}. \end{equation}

Here, $c$ is the speed of light in a vacuum, v is the velocity of the particle. The toroidal electric field is given as $\boldsymbol {E}=E_l R_0/R\boldsymbol {e}_\phi$. The major radius of CFETR is $R_0=7.2\,\mathrm {m}$, the central magnetic field is $B_0=6.5\,\mathrm {T}$ and $E_l$ is given as $E_l=10.0\,\mathrm {V}\,\mathrm {m}^{-1}$. Now we include a test energetic electron. The initial physical location is $R=8.0\,\mathrm {m},\ \phi = 0, Z=0\,\mathrm {m}$, and the initial momentum is $p_{\parallel 0}=5m_0 c, p_{\bot 0}=m_0 c$. The time step of the simulation is set to be $\delta T =0.01T_{ce}$, $1\,\%$ of the gyro-period $T_{ce}$.

Figure 1 shows the phase space of the RE. In a real tokamak, the evolution trend is consistent with the results in an ideal magnetic geometry in Liu et al. (2016). Without radiation, figure 1(a) depicts the increase of the parallel momentum due to the electric field. Figure 1(b) demonstrates the evolution trend of the perpendicular momentum. Zooming in, we see that the time scale is approximately $10^{-8}$ s, and it is clear that the interval between the two maximum values of $p_\bot$ is approximately $10^{-9}$ s. Hence, the oscillation of $p_\bot$ can be regarded as a fast oscillation. The amplitude increases until the minimum of the perpendicular momentum reaches a zero point. As time goes by, the perpendicular momentum becomes greater than zero and oscillates at an almost fixed amplitude. When considering the magnetic fluctuations or the ripple fields in a real tokamak, this changes the radiation effects of REs (Martín-Solís, Sánchez & Esposito 1999; Liu et al. 2016). Figure 1(c) depicts that the lowest magnetic moment $\mu _0$ increases and oscillates, which is against the assumptions of the first-order guiding-centre model. The lowest-order magnetic moment is no longer conservative for REs.

Figure 1. Evolution of (a) the parallel momentum, (b) the perpendicular momentum and (c) the lowest magnetic moment of the electron.

The asymptotic expansion is $\mu =\mu _0 + \epsilon \mu _1 + \cdots$, where $\mu _1$ represents the first-order correction that explicitly involves the non-uniformity of the magnetic field (Burby, Squire & Qin 2013). In Liu et al. (2018b), keeping to the second order, a new expression was derived under the limit $p_\parallel \gg p_\bot$. The expression consists of three terms

(2.3)\begin{gather} \mu=\frac{|\boldsymbol{p}_\bot+p_\parallel^2 \boldsymbol{\kappa}\times \boldsymbol{b} /\left(q B\right)|^2}{\left(2 m_0 B\right)}=\mu_0+\mu_1+\mu_2, \end{gather}

(2.4a–c)\begin{gather}\mu_0=\frac{p_\bot^2}{2 m_0 B},\quad \mu_1=\frac{p_\parallel^2\boldsymbol{p}_\bot\boldsymbol{\cdot}\boldsymbol{\kappa}\times\boldsymbol{b}}{q m_0 B^2},\quad \mu_2=\frac{p_\parallel^4}{q^2 B^2}\frac{\left|\boldsymbol{\kappa}\times\boldsymbol{b}\right|^2}{2 m_0 B}, \end{gather}

where $\boldsymbol {\kappa }$ is the curvature vector defined by $\boldsymbol {\kappa }=\boldsymbol {b}\boldsymbol {\cdot }\boldsymbol {\nabla }\boldsymbol {b}$, and the unit vector $\boldsymbol {b}=\boldsymbol {B}/B$. Furthermore, an approximate value $\widetilde {p_\bot }$ is given by solving (2.3) and ignoring the cross-product,

(2.5)\begin{equation} \widetilde{p_\bot}=\sqrt{\left(p_\parallel^2\boldsymbol{\kappa}\times \frac{\boldsymbol{b}}{q B}\right)^2+2\mu m_0 B}. \end{equation}

Figure 2(a) shows that $\mu _0$ is no longer conservative, but $\mu$ still holds as an adiabatic invariant; figure 2(b) depicts the advancement of $p_\bot ^2$ and $\widetilde {p_\bot }^2$, which demonstrates that the evolution trend is consistent with ‘collisionless pitch-angle scattering’.

Figure 2. Evolution of (a) the lowest magnetic moment $\mu _0$ and the new magnetic moment $\mu$, and (b) $p_\bot ^2$ and $\widetilde {p_\bot }^2$.

In (2.3), the fraction includes the vectors $\boldsymbol {p}_\bot$ and $\boldsymbol {p}_C=-p_\parallel ^2\boldsymbol {\kappa }\times \boldsymbol {b}/q B$. The latter is the curvature drift. To better understand the collisionless pitch-angle scattering, we can consider an ideal and simple physical model. We can treat the complete perpendicular momentum vector $\boldsymbol {p}_\bot$ as the sum of $\boldsymbol {p}_{\bot 0}$ and $\boldsymbol {p}_C$, where $\boldsymbol {p}_{\bot 0}$ spans a symmetric circle around $\boldsymbol {b}$. (While the value of the parallel momentum is not so large, $\boldsymbol {p}_{\bot 0}$ is almost the initial perpendicular momentum.) The conservative magnetic moment becomes $\mu =|\boldsymbol {p}_\bot -\boldsymbol {p}_C|^2/(2 m_0 B)=|\boldsymbol {p}_{\bot 0}|^2/(2 m_0 B)$. With the increase of the parallel momentum, $\boldsymbol {p}_C$ grows and becomes comparable to $\boldsymbol {p}_{\bot 0}$, so we cannot merely regard $\boldsymbol {p}_C$ as the result of drift precession; it greatly changes the perpendicular momentum all the time and breaks the conservation of the lowest magnetic moment.

In CFETR, $\boldsymbol {\kappa }\times \boldsymbol {b}$ mainly tilts towards the $Z$-direction. The perpendicular momentum is almost in the $R$–$Z$ plane. In the beginning, the minimum is decided by $|\boldsymbol {p}_{\bot }|_{\min }=|(\boldsymbol {p}_{\bot 0}+\boldsymbol {p}_C)|_{\min }=|\boldsymbol {p}_{\bot 0}|-|\boldsymbol {p}_C|$ and the maximum is $|\boldsymbol {p}_{\bot }|_{\max }=|(\boldsymbol {p}_{\bot 0}+\boldsymbol {p}_C)|_{\max }=|\boldsymbol {p}_{\bot 0}|+|\boldsymbol {p}_C|.$ The minimum and the maximum of the perpendicular momentum also mainly tilt towards the $Z$-direction. The amplitude $A_{mp}$ is decided by $A_{mp}=|\boldsymbol {p}_{\bot }|_{\max }-|\boldsymbol {p}_{\bot }|_{\min }=2|\boldsymbol {p}_C|$. With the increase of the parallel momentum, $|\boldsymbol {p}_{\bot }|_{\min }$ decreases, $|\boldsymbol {p}_{\bot }|_{\max }$ and $A_{mp}$ continues to grow until $|\boldsymbol {p}_{\bot }|_{\min }$ reaches zero and $|\boldsymbol {p}_C|$ equals $|\boldsymbol {p}_{\bot 0}|$. As time goes by, $|\boldsymbol {p}_{\bot }|_{\min }$ increases and becomes greater than zero. $|\boldsymbol {p}_{\bot }|_{\min }$ becomes $|\boldsymbol {p}_C|-|\boldsymbol {p}_{\bot 0}|$ and then $|\boldsymbol {p}_{\bot }|_{\min }$ and $|\boldsymbol {p}_{\bot }|_{\max }$ increase and the amplitude $A_{mp}$ becomes almost a constant, i.e. $A_{mp}=2|\boldsymbol {p}_{\bot 0}|$. In Wang, Qin & Liu (2016), when the pitch angle is in the range $[-0.5,0.5]$, $A_{mp}$ is almost a constant, this is because $A_{mp}$ is nearly proportional to $p_{\perp 0}$, which changes a little,varying in the range around $[0.87,1]$ $p_0$. In figure 3 we reveal the diagram for the rotations of the momentum when the parallel momentum of the particle is $200 m_0 c$. The red solid line is the motion of the cyclotron centre of the electron, and it moves from point A to point B. The coloured line denotes the endpoint of the vector of the momentum, which starts from the red line. The distance of the corresponding points between the two lines represents the relative value of the momentum. Figure 3(a,b) plots $\boldsymbol {p}_\bot$ and panels (c,d) plot $(\boldsymbol {p}_\bot -\boldsymbol {p}_C)$. Clearly, the amplitude is approximately $A_{mp}= 2p_{\bot 0}=2 m_0 c$. In figure 3(c,d), the track of the vector $(\boldsymbol {p}_\bot -\boldsymbol {p}_C)$ spans a symmetric circle around $\boldsymbol {b}$, which is the same as the rotation of $\boldsymbol {p}_{\bot 0}$. In figure 4(a), outputting data every $10^9$ steps shows that the $Z$-component of the perpendicular momentum keeps growing, which is the result of the increase of $p_\parallel$. In figure 4(b), comparing the absolute value with the $Z$-component value of the perpendicular momentum, we demonstrate that the ranges of the perpendicular momentum and the $Z$-component are almost the same, which means the extremum appears while $\boldsymbol {p}_\bot$ mainly slants towards the $Z$-direction.

Figure 3. Diagram for the rotations of the momentum vector. The red solid line is the motion of the cyclotron centre, and it moves from point A to point B. The coloured line denotes the endpoint of the vector, which starts from the red line. The distance of the corresponding points between the two lines represents the relative value of the momentum.

Figure 4. Evolution of (a) the $Z$-component of the perpendicular momentum, and (b) the perpendicular momentum and the absolute value of the $Z$-component of the perpendicular momentum.

2.2 High-order guiding-centre orbit model with synchrotron radiation

Synchrotron radiation plays a key part in REs during their entire life cycles. More than that, in the generation of REs, experiments (Martín-Solís, Sánchez & Esposito 2010) show that the synchrotron radiation increases the threshold electric field. For a relativistic particle, the radiation-reaction force (RR-force) (Pauli 1958) is

(2.6)\begin{equation} \boldsymbol{K} = \frac{q^2\gamma^2}{6{\rm \pi} \epsilon_0 c^3}\left[\ddot{\boldsymbol{v}}+ \frac{3\gamma^2}{c^2}(\boldsymbol{v}\boldsymbol{\cdot}\dot{\boldsymbol{v}})\dot{\boldsymbol{v}}+ \frac{\gamma^2}{c^2}\left(\boldsymbol{v}\boldsymbol{\cdot}\ddot{\boldsymbol{v}}+ \frac{3\gamma^2}{c^2}(\boldsymbol{v}\boldsymbol{\cdot}\dot{\boldsymbol{v}})^2\right)\boldsymbol{v}\right], \end{equation}

and for the relativistic full equation, the RR-force can be simplified to (Hirvijoki et al. 2015)

(2.7)\begin{equation} \boldsymbol{K}={-}\nu_r\left(\boldsymbol{p}_\bot+ \frac{p_\bot^2}{(m_0 c)^2}\boldsymbol{p}\right)-\nu_r \frac{{\dot{\boldsymbol{B}}\times\boldsymbol{p}}}{\varOmega B},\end{equation}

where the Larmor frequency is $\varOmega =eB/\gamma m_0$. The characteristic time for the RR-force is

(2.8)\begin{equation} \nu_r^{{-}1}=\frac{6{\rm \pi} \epsilon_0 \gamma (m_0 c)^3}{q^4 B^2}, \end{equation}

where $\epsilon _0$ is the permittivity of a vacuum. After a series of algebraic treatments of the formula, the components for the guiding-centre radiation reaction force in magnetic-field non-uniformity are derived (Hirvijoki et al. 2015).

Based on the non-canonical Hamiltonian mechanics approach using the Lie perturbation method, the high-order relativistic guiding-centre equations (Liu et al. 2018b) are derived. Besides, in Appendix A, by separating the perpendicular momentum, the equations are derived from some other physical perspectives. The high-order relativistic guiding-centre equations with synchrotron radiation are

(2.9)\begin{equation} \left.\begin{gathered} \dot{\boldsymbol{X}}=\frac{p_\parallel^\ast}{\gamma_r m_0} \frac{\boldsymbol{B}^\ast}{B_\parallel^\ast}+\frac{\boldsymbol{b}^\ast}{e B_\parallel^\ast}\times \left(\boldsymbol{\nabla} H+\frac{\partial\boldsymbol{A}^\ast}{\partial t}\right)+K^{\boldsymbol{X}}, \\ {\dot{p}}_\parallel={-}\frac{\boldsymbol{B}^\ast}{B_\parallel^\ast}\boldsymbol{\cdot} \left(\boldsymbol{\nabla} H+\frac{\partial\boldsymbol{A}^\ast}{\partial t}\right)+K^{p_\parallel},\\ \dot{\mu}=K^\mu, \end{gathered}\right\} \end{equation}

where

(2.10)\begin{equation} \left.\begin{gathered} \boldsymbol{A}^\ast=\boldsymbol{A}+\frac{p_\parallel}{q}\boldsymbol{b}-p_\parallel^2 \boldsymbol{\kappa}\times\frac{\boldsymbol{b}}{q^2 B}, \\ \boldsymbol{B}^\ast=\boldsymbol{\nabla}\times\boldsymbol{A}^\ast, \\ \boldsymbol{b}^\ast=\boldsymbol{b}-2 p_\parallel\boldsymbol{\kappa}\times \frac{\boldsymbol{b}}{q B}, \\ B_\parallel^\ast=\boldsymbol{B}^\ast\boldsymbol{\cdot}\boldsymbol{b}^\ast,\\ \boldsymbol{\nabla} H=\frac{\mu}{\gamma_r}\boldsymbol{\nabla} B+\left(\frac{1}{\gamma_r} \frac{p_\parallel^4}{2 m_0 \textit{q}^2}\boldsymbol{\nabla}\frac{\kappa^2}{B^2}\right)+q\boldsymbol{\nabla}\varPhi,\\ p_\parallel^\ast=p_\parallel+{2p}_\parallel^3\frac{\kappa^2}{{q^2B}^2},\\ \gamma_r=\sqrt{1+\frac{p_\parallel^2\left[1+p_\parallel^2\kappa^2/\!\left(q B\right)^2+2\mu m_0 B\right]}{(m_0 c)^2}}, \end{gathered}\right\} \end{equation}

and the components for the guiding-centre RR-force are

(2.11)\begin{equation} \left.\begin{gathered} K^{\boldsymbol{X}}={-}\frac{\nu_r}{\varOmega_\parallel^\ast} \frac{p_\bot^2 }{(m_0 c)^2}(\boldsymbol{b}\times\dot{\boldsymbol{X}}+3v_\parallel \rho_\parallel \boldsymbol{\kappa}), \\ K^{p_\parallel}={-}\nu_r p_\parallel \frac{p_\bot^2}{2(m_0 c)^2}(2+\rho_\parallel \tau_B)-\nu_r \frac{p_\bot \gamma^2}{2}\rho_\bot\tau_B, \\ K^\mu={-}\nu_r\mu\left(1+\frac{p_\bot^2}{(m_0 c)^2}\right)(2+\rho_\parallel \tau_B). \end{gathered}\right\} \end{equation}

The magnetic field-line twist parameter $\tau _B$ is $\tau _B=\boldsymbol {b}\boldsymbol {\cdot }\boldsymbol {\nabla }\times \boldsymbol {b}$. The parallel and perpendicular gyro-radii and the modified gyro-frequency are separately defined as

(2.12a–c)\begin{equation} \rho_\parallel=\frac{p_\parallel}{q B},\quad \rho_\bot=\frac{p_\bot}{q B},\quad \varOmega_\parallel^\ast=\frac{q B_\parallel^\ast}{m_0\gamma}.\end{equation}

We perform the simulation of the test particle under three models without radiation. The initial values are the same as section § 2.1. The energy of the electron in the snapshot of figure 5 is approximately 70 MeV after the acceleration. This signifies the comparison of the orbital trajectory of an electron with different models in one bounce period. The high-order guiding-centre model is consistent with the full orbit. Figure 6 depicts the comparison of the momentum evolution. It demonstrates that HGSA reproduces the ‘collisionless pitch-angle scattering’ phenomenon. The results confirm the applicability of the high-order guiding-centre model in a real magnetic configuration of the tokamak. Note that the evolution of the perpendicular momentum has a great influence on the process of the energy balance because of the radiation, and the modification of the perpendicular momentum is significant. Thus, for the simulation of REs, the high-order guiding-centre model is more appropriate.

Figure 5. Comparison of the orbital trajectory of an electron. The light blue solid line plots the relativistic full orbit, the dark blue dashed line plots the high-order guiding-centre orbit and the red solid line plots the first-order guiding-centre orbit.

Figure 6. Comparison of the evolution of the momentum of an electron. The light blue line plots the relativistic full orbit, the dark blue circle (hollow) points plot the high-order guiding-centre orbit and the red solid points plot the first-order guiding-centre orbit.

The perpendicular momentum is dominant for the radiation effect at a finite pitch angle (Martín-Solís et al. 1998). To see a clearer radiation effect, we perform a test particle simulation in which the initial momentum is set to be $p_{\bot 0}=1.7m_0 c$ and $p_{\parallel }=200m_0 c$ in the high-order guiding-centremodel, while $|\boldsymbol {p}_{\bot }|=|\boldsymbol {p}_{\bot 0}+\boldsymbol {p}_C|=2.4m_0 c$ and $p_{\parallel }=200m_0 c$ in the full orbit model. The parallel radiation could be considered as an effective electric field. Figure 7 indicates the comparison of the advancement of (a) the perpendicular momentum, (b) the equivalent electric field of parallel radiation and (c) the kinetic energy loss of an electron in the equilibrium magnetic field of the CFETR without the electric field, where the kinetic energy is $T=(\gamma -1)m_0c^2$. In figure 7(b), in detail, the real radiation of the electron motion has zero points, but the high-order guiding-centre model aims to decouple the motion, such that it cannot accurately simulate the details. Figure 7(c) demonstrates the kinetic energy loss separately through the full orbit model and the high-order guiding-centre model. It is clear that, under the two models, the loss is approximately equal.

Figure 7. Comparison of the evolution of (a) the perpendicular momentum, (b) the equivalent electric field of parallel radiation, (c) the kinetic energy loss of an electron in the equilibrium magnetic field of the CFETR without an electric field. The initial momentum is $p_{\bot 0}=1.7m_0 c,\ p_{\bot }=2.4m_0 c$ and $p_{\parallel }=200m_0 c$.

Now we carry out the simulation work including the radiation effect and the accelerating electric field. The parameter of the loop electric field is $E_l=10.0\,\textrm {V}\,\textrm {m}^{-1}$. For a RE in CFETR, the initial momentum of the RE is $p_{\parallel 0} = 500m_0 c$ and $p_{\bot 0} = 7.5m_0 c$. The initial physical location is $R=8.0\,\mathrm {m},\ \phi = 0,\ Z=0\,\mathrm {m}$. Figure 8 reveals the evolution of (a) the parallel momentum, (b) the equivalent electric field of parallel radiation and (c) the kinetic energy of an electron in the equilibrium magnetic field of the CFETR and the loop electric field. It shows that the basic physical process of the energy balance is consistent under the two models.

Figure 8. Comparison of the evolution of (a) the parallel momentum, (b) the equivalent electric field of parallel radiation, (c) the kinetic energy of an electron in the equilibrium magnetic field of the CFETR with an electric field. The initial momentum of the RE is $p_{\parallel 0} = 500m_0 c$ and $p_{\bot 0} = 7.5m_0 c$.

In the simulation work, generally, we need to set up a simulation stop time first, and this greatly affects the calculation cost. It is meaningful if the energy balance time of the electron can be calculated. For REs, we set $p_\bot /p_\parallel =\epsilon$. We keep to the zero order in the following calculation work.

In (2.9) and (2.11), only considering energy gain and radiation dissipation, we can obtain

(2.13)\begin{equation} {\dot{p}}_\parallel= E_\parallel q+K^{p_\parallel}.\end{equation}

For the radiation dissipation term, we divide it into three parts,

(2.14)\begin{equation} \left.\begin{gathered} K^{p_\parallel}=K^{p_\parallel}_a+K^{p_\parallel}_b+K^{p_\parallel}_c,\\ K^{p_\parallel}_a ={-}\nu_r p_\parallel \frac{p_\bot^2}{(m_0 c)^2},\quad K^{p_\parallel}_b ={-}\nu_r p_\parallel \frac{p_\bot^2}{2(m_0 c)^2}\rho_\parallel \tau_B,\quad K^{p_\parallel}_c ={-}\nu_r\frac{p_\bot \gamma^2}{2}\rho_\bot\tau_B. \end{gathered}\right\} \end{equation}

First of all, we have

(2.15)\begin{equation} \left.\begin{gathered} \frac{K^{p_\parallel}_b}{K^{p_\parallel}_a} = \frac{\rho_\parallel \tau_B}{2} = \frac{\rho_\parallel}{2L} \sim O(\epsilon), \\ \frac{K^{p_\parallel}_c}{K^{p_\parallel}_b} = \frac{p_\perp \rho_\perp \gamma^2 }{p_\parallel \rho_\parallel p_\perp^2/{(m_0c)^2} } = \left(\frac{\gamma}{p_\parallel{/}{m_0c}}\right)^2\sim O(1), \end{gathered}\right\} \end{equation}

where $L$ is the magnetic-field non-uniformity length scale. When the radiation effect becomes dominant, the parallel momentum is far more than $m_0c$, giving the approximation $\gamma =\sqrt {1+(p/(m_0c))^2}\approx p_\parallel /(m_0c)$. We keep to zero order, i.e. keeping $part\ a$ of the radiation term. The approximate perpendicular momentum is given in (2.5). The equation (2.13) turns out to be

(2.16)\begin{equation} {\dot{p}}_\parallel= E_\parallel q -\frac{q^4B^2}{6{\rm \pi} \epsilon_0(m_0c)^4} \left(\frac{p_\parallel^4 \kappa^2}{q^2B^2}+2\mu m_0B \right).\end{equation}

Assuming that $E = E_l,B = B_0$ are constants, we can obtain the analytic solution in mathematics

(2.17)\begin{equation} \left.\begin{gathered}\frac{{\rm d} (p_\parallel{/}(m_0c))}{\rm{d} t} =\frac{\rm{d} p_\parallel}{\rm{d} t} = c_1 -c_2{p_\parallel}^4, \\ c_1 =\frac{E_\parallel q}{m_0c}-\frac{q^4B^2}{6{\rm \pi} \epsilon_0(m_0c)^5}2\mu m_0B,\quad c_2 =\frac{q^4B^2}{6{\rm \pi} \epsilon_0m_0c}\frac{\kappa^2}{q^2B^2}. \end{gathered}\right\} \end{equation}

The solution of (2.17) is

(2.18)\begin{equation} \left(\frac{1}{2}\ln{\frac{\sqrt{c_1}-\sqrt{c_2}p_\parallel^2}{(c_1^{{1}/{4}}-c_2^{{1}/{4}}p_\parallel)^2}}+ \arctan{\left(\left(\frac{c_2}{c_1}\right)^{{1}/{4}}p_\parallel\right)}\right) \frac{1}{2\sqrt{c_1}(c_1c_2)^{{1}/{4}}}=t-t_0,\end{equation}

where $p_\parallel$ is a number that is a multiple of $m_0c$, and $t_0$ is the integration constant, with $t = 0,\ p_\parallel = p_{\parallel {0}}$. The parameters $c_1$ and $c_2$ include information on the electric and magnetic fields, the geometric parameters of the device and so on. According to the asymptotic theory, taking $t \rightarrow +\infty$, we can figure out that the limit of $p_{\parallel \max }$ is approximately $p_\parallel = c_1^{{1}/{4}}/c_2^{{1}/{4}}$. The energy saturation time can be solved by giving $p_\parallel = f p_{\parallel \max }$, where $f$ is an estimation factor based on computation needs, say $f=0.9$ or $f=0.99$. In figure 9, we show the comparison of the analytical solutions and simulation results under two electric fields, and they are consistent with each other. Hence, (2.18) can be used to calculate the energy saturation time and other related physical quantities.

Figure 9. The comparison of the momentum (kinetic energy) balance process between the simulation and analysis. The different colours show the electric fields. The points are the simulation results, and the lines represent the analytical solutions.

2.3 Applicability

This section separately focuses on the greatest variation of the following three quantities: the magnetic field $\delta \boldsymbol {B}=|(\boldsymbol {B}_{t=t_0}-\boldsymbol {B}_{t=0})|/|\boldsymbol {B}_{t=0}|$, the maximum amplitude of the lowest magnetic moment $\delta \mu _0 =(\mu _{0t=t_0}-\mu _{0t=0})/\mu _{0t=0}$, i.e. the conservation of $\mu _0$, and the maximum amplitude of the magnetic moment $\delta \mu =(\mu _{t=t_0}-\mu _{t=0})/\mu _{t=0}$ in one gyro-period. From the above, we know that the collisionless pitch-angle scattering is greatly influenced by the curvature of the device. For a typical configuration for the tokamak field, the smaller the major radius $R_0$ is, the bigger the curvature is, and thus the more obvious the phenomenon is. In this section, we carry out the simulation in Experimental Advanced Superconducting Tokamak, in which the major radius is $R_0= 1.85\,\mathrm {m}$. The three sampling kinetic energy values are chosen as $T= 168.4\,\mathrm {MeV}(\gamma = 400)$, $T= 126.3\,\mathrm {MeV}(\gamma = 300)$ and $T= 41.8\,\mathrm {MeV}(\gamma = 100)$. The pitch angle $p_\parallel /p_0$ is in the range $[0,1]$. In figure 10, $\delta \boldsymbol {B}$ and $\delta \mu _0$ increase as pitch angle grows. This means that, when the parallel motion scale is close to the perpendicular motion scale (pitch $\approx 0$ and $\approx 1$), the change in the magnetic field caused by the former is more pronounced. The value of $\delta \mu$ increases with energy, but has little to do with pitch angles. For the highly energetic particles ($\gamma = 400$ and 300), the $\delta \boldsymbol {B}$ values are all above $30\,\%$, and in our perception, this breaks the condition of conservation of the lowest magnetic moment. However, the results turn out to not quite fit with the theoretical analysis in figure 11. For the two points 0 and 1 in which the $\delta \boldsymbol {B}$ values are both almost $40\,\%$, the value of $\delta \mu _0$ is just a few per cent when the pitch approaches zero (point 0), i.e. $\mu _0$ is still an adiabatic invariant; the value of $\delta \mu _0$ reaches $20\,\%$ when the pitch is around $0.5$ (point 1), breaking the conservation of $\mu _0$. In fact, although the perpendicular motion has a great influence on the change of magnetic field (point 0), it has little impact on the conservation of the lowest magnetic moment. From the above, it is clear that the curvature drift leads to the fact that the parallel motion has an obvious effect on $\mu _0$. To sum up, in a small size tokamak, for REs, in the general energy range (probably $100\,\mathrm {MeV}$), it is appropriate to use the high-order guiding-centre equations; for the very highly energetic electrons ($150\,\mathrm {MeV}$ or more), if the pitch angle is nearly zero, the first-order guiding-centre equations are still applicable; for the case where the parallel and perpendicular momenta are both so large that the conservation of $\mu _0$ is broken and the condition of the high-order guiding-centre equations cannot be met, the full model is more accurate.

Figure 10. The greatest variation of magnetic field $\delta \boldsymbol {B}$, the maximum amplitude of the lowest magnetic moment $\delta \mu _0$ and the maximum amplitude of the magnetic moment $\delta \mu$ during one gyro-period. The different colours stand for different energies. The light blue points represent $p_0 =400m_0c$ or $T=168.4\,\mathrm {MeV}$, the dark blue points represent $p_0 =300m_0c$ or $T =126.3\,\mathrm {MeV}$ and the light coral points represent $p_0 =100m_0c$ or $T =41.8\,\mathrm {MeV}$.

Figure 11. Theoretical analysis of $\delta \boldsymbol {B}$ versus $\delta \mu _0$. The different colours stand for different energies. The light blue points represent $p_0 =400m_0c$ or $T=168.4\,\mathrm {MeV}$, the dark blue points represent $p_0 =300m_0c$ or $T =126.3\,\mathrm {MeV}$.

3 Conclusion

In this paper, focusing on synchrotron radiation, extensive numerical verifications of the validity of the high-order guiding-centre theory are given for simulations involving REs. The electron energy keeps growing until the radiation loss is balanced by the electric field acceleration. First, we push a test particle by applying the full orbit model in CFETR. The non-conservation of the lowest magnetic moment $\mu _0$ means the breakdown of the first-order guiding-centre model, but this has been addressed by considering the higher-order correction terms depending on $p_\parallel$ and the magnetic-field curvature. Then, we mainly study the magnetic moment and the perpendicular momentum. The results show that the new expression $\mu$ is still a good adiabatic invariant. After the analysis of the relation between the variation of the magnetic field $\delta \boldsymbol {B}$ and the conservation of the magnetic moment, we show that the breakdown of the conservation of $\mu _0$ is not merely the result of the large variation of $\delta \boldsymbol {B}$, but is the combination of $p_\parallel$ and the large variation of $\delta \boldsymbol {B}$. We analyse the applicable conditions of the high-order guiding-centre equations. By dividing the perpendicular momentum into $\boldsymbol {p}_{\bot }=\boldsymbol {p}_{\bot 0}+\boldsymbol {p}_C$, where $\boldsymbol {p}_{\bot 0}$ is the vector spanning a symmetric circle around $\boldsymbol {b}$ and $\boldsymbol {p}_C$ is the curvature drift, we explain the details of the advancement of the perpendicular momentum, especially the extremum and the amplitude, and derive the high-order guiding-centre equations in Appendix A.

The significant reason why we choose the high-order guiding-centre equations instead of the first-order guiding-centre equations lies in the radiation of the gyro-motion. Taking the synchrotron radiation into consideration, we simulate the energy saturation process of an electron under the high-order guiding-centre model. It is found that the HGSA can effectively produce the collisionless pitch-angle scattering, which dominantly affects synchrotron radiation. Synchrotron radiation plays an essential role in the life cycle of REs. The energy of REs first increases and then becomes saturated until the electric field acceleration is balanced by the radiation dissipation, and HGSA can precisely represent this process. Since the time scale of the energy saturation of REs is close to seconds, the computational cost becomes significant. In order to save costs, an analytical estimate is derived for the time it takes for synchrotron drag to balance an accelerating electric field and the provided formula has been numerically verified. Test calculations reveal that HGSA is favourable for exploiting the dynamics of REs in tokamak plasmas.