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. Reference Sommariva, Nardon, Beyer, Hoelzl, Huijsmans and van Vugt2017; Hao, Chen & Cai Reference Hao, Chen, Cai, Li, Wang, Wu, Wang, Chen, Wang and Gao2020; Zhao, Wang & Wang Reference Zhao, Wang, Wang, Hao and CFETR2020). 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. Reference Liu, Wang and Qin2016).

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. (Reference Liu, Wang and Qin2016). 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 Reference Martín-Solís, Sánchez and Esposito1999; Liu et al. Reference Liu, Qin, Wang, Yang, Zheng, Yao, Zheng, Liu and Liu2016). 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 Reference Burby, Squire and Qin2013). In Liu et al. (Reference Liu, Qin, Hirvijoki, Wang and Liu2018b), 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 (Reference Wang, Qin and Liu2016), 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 Reference Martín-Solís, Sánchez and Esposito2010) show that the synchrotron radiation increases the threshold electric field. For a relativistic particle, the radiation-reaction force (RR-force) (Pauli Reference Pauli1958) 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. Reference Hirvijoki, Decker, Brizard and Embréus2015)

(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. Reference Hirvijoki, Decker, Brizard and Embréus2015).

Based on the non-canonical Hamiltonian mechanics approach using the Lie perturbation method, the high-order relativistic guiding-centre equations (Liu et al. Reference Liu, Qin, Hirvijoki, Wang and Liu2018b) 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. Reference Martín-Solís, Alvarez, Sánchez and Esposito1998). 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}$.