FOXTAIL: Modeling the nonlinear interaction between
Alfvén eigenmodes and energetic particles in tokamaks

The equations used in FOXTAIL are based on a Lagrangian formulation of the wave–particle system [4]. Each simulation is formulated as an initial value problem, starting from an energetic particle distribution in a background equilibrium plasma and a set of eigenmodes with a static spatial structure and a dynamical amplitude and phase. Eigenmodes are treated as weak perturbations of the equilibrium, excluding direct mode–mode interaction. The nonlinear coupling of eigenmodes is taken into account only by the energy and momentum exchange between the modes via the energetic particles.

The physical process of the considered wave–particle interactions is essentially absorption and stimulated emission, and the interactions are energy and momentum conserving. The total energy of the eigenmode is proportional to the amplitude squared. Energy conservation is ensured by the equation

$$\sum_{\mathrm{particles}}\dot{W}_{\mathrm{wave}}+\sum_{\mathrm{eigenmodes}}C\operatorname{Re}(A\dot{A}^{*})=0,$$

(1)

where $\dot{W}_{\mathrm{wave}}$ is the time derivative of the kinetic particle energy, as accelerated by the wave field, $A$ is the complex amplitude of the eigenmode, and $C$ is the ratio between the eigenmode energy and $|A|^{2}/2$. A Lagrangian formulation of the wave–particle interaction is presented in section 2.7, which is shown to be consistent with eq. (1).

The acceleration of a particle by an Alfvén eigenmode convects the momentum of the particles along curves in the phase space $(W,P_{\phi},\mu)$ according to

$$\frac{\mathrm{d}W}{\omega}=\frac{\mathrm{d}P_{\phi}}{n},\leavevmode\nobreak\ \mathrm{d}\mu=0,$$

(2)

where $W$ is the particle kinetic energy, $P_{\phi}$ is the toroidal canonical momentum, $\mu$ is the magnetic moment, $n$ is the toroidal mode number of the eigenmode and $\omega$ is the eigenmode frequency. The magnetic moment is unperturbed, since the Alfvén eigenmodes are low frequency waves ($\omega\ll\omega_{\mathrm{c}}$). The curves in $W,P_{\phi}$-space specified by eq. (2) are referred to as the characteristic curves of wave–particle interaction. Within certain parameter limits, FOXTAIL is equivalent with a one-dimensional bump-on-tail model describing the wave–particle interaction (cf. Ref. [12]). Different characteristic curves are indistinguishable in this limit, and the momentum of the 1D model then corresponds to a variable indexing the location of the particle along the characteristic curves. From the theory of the 1D bump-on-tail model, it is apparent that an inverted energy distribution along the characteristic curves at the wave–particle resonance can excite eigenmodes via a process analogous to Landau damping.

The FOXTAIL code is essentially split into two parts, “FOX” and “TAIL”, as illustrated in Fig. 1. TAIL is the numerical dynamics solver of the eigenmodes and the energetic particles, solving the wave–particle system as an initial value problem. FOX can be viewed as a preprocessor of TAIL, calculating orbital, eigenmode and interaction related data for a given set of eigenmodes and energetic particle species in a defined equilibrium configuration. The distribution of energetic particles is represented using a finite set of markers with predefined weights.

The phase space of markers in TAIL is described using action–angle coordinates of the equilibrium system [3]. The equations of motion of the system become particularly simple in this choice of coordinates, which contributes to the computational efficiency of the solver. In the absence of wave field perturbations, the “action” coordinates (i.e. the momentum coordinates of the canonical action–angle coordinate system) are constants of motion, whereas the “angle” coordinates (configuration space coordinates) evolve with constant velocities. The presence of eigenmodes perturbs this simple dynamics, and the adiabatic invariants are convected according to eq. (2).

The angle coordinates of the system are given by

$$\left\{\begin{array}[]{*3{>{\displaystyle\vspace{1mm}}l}}\vskip 2.84526pt\tilde{\alpha}=\alpha+\int_{0}^{\theta}\mathrm{d}\theta^{\prime}\frac{\bar{\omega}_{\mathrm{c}}(\bm{P})-\omega_{\mathrm{c}}(\bm{P},\theta^{\prime})}{\dot{\theta}(\bm{P},\theta^{\prime})},\\ \vskip 2.84526pt\tilde{\theta}=\int_{0}^{\theta}\mathrm{d}\theta^{\prime}\frac{\omega_{\mathrm{B}}(\bm{P})}{\dot{\theta}(\bm{P},\theta^{\prime})},\\ \vskip 2.84526pt\tilde{\phi}=\phi+\int_{0}^{\theta}\mathrm{d}\theta^{\prime}\frac{\omega_{\mathrm{p}}(\bm{P})-\dot{\phi}(\bm{P},\theta^{\prime})}{\dot{\theta}(\bm{P},\theta^{\prime})},\end{array}\right.$$

(3)

where $\alpha$ is the gyro-angle, $\theta$ and $\phi$ are the poloidal and toroidal angles, respectively, $\bm{P}\equiv(P_{\alpha},P_{\theta},P_{\phi})$ are the action coordinates, canonical to the angles $(\tilde{\alpha},\tilde{\theta},\tilde{\phi})$. Furthermore, $\omega_{\mathrm{c}}$ ($\bar{\omega}_{\mathrm{c}}$) is the (time averaged) gyro-frequency and $\omega_{\mathrm{B}}$ and $\omega_{\mathrm{p}}$ are the bounce and precession frequencies, respectively. The integrals of eq. (3) are evaluated on the poloidal coordinates along the guiding center orbits.

FOXTAIL uses the adiabatic invariants $\bm{J}\equiv(\mu,\Lambda,P_{\phi})$ as momentum coordinates ($\Lambda=\mu B_{0}/W$ is the normalized magnetic moment, where $B_{0}$ is the on-axis magnetic field strength). These momentum coordinates can be expressed as functions of $\bm{P}$. The angle coordinates $(\tilde{\alpha},\tilde{\theta},\tilde{\phi})$ are referred to as the transformed gyro-angle, poloidal angle and toroidal angle, respectively. In the equilibrium system, where the momentum coordinates $\bm{J}$ are constant, the transformed angles evolve at a constant velocity $(\bar{\omega}_{\mathrm{c}},\omega_{\mathrm{B}},\omega_{\mathrm{p}})$. The dynamics in the FOXTAIL model is averaged over gyration time scales, and consequently $\tilde{\alpha}$ is an ignorable coordinate of the system.

The transformed poloidal angle, $\tilde{\theta}$, can be viewed as an index of the location of the particle in the guiding center orbit, where $\tilde{\theta}:0\to 2\pi$ is a complete period of the guiding center orbit in the poloidal plane ($\tilde{\theta}=0$ is defined as the point where the outer leg of the orbit intersects with the equatorial plane). In section 2.5, it is shown that a Fourier expansion of the instant acceleration of the particle in the wave field is a convenient representation of the wave–particle interactions.

Complications arise for the $\tilde{\theta}$ coordinate close to the boundary $\omega_{\mathrm{B}}=0$, where the particle asymptotically approaches one of the turning points. In this limit, all points along the orbit besides the turning points are represented by infinitely narrow intervals in $\tilde{\theta}$, and the representation of the wave–particle interaction using Fourier series expansions in $\tilde{\theta}$-space becomes invalid. These complications can potentially be resolved either by ad hoc boundary conditions or by more sophisticated coordinate transformations close to this boundary. However, all of the numerical studies presented in this paper consider scenarios where the simulated energetic particle distributions are sufficiently far from the boundary.

FOX contains subroutines that solves the 3D motion of particles in a given equilibrium field configuration on a grid in $\bm{J}$-space. The equilibrium configuration and the particle orbits are described in $\psi,\theta,\phi$-space, where $\psi$ is the poloidal magnetic flux per radian. For each orbit on the $\bm{J}$-grid, the time evolution of $\psi$, $\theta$ and $\phi_{\mathrm{s}}$ is calculated, where $\phi_{\mathrm{s}}\equiv\phi-\tilde{\phi}+\omega_{\mathrm{p}}\tilde{\theta}/\omega_{\mathrm{B}}$ is the shifted toroidal coordinate ($\phi_{\mathrm{s}}$ is $\phi$ shifted such that $\phi_{\mathrm{s}}=0$ coincides with $\tilde{\theta}=0$). The equilibrium configuration is taken to be axisymmetric. Assuming MHD force balance and nested magnetic flux surfaces, the equilibrium magnetic field can be expressed as

$$\bm{B}=F(\psi)\nabla\phi+\nabla\phi\times\nabla\psi.$$

(4)

The guiding center motion consists of a parallel motion and a drift motion, where the drift is given by the combined $\bm{E}\times\bm{B}$, $\nabla B$ and curvature drifts according to

$$\bm{v}_{\mathrm{d}}=\frac{\bm{E}\times\bm{B}}{B^{2}}+\frac{\mu(2B_{0}-\Lambda B)}{Ze\Lambda B^{3}}\bm{B}\times\nabla B.$$

(5)

By combining $W=mv_{\parallel}^{2}/2+\mu B$ with $P_{\phi}=mv_{\parallel}F/B+mv_{\mathrm{d},\phi}-Ze\psi$ and eliminating $v_{\parallel}$, it can be shown that the coordinates of the guiding center orbits in the poloidal plane follow the equation

$$f(\psi,\theta,\bm{J})=0,$$

(6)

where

	$\displaystyle f(\psi,\theta,\bm{J})$	$\displaystyle\equiv 1-\frac{\Lambda B(\psi,\theta)}{B_{0}}-\frac{\Lambda B^{2}(\psi,\theta)}{2m\mu B_{0}F^{2}(\psi)}$
		$\displaystyle\quad\,\times[P_{\phi}+Ze\psi-mv_{\mathrm{d},\phi}(\psi,\theta,\mu,\Lambda)]^{2},$		(7)

	$\displaystyle v_{\mathrm{d},\phi}(\psi,\theta,\mu,\Lambda)$	$\displaystyle=\frac{E^{\psi}(\psi,\theta)}{B^{2}(\psi,\theta)}-\frac{\mu[2B_{0}-\Lambda B(\psi,\theta)]}{Ze\Lambda B^{3}(\psi,\theta)}$
		$\displaystyle\quad\,\times[\nabla B(\psi,\theta)]^{\psi}$		(8)

is the covariant toroidal component of the drift velocity, $E^{\psi}$ and $[\nabla B]^{\psi}$ are the contravariant $\psi$-components of $\bm{E}$ and $\nabla B$, respectively, and $m$ is the particle mass.

Once the projections of the orbits on the poloidal plane are solved on the chosen $\bm{J}$-grid using eq. (6), the corresponding time coordinates are calculated according to

$$t(\psi)=\int_{\psi_{0}}^{\psi}\frac{\mathrm{d}\psi}{\dot{\psi}},$$

(9)

where

$$\dot{\psi}=\frac{JFE_{\theta}-g_{\theta\theta}E_{\phi}}{J^{2}B^{2}}-\frac{\mu F(2B_{0}-\Lambda B)}{Ze\Lambda JB^{3}}\frac{\partial B}{\partial\theta},$$

(10)

$$J\equiv\left(\frac{\partial\bm{r}}{\partial\psi}\times\frac{\partial\bm{r}}{\partial\theta}\right)\cdot\frac{\partial\bm{r}}{\partial\phi},$$

(11)

and $\psi(t=0)=\psi_{0}$ is defined as the point where the outer leg of the guiding center orbit intersects the equatorial plane ($\tilde{\theta}=0$). Similarly, the $\phi_{\mathrm{s}}$ coordinates are calculated from the toroidal velocity according to

$$\phi_{\mathrm{s}}(t)=\int_{0}^{t}\mathrm{d}t^{\prime}\>\dot{\phi}(\bm{J},t^{\prime}),$$

(12)

where

$$\dot{\phi}(\bm{J},t)=\frac{P_{\phi}+Ze\psi(\bm{J},t)}{mR^{2}(\bm{J},t)}.$$

(13)

When calculating the wave–particle interaction coefficients, the poloidal velocity is required at each point of the orbit (see eq. (22)). It is given by

	$\displaystyle\dot{\theta}$	$\displaystyle=\frac{P_{\phi}+Ze\psi-mv_{\mathrm{d},\phi}}{mJF}-\frac{JFE_{\psi}+g_{\psi\theta}E_{\phi}}{J^{2}B^{2}}$
		$\displaystyle\quad\,+\frac{\mu F(2B_{0}-\Lambda B)}{Ze\Lambda JB^{3}}\frac{\partial B}{\partial\psi}.$		(14)

The present version of FOXTAIL describes the dynamics of low frequency, shear eigenmodes, such as the toroidal Alfvén eigenmodes (TAEs), but the model can be extended to describe eigenmodes with, e.g., compressible components as well.

Neglecting plasma resistivity, the general electric wave field can be represented by the two scalar potentials $\Phi$ and $\Psi$ according to

$$\delta\bm{E}=-\nabla_{\perp}\Phi+\frac{\bm{B}\times\nabla\Psi}{B},$$

(15)

where the first term is associated with magnetic shear, and the second term is associated with magnetic compression. When parallel gradients in the plasma are negligible in comparison to perpendicular gradients, the excitation of the two scalar potentials $\Phi$ and $\Psi$ is almost decoupled, and it is sufficient to describe the shear Alfvén wave using the first term [4].

In an axisymmetric toroidal plasma, the scalar potential of each eigenmode (indexed by $i$) can be written on the form

	$\displaystyle\Phi_{i}(\bm{r},t)$	$\displaystyle=\operatorname{Re}\sum_{m}C_{i}(t)\mathrm{e}^{\mathrm{i}\chi_{i}(t)}\Phi_{i,m}(\psi)$
		$\displaystyle\quad\,\times\mathrm{e}^{\mathrm{i}(n_{i}\phi-m\theta-\omega_{i}t)},$		(16)

where $C_{i}$ and $\chi_{i}$ are the slowly varying amplitude and phase of the eigenmode, respectively ($\dot{C}_{i}/C_{i}\ll\omega_{i}$, $\dot{\chi_{i}}\ll\omega_{i}$), $n_{i}$ is the toroidal mode number and $\omega_{i}$ is the eigenmode frequency. The electric wave field therefore given by

$$\delta\bm{E}=\operatorname{Re}\sum_{i}C_{i}\mathrm{e}^{\mathrm{i}\chi_{i}}\tilde{\bm{E}}_{i}\mathrm{e}^{\mathrm{i}(n_{i}\phi-\omega_{i}t)},$$

(17)

where

	$\displaystyle\tilde{\bm{E}}_{i}$	$\displaystyle=\sum_{m}\mathrm{e}^{-\mathrm{i}m\theta}\bigg{(}\left[\mathrm{i}\Phi_{i,m}g_{\psi\theta}G_{i,m}-\frac{\mathrm{d}\Phi_{i,m}}{\mathrm{d}\psi}\right]\nabla\psi+\mathrm{i}\Phi_{i,m}$
		$\displaystyle\quad\,\times\Big{[}(g_{\theta\theta}G_{i,m}+m)\nabla\theta+(JFG_{i,m}-n_{i})\nabla\phi\Big{]}\bigg{)},$		(18)

$G_{i,m}\equiv(n_{i}JF-mR^{2})/(J^{2}R^{2}B^{2})$.

The acceleration of the energetic particle in the wave field is described by the equation

$$\dot{W}=Ze\bm{v}\cdot\delta\bm{E}\approx Ze\langle\bm{v}\cdot\delta\bm{E}\rangle_{\mathrm{g}},$$

(19)

where $\langle\cdot\rangle_{\mathrm{g}}$ averages over the gyro-motion. Initial versions of FOXTAIL only consider the lowest order averaging over the gyro-motion, but gyro-kinetic corrections to the averaging may be included in later versions. The averaged $\bm{v}\cdot\delta\bm{E}$ can be written as

$$\langle\bm{v}\cdot\delta\bm{E}\rangle_{\mathrm{g}}=\langle\bm{v}\rangle_{\mathrm{g}}\cdot\langle\delta\bm{E}\rangle_{\mathrm{g}}+\frac{\mu}{Ze}\frac{\partial\langle\delta B_{\parallel}\rangle_{S_{\mathrm{g}}}}{\partial t},$$

(20)

where $\langle\delta B_{\parallel}\rangle_{S_{\mathrm{g}}}$ averages the parallel magnetic wave field over the surface $S_{\mathrm{g}}$ enclosed by the gyro-ring (generated by the span of the gyro-angle. For shear waves, the $\delta B_{\parallel}$ term can be neglected, and we are left with

$$\dot{W}=Ze\langle\bm{v}\rangle_{\mathrm{g}}\cdot\langle\delta\bm{E}\rangle_{\mathrm{g}}=\operatorname{Re}\sum_{i}C_{i}\mathrm{e}^{\mathrm{i}\chi_{i}}V_{i}\mathrm{e}^{\mathrm{i}(n_{i}\phi-\omega_{i}t)},$$

(21)

where

$$V_{i}(\bm{J},t)=Ze(\dot{\psi}\tilde{E}_{\psi}+\dot{\theta}\tilde{E}_{\theta}+\dot{\phi}\tilde{E}_{\phi}),$$

(22)

$\dot{\psi}(\bm{J},t)$, $\dot{\theta}(\bm{J},t)$ and $\dot{\phi}(\bm{J},t)$ are the guiding center velocities, given by eqs. (10), (2.4) and (13), respectively, and $\tilde{E}_{\psi}(\psi(t),\theta(t))$, $\tilde{E}_{\theta}(\psi(t),\theta(t))$ and $\tilde{E}_{\phi}(\psi(t),\theta(t))$ are given by eq. (18).

All guiding center coordinates $(\psi,\theta,\phi_{\mathrm{s}})$ and their velocities can be written as functions of $\bm{J}$ and $\tilde{\theta}$. It is then possible to write $\dot{W}$ on a very compact form using action–angle coordinates:

$$\dot{W}=\operatorname{Re}\sum_{i,\ell}C_{i}\mathrm{e}^{\mathrm{i}\chi_{i}}V_{i,\ell}\mathrm{e}^{\mathrm{i}(\ell\tilde{\theta}+n_{i}\tilde{\phi}-\omega_{i}t)},$$

(23)

where

	$\displaystyle V_{i,\ell}(\bm{J})$	$\displaystyle=\frac{Ze}{2\pi}\int_{0}^{2\pi}\mathrm{d}\tilde{\theta}\>V_{i}(\tilde{\theta},\bm{J})$
		$\displaystyle\quad\,\times\exp\bigg{[}\mathrm{i}\bigg{(}n_{i}\left[\phi_{\mathrm{s}}(\tilde{\theta},\bm{J})-\frac{\omega_{\mathrm{p}}(\bm{J})}{\omega_{\mathrm{B}}(\bm{J})}\tilde{\theta}\right]-\ell\tilde{\theta}\bigg{)}\bigg{]}$		(24)

The coefficients $V_{i,\ell}$ are the Fourier series expansion of the wave–particle interaction that named the FOX code. It can be seen in eq. (23) that a wave–particle resonance is characterized by the condition $\ell\omega_{\mathrm{B}}+n_{i}\omega_{\mathrm{p}}\approx\omega_{i}$. The model is very efficient in the sense that one can select the relevant modes $i$ and coefficients $\ell$ that are close to resonant with an ensemble of energetic particles and neglect all other modes and non-resonant Fourier coefficients. If the relevant coefficients only have $\ell=0$, $\tilde{\theta}$ becomes an ignorable coordinate of the system, and the shortest time scales that need to be resolved in simulations using TAIL are the precession time scales, $\omega_{\mathrm{p}}^{-1}$.

As was mentioned in section 2.3, complications arise for the choice of action–angle coordinates when describing the wave–particle interaction in regions where $\omega_{\mathrm{B}}\approx 0$. This is explicitly seen in eq. (24), where $\omega_{\mathrm{p}}/\omega_{\mathrm{B}}\to\infty$, and the integrand oscillates infinitely fast in $\tilde{\theta}$. $V_{i}(\tilde{\theta},\bm{J})$ is bounded, and constant for all $\tilde{\theta}$ in this limit except for the infinitely narrow intervals that do not represent the turning points. For this reason, it can be understood that $V_{i,\ell}$ tends to zero towards the $\omega_{\mathrm{B}}=0$ boundary surface. Note that the sum of eq. (23) over all integers $\ell\in\mathbb{Z}$ should remain finite, given that particles can be accelerated across these surfaces by wave fields.

A model for describing the dynamics of the momentum variables $(\mu,\Lambda,P_{\phi})$ and the amplitudes and phases of the eigenmodes remains to be found. Such a model can be derived from a Lagrangian formulation of the wave–particle system. We consider additions to the equilibrium Lagrangian by power series of the eigenmode amplitude $C_{i}$.

Expressed in action–angle variables, the zeroth order Lagrangian reads [13]

$$\mathcal{L}_{0}=\sum_{k}\left[P_{\alpha,k}\dot{\tilde{\alpha}}_{k}+P_{\theta,k}\dot{\tilde{\theta}}_{k}+P_{\phi,k}\dot{\tilde{\phi}}_{k}-\mathcal{H}_{0}(\bm{P}_{k})\right],$$

(25)

where $k$ is the particle index, and $\mathcal{H}_{0}$, satisfying

$$\frac{\partial\mathcal{H}_{0}}{\partial\bm{P}}=(\bar{\omega}_{\mathrm{c}},\omega_{\mathrm{B}},\omega_{\mathrm{p}}),$$

(26)

is the equilibrium Hamiltonian. A first order Lagrangian consistent with eq. (23) is

$$\mathcal{L}_{1}=\operatorname{Im}\sum_{i,k,\ell}\frac{C_{i}}{\omega_{i}}\mathrm{e}^{\mathrm{i}\chi_{i}}V_{i,\ell}(\bm{J}_{k})\mathrm{e}^{\mathrm{i}(\ell\tilde{\theta}_{k}+n_{i}\tilde{\phi}_{k}-\omega_{i}t)}.$$

(27)

In a low-$\beta$ plasma, the second order Lagrangian for shear Alfvén eigenmodes can be expressed on the form [4]

$$\mathcal{L}_{2}=-\sum_{i}\frac{\dot{\chi}_{i}C_{i}^{2}}{2\mu_{0}\omega_{i}}\int\mathrm{d}V\>\frac{|\tilde{\bm{E}}_{i}(\bm{r})|^{2}}{v_{\mathrm{A}}^{2}(\bm{r})}$$

(28)

when neglecting rapidly oscillating terms (on time scales $\omega_{i}^{-1}$), where $\mu_{0}$ is the vacuum permeability.

Note that the system is invariant under the transformation $C_{i}\mathrm{e}^{\mathrm{i}\chi_{i}}\to\kappa_{i}C_{i}\mathrm{e}^{\mathrm{i}\chi_{i}}$, $\Phi_{i,m}\to\Phi_{i,m}/\kappa_{i}$ for any constant $\kappa_{i}$. A normalization can be imposed by the condition

$$\int\mathrm{d}V\>\frac{|\tilde{\bm{E}}_{i}(\bm{r})|^{2}}{\mu_{0}v_{\mathrm{A}}^{2}(\bm{r})}=1$$

(29)

for each mode $i$. Then the second order Lagrangian reduces to

$$\mathcal{L}_{2}=-\sum_{i}\frac{\dot{\chi}_{i}C_{i}^{2}}{2\omega_{i}}.$$

(30)

Assuming that the amplitudes $C_{i}$ are small enough, such that $|\mathcal{L}_{1}|\ll|\mathcal{H}_{0}|$, the corresponding Hamiltonian for the wave–particle system is $\mathcal{H}=\mathcal{H}_{0}-\mathcal{L}_{1}$, with the canonically conjugate pairs $(\tilde{\alpha}_{k}\,;P_{\alpha,k})$, $(\tilde{\theta}_{k}\,;P_{\theta,k})$, $(\tilde{\phi}_{k}\,;P_{\phi,k})$ and $(-\chi_{i}\,;C_{i}^{2}/2\omega_{i})$. From this, one can derive all the remaining equations of motion of the wave–particle system self-consistently.

Defining the complex amplitude $A_{i}\equiv C_{i}\mathrm{e}^{\mathrm{i}\chi_{i}}$ and using that $\mu=ZeP_{\alpha}/m$, the equations of motion of the wave–particle system is

$$\left\{\begin{array}[]{*3{>{\displaystyle\vspace{1mm}}l}}\vskip 2.84526pt\dot{\mu}_{k}=0,&\vskip 2.84526pt\dot{\tilde{\theta}}_{k}=\omega_{\mathrm{B}}(\bm{J}_{k}),\\ \vskip 2.84526pt\dot{\Lambda}_{k}=-\frac{\Lambda_{k}^{2}}{\mu_{k}B_{0}}\operatorname{Re}\sum_{i}A_{i}U_{i,k},&\vskip 2.84526pt\dot{\tilde{\phi}}_{k}=\omega_{\mathrm{p}}(\bm{J}_{k}),\\ \vskip 2.84526pt\dot{P}_{\phi,k}=\operatorname{Re}\sum_{i}\frac{n_{i}}{\omega_{i}}A_{i}U_{i,k},&\vskip 2.84526pt\dot{A}_{i}=-\sum_{k}U_{i,k}^{*},\end{array}\right.$$

(31)

where

$$U_{i,k}\equiv\sum_{\ell}V_{i,\ell}(\bm{J}_{k})\mathrm{e}^{\mathrm{i}(\ell\tilde{\theta}_{k}+n_{i}\tilde{\phi}_{k}-\omega_{i}t)}.$$

(32)

The effects of particle collisions are not covered by the theory presented in the preceding sections. These effects can be added explicitly to the system, e.g. by using a diffusion operator in momentum space [14] combined with a momentum drag [15] derived from the Fokker–Planck operator, which act directly on the energetic particle distribution. Adding diffusion to the system transforms the system of ordinary differential equations in eq. (31) to a system of stochastic differential equations. The general drag–diffusion operator can be written on Itō form according to

$$\left\{\begin{array}[]{*3{>{\displaystyle\vspace{1mm}}l}}\vskip 2.84526pt\mathrm{d}\mu_{\mathrm{c}}=a_{\mu}\mathrm{d}t+\bm{b}_{\mu}\cdot\mathrm{d}\bm{W}_{t},\\ \vskip 2.84526pt\mathrm{d}\Lambda_{\mathrm{c}}=a_{\Lambda}\mathrm{d}t+\bm{b}_{\Lambda}\cdot\mathrm{d}\bm{W}_{t},\\ \vskip 2.84526pt\mathrm{d}P_{\phi,\mathrm{c}}=a_{P_{\phi}}\mathrm{d}t+\bm{b}_{P_{\phi}}\cdot\mathrm{d}\bm{W}_{t},\end{array}\right.$$

(33)

where $a_{\bm{J}}$ is associated with drag, $\bm{b}_{\bm{J}}$ is associated with diffusion and the components of $\bm{W}_{t}$ are independent Wiener processes in time. Both $a_{\bm{J}}$ and $\bm{b}_{\bm{J}}$ are functions of $\bm{J}$ and $\tilde{\theta}$ in general, but orbit averaged (i.e. $\tilde{\theta}$ averaged) versions of the operators may be used for simplicity (see e.g. Ref. [16]).

There are other processes external to the energetic particle–Alfvén eigenmode system which may be included. Depending on the source of the energetic particle distribution (e.g. neutral beam injection, cyclotron resonance heating or nuclear fusion), operators can be added that continuously supply particles and/or energy to the system. Particle sources are typically modeled by dynamical weights and statistical redistribution of markers. Losses from Bremsstrahlung and cyclotron radiation can be modeled by additions to the $a_{\bm{J}}$ terms in eq. (33).

An additional $-\gamma_{\mathrm{d},i}A_{i}$ term can be added to the amplitude equation ($\dot{A}_{i}$ in eq. (31)) in order to model various damping mechanisms on the eigenmode amplitudes, where $\gamma_{\mathrm{d},i}$ is the damping rate of the $i$:th mode. This damping can, e.g., come from Landau damping in the interaction with thermal particles or damping due to mode conversion. The $-\gamma_{\mathrm{d},i}A_{i}$ term of the amplitude equation is here referred to as explicit wave damping, unlike the Landau damping coming from the interaction with the energetic particle distribution, which arises implicitly from the model equations.

When going from the Lagrangians of eqs. (25), (27) and (30) to eq. (31), it was assumed that $|\mathcal{L}_{1}|\ll|\mathcal{H}_{0}|$, which allowed one to neglect corrections of the canonical momenta coming from the implicit dependence of $\mathcal{L}_{1}$ with respect to $\dot{\tilde{\theta}}$ and $\dot{\tilde{\phi}}$. For a small enough $\mathcal{L}_{1}$, the final equations for $\dot{\tilde{\theta}}$ and $\dot{\tilde{\phi}}$ are given only by the derivatives of the equilibrium Hamiltonian $\mathcal{H}_{0}$ with respect to $P_{\theta}$ and $P_{\phi}$, respectively. However, the above assumptions might not be valid for large amplitude eigenmodes and for processes affecting the energetic particles or the eigenmodes that are not direct wave–particle interaction, and one may have to consider the lowest order correction to $\dot{\tilde{\theta}}$ and $\dot{\tilde{\phi}}$ due to momentum perturbation. The correction can be understood from the fact that $\tilde{\theta}$ maps differently to locations on the guiding center orbit for different $\bm{J}$. Without the correction, one pushes the particle forwards or backwards along the orbit due to different mapping of $\tilde{\theta}$.

Assuming a minimal perturbation of the guiding center position in $R,Z$-space while perturbing $\bm{J}$ by the amount $\mathrm{d}\bm{J}$, it can be shown that $\tilde{\theta}$ should be corrected by the amount $\partial\tilde{\theta}/\partial\bm{J}\cdot\mathrm{d}\bm{J}$, where

$$\frac{\partial\tilde{\theta}}{\partial\bm{J}}=-\frac{\omega_{\mathrm{B}}}{\dot{R}^{2}+\dot{Z}^{2}}\left(\dot{R}\frac{\partial R}{\partial\bm{J}}+\dot{Z}\frac{\partial Z}{\partial\bm{J}}\right),$$

(34)

$R$ is the distance from the symmetry axis and $Z$ is the vertical guiding center position. Using that $\tilde{\phi}=\phi-\phi_{\mathrm{s}}+\omega_{\mathrm{p}}\tilde{\theta}/\omega_{\mathrm{B}}$, it can be shown that

	$\displaystyle\frac{\partial\tilde{\phi}}{\partial\bm{J}}$	$\displaystyle=\frac{\dot{\phi}-\omega_{\mathrm{p}}}{\dot{R}^{2}+\dot{Z}^{2}}\left(\dot{R}\frac{\partial R}{\partial\bm{J}}+\dot{Z}\frac{\partial Z}{\partial\bm{J}}\right).$
		$\displaystyle\quad\,+\frac{\partial}{\partial\bm{J}}\left(\frac{\omega_{\mathrm{p}}}{\omega_{\mathrm{B}}}\tilde{\theta}-\phi_{\mathrm{s}}\right).$		(35)

All the derivatives with respect to $\bm{J}$ on the right hand sides of eqs. (34) and (2.9) are evaluated while keeping $\tilde{\theta}$ constant.

When the $\bm{J}$ perturbations of any process is stochastic ($\bm{b}_{\bm{J}}$ of eq. (33) is nonzero), the angle coordinates become stochastic as well when including the $\partial\tilde{\theta}/\partial\bm{J}$ and $\partial\tilde{\phi}/\partial\bm{J}$ corrections. This is typically how phase decorrelation is introduced to the system. Such an operator was analyzed for the one-dimensional bump-on-tail model in Ref. [17].

To summarize, the model equations used by TAIL are

$$\left\{\begin{array}[]{*3{>{\displaystyle\vspace{1mm}}l}}\vskip 2.84526pt\mathrm{d}\tilde{\theta}_{k}=\omega_{\mathrm{B}}(\bm{J}_{k})\mathrm{d}t+\frac{\partial\tilde{\theta}(\tilde{\theta}_{k},\bm{J}_{k})}{\partial\bm{J}}\cdot\mathrm{d}\bm{J}_{k},\\ \vskip 2.84526pt\mathrm{d}\tilde{\phi}_{k}=\omega_{\mathrm{p}}(\bm{J}_{k})\mathrm{d}t+\frac{\partial\tilde{\phi}(\tilde{\theta}_{k},\bm{J}_{k})}{\partial\bm{J}}\cdot\mathrm{d}\bm{J}_{k},\\ \vskip 2.84526pt\mathrm{d}\mu_{k}=\mathrm{d}\mu_{\mathrm{c}}(\bm{J}_{k}),\\ \vskip 2.84526pt\mathrm{d}\Lambda_{k}=-\frac{\Lambda_{k}^{2}}{\mu_{k}B_{0}}\operatorname{Re}\sum_{i}A_{i}U_{i,k}\mathrm{d}t+\mathrm{d}\Lambda_{\mathrm{c}}(\bm{J}_{k}),\\ \vskip 2.84526pt\dot{P}_{\phi,k}=\operatorname{Re}\sum_{i}\frac{n_{i}}{\omega_{i}}A_{i}U_{i,k}+\mathrm{d}P_{\phi,\mathrm{c}}(\bm{J}_{k}),\\ \vskip 2.84526pt\dot{A}_{i}=-\sum_{k}w_{k}U_{i,k}^{*}-\gamma_{\mathrm{d},i}A_{i},\end{array}\right.$$

(36)

where $k$ is now the marker index with weight $w_{k}$, and $(\mathrm{d}\mu_{\mathrm{c}},\mathrm{d}\Lambda_{\mathrm{c}},\mathrm{d}P_{\phi,\mathrm{c}})$ represents additional differential operators acting on the momentum space of markers. The total wave–particle energy of the system can be defined as

	$\displaystyle W_{\mathrm{tot}}$	$\displaystyle\equiv\sum_{k}w_{k}W_{k}+\sum_{i}\frac{|A_{i}|^{2}}{2}$
		$\displaystyle=\sum_{k}w_{k}\frac{\mu_{k}B_{0}}{\Lambda_{k}}+\sum_{i}\frac{|A_{i}|^{2}}{2}.$		(37)

In the absence of explicit wave damping and particle sources and sinks, it can easily be shown that

$$\dot{W}_{\mathrm{tot}}=-\sum_{k}\frac{w_{k}\mu_{k}B_{0}}{\Lambda_{k}^{2}}\dot{\Lambda}_{k}+\operatorname{Re}\sum_{i}A_{i}\dot{A}_{i}^{*}=0,$$

(38)

which is consistent with the condition for energy conservation in eq. (1).

One of the possible applications of FOXTAIL is to determine whether a one-dimensional bump-on-tail approximation of the system is sufficient to capture the essential wave–particle dynamics, such as growth rate and saturation amplitude. A higher computational efficiency typically follows from the lower dimensionality of the approximation. However, the 1D bump-on-tail model is only valid within certain parameter regimes. Section 4.1 presents a quantitative numerical study of these regimes.

The presence of multiple eigenmodes proved to have negligible effect on the system in the scenarios presented in section 4.1.3 due to the low interaction coefficient of the second mode in the considered part of $\Lambda,P_{\phi}$-space, compared to the interaction coefficient of the first mode. Scenarios with significant multimode dynamics can be constructed by adding an ad hoc scaling factor to the interaction coefficient of the second mode. Multiplying $V_{2,1}$ by a factor of 10.1 gives approximately the same linear growth rate of the two modes. This has been done for cases #5 and #6, presented in Fig. 7, along with the previous case #2. The three scenarios are the same, except that in cases #2 and #6 the eigenmodes are simulated individually, whereas in case #5 both modes are included to the system. Such a comparison allows one to isolate the nonlinear effects of mode interaction via the energetic particle distribution.

Comparing the multimode scenario with the scenarios where the modes are included individually, it can be seen that the indirect interaction between the modes via the energetic particles has significant macroscopic effects on the system. This can partly be understood as a consequence of stochastization of particle trajectories in phase space due to resonance overlap of the two eigenmodes. Stochastization of trajectories causes a locally enhanced transport of energetic particles around the resonances, allowing the eigenmodes to exhaust more energy from the energetic particle distribution (see e.g. stochastization from resonance overlap, [29, 30], and from phase decorrelation, [17]). This results in a wider portion of the energetic particle distribution being flattened around the resonances, compared to the individual eigenmode cases, as seen in Fig. 7.b.

The resonance width of an eigenmode can be estimated as the separatrix width of the unperturbed mode (i.e. in the absence of other modes) along the characteristic curve in $W$-space, using the 1D bump-on-tail approximation. When the resonance widths of the two modes are comparable, the resonance-overlap parameter [31] (English translation: Ref. [32]), estimated as the average resonance width of the two modes divided by their distance in phase space, is an approximate measure of the level of stochastization of particle trajectories. The Chirikov criterion for stochastization is satisfied when the resonance-overlap parameter is larger than unity. The full separatrix width in $W$-space, $W_{\mathrm{sep}}$, is $4\omega_{\mathrm{b}}/\Omega^{\prime}(\bm{K}_{\mathrm{res}})$, with $\omega_{\mathrm{b}}$ given by eq. (47). As seen in Fig. 7.c, the Chirikov criterion is well satisfied for case #5 after $t\gtrsim$1.3\text{\,}\mathrm{m}\mathrm{s}$$.

Slightly different scenarios are tested in the simulations presented in Fig. 8. Cases #7, #8 and #9 are the same as cases #2, #5 and #6, respectively, except that the frequency separation between the two eigenmodes is artificially increased by a factor of 3, and the width of the energetic distribution function along the characteristic curves of the first mode is increased to $\sigma=$4\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}$$ instead of $3\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}$. The scaling factor of the interaction coefficient is also adjusted such that the linear growth rates of the two modes approximately match. As can be seen in Fig. 8.c, the Chirikov criterion is never satisfied for case #8, although the resonance-overlap parameter is of the order of unity. Comparing the amplitude evolutions in Figs. 7.a and 8.a, the two modes of case #8 evolves more similarly to the corresponding individual mode scenarios than case #5 does. This is especially seen in Fig. 8.b, where the modes of case #8 flattens two separate regions of the energetic distribution function, matching the flattening regions of the individual mode scenarios.