2 The gyrokinetic model

The EGAM problem is investigated here with the global nonlinear GK particle-in-cell code ORB5. ORB5 was written for studying electrostatic turbulence in tokamak plasmas (Jolliet et al. Reference Jolliet, Bottino, Angelino, Hatzky, Tran, Mcmillan, Sauter, Appert, Idomura and Villard2007), and extended to treat multiple kinetic species (i.e. thermal ions, electrons, EP, impurities, etc.) and electromagnetic perturbations (Bottino et al. (Reference Bottino, Vernay, Scott, Brunner, Hatzky, Jolliet, McMillan, Tran and Villard2011), Bottino & Sonnendrücker (Reference Bottino and Sonnendrücker2015)). A collision operator is also implemented in ORB5, for the linearized inclusion of inter-species and like-species collisions. In this article, electrostatic collisionless simulations are performed with ORB5. Only the wave–particle nonlinearity is considered, by filtering out all $n\neq 0$ modes and pushing only the EP species along the perturbed trajectories, whereas the bulk ions and electrons follow the unperturbed trajectories) (see (2.1)–(2.4)).

The model equations of the electrostatic version of ORB5 are the marker trajectories, and the gyrokinetic Poisson law regulating the evolution in time of the scalar potential $\unicode[STIX]{x1D719}$ . These equations are derived in a Lagrangian formulation (Tronko, Bottino & Sonnendrücker Reference Tronko, Bottino and Sonnendrücker2016). The equations for the marker trajectories for the thermal ions and fast ions (in the electrostatic version of the code) are (Bottino & Sonnendrücker Reference Bottino and Sonnendrücker2015):

(2.1) $$\begin{eqnarray}\displaystyle & \displaystyle \dot{\boldsymbol{R}}=\frac{1}{m_{s}}p_{\Vert }\frac{\boldsymbol{B}^{\ast }}{B_{\Vert }^{\ast }}+\frac{c}{q_{s}B_{\Vert }^{\ast }}\boldsymbol{b}\times [\unicode[STIX]{x1D707}\unicode[STIX]{x1D735}B+q_{s}\unicode[STIX]{x1D735}\tilde{\unicode[STIX]{x1D719}}] & \displaystyle\end{eqnarray}$$

(2.2) $$\begin{eqnarray}\displaystyle & \displaystyle \dot{p_{\Vert }}=-\frac{\boldsymbol{B}^{\ast }}{B_{\Vert }^{\ast }}\boldsymbol{\cdot }[\unicode[STIX]{x1D707}\unicode[STIX]{x1D735}B+q_{s}\unicode[STIX]{x1D735}\tilde{\unicode[STIX]{x1D719}}] & \displaystyle\end{eqnarray}$$

(2.3) $$\begin{eqnarray}\displaystyle & \displaystyle \dot{\unicode[STIX]{x1D707}}=0. & \displaystyle\end{eqnarray}$$

The set of coordinates used for the phase space is $(\boldsymbol{R},p_{\Vert },\unicode[STIX]{x1D707})$ , i.e. respectively the gyrocentre position, canonical parallel momentum $p_{\Vert }=m_{s}v_{\Vert }$ and magnetic momentum $\unicode[STIX]{x1D707}=m_{s}v_{\bot }^{2}/(2B)$ (with $m_{s}$ and $q_{s}$ being the mass and charge of the species); $v_{\Vert }$ and $v_{\bot }$ are respectively the parallel and perpendicular component of the particle velocity. We label the gyroaverage operator with the tilde symbol $\tilde{\,}$ . $\boldsymbol{B}^{\ast }=\boldsymbol{B}+(c/q_{s})\unicode[STIX]{x1D735}\times (\boldsymbol{b}p_{\Vert })$ , where $\boldsymbol{B}$ and $\boldsymbol{b}$ are the equilibrium magnetic field and magnetic unitary vector.

As we are not interested here in comparing the EGAM damping/growth rates with experimental observations, but only in investigating a specific piece of the nonlinear physics of EGAMs, we neglect kinetic effects of the electrons. This is done by calculating the electron gyrocentre density directly from the value of the scalar potential as (Bottino & Sonnendrücker Reference Bottino and Sonnendrücker2015):

(2.4) $$\begin{eqnarray}\displaystyle n_{e}(\boldsymbol{R},t)=n_{e0}+\frac{en_{e0}}{T_{e0}}(\unicode[STIX]{x1D719}-\bar{\unicode[STIX]{x1D719}}), & & \displaystyle\end{eqnarray}$$

where $\bar{\unicode[STIX]{x1D719}}$ is the flux-surface-averaged potential, instead of treating the electrons with markers evolved with (2.1)–(2.3).

We are interested here in the dynamics of zonal perturbations, and we filter out all non-zonal components. We are not interested here in the wave–wave coupling, and therefore we push the bulk-ion and electron markers along unperturbed trajectories. This means that, in (2.1)–(2.3) for the bulk ions, the last terms, proportional to the EGAM electric field, are dropped. On the other hand, we push the EP markers along the trajectories which include perturbed terms associated with the EGAM electric field. This means that the EP markers are evolved with (2.1)–(2.3) where the terms proportional to the EGAM electric field are retained. This retains the wave–particle nonlinearity of the EP.

Finally, the gyrokinetic Poisson equation is (Bottino & Sonnendrücker (Reference Bottino and Sonnendrücker2015)):

(2.5) $$\begin{eqnarray}\displaystyle -\mathop{\sum }_{s\neq e}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\frac{n_{0s}m_{s}c^{2}}{B^{2}}\unicode[STIX]{x1D735}_{\bot }\unicode[STIX]{x1D719}=\mathop{\sum }_{s\neq e}\int \text{d}W_{s}q_{s}\tilde{\unicode[STIX]{x1D6FF}f_{s}}+q_{e}n_{e}(\boldsymbol{R},t), & & \displaystyle\end{eqnarray}$$

with $n_{0s}=\int \text{d}W_{s}f_{0s}$ . The summation over the species is performed for the bulk ions and for the EP, whereas the electron contribution is given by $q_{e}n_{e}(\boldsymbol{R},t)$ . Here $\unicode[STIX]{x1D6FF}f_{s}=f_{s}-f_{0s}$ is the gyrocentre perturbed distribution function, with $f_{s}$ and $f_{0s}$ being the total and equilibrium (i.e. independent of time, assumed here to be a Maxwellian) gyrocentre distribution functions. The integrals are over the phase-space volume. The phase-space volume element is $\text{d}W_{s}=(2\unicode[STIX]{x03C0}/m_{s}^{2})B_{\Vert }^{\ast }\,\text{d}p_{\Vert }\,\text{d}\unicode[STIX]{x1D707}$ , being the velocity-space infinitesimal volume element. In this article, finite-Larmor-radius effects are considered, for both thermal and fast ions.

3 Equilibrium and linear EGAM dynamics

We consider here the same tokamak configuration adopted in Biancalani et al. (Reference Biancalani, Chavdarovski, Qiu, Bottino, Del Sarto, Ghizzo, Gürcan, Morel and Novikau2017), where the scalings of the EGAM nonlinear saturation levels were studied. The tokamak magnetic equilibrium is defined by major and minor radii of $R_{0}=1~\text{m}$ and $a=0.3125~\text{m}$ , a magnetic field on axis of $B_{0}=1.9~\text{T}$ , a flat safety factor radial profile, with $q=2$ , and circular flux surfaces (with no Grad–Shafranov shift). Flat temperature and density profiles are considered at the equilibrium. The bulk plasma temperature is defined by $\unicode[STIX]{x1D70C}^{\ast }=\unicode[STIX]{x1D70C}_{s}/a$ , with $\unicode[STIX]{x1D70C}_{s}=c_{s}/\unicode[STIX]{x1D6FA}_{i}$ , with $c_{s}=\sqrt{T_{e}/m_{i}}$ being the sound speed. We choose $\unicode[STIX]{x1D70C}^{\ast }=1/128=0.0078$ ( $\unicode[STIX]{x1D70F}_{e}=T_{e}/T_{i}=1$ for all cases described in this article), corresponding to $2/\unicode[STIX]{x1D70C}^{\ast }=256$ .

In the case of a hydrogen plasma, we get a value of the ion cyclotron frequency of $\unicode[STIX]{x1D6FA}_{i}=1.82\times 10^{8}~\text{rad}~\text{s}^{-1}$ and a temperature of $T_{i}=2060~\text{eV}$ . The sound frequency is defined as $\unicode[STIX]{x1D714}_{s}=2^{1/2}v_{ti}/R$ (with $v_{ti}=\sqrt{T_{i}/m_{i}}$ , which for $\unicode[STIX]{x1D70F}_{e}=1$ reads $v_{ti}=c_{s}$ ). We obtain $c_{s}=4.44\times 10^{5}~\text{m}~\text{s}^{-1}$ . This corresponds to the following value of the sound frequency: $\unicode[STIX]{x1D714}_{s}=6.28\times 10^{5}~\text{rad}~\text{s}^{-1}$ .

The energetic-particle distribution function is a double bump-on-tail, with two bumps at $v_{\Vert }=\pm v_{\text{bump}}$ , as in Biancalani et al. (Reference Biancalani, Bottino, Lauber and Zarzoso2014, Reference Biancalani, Chavdarovski, Qiu, Bottino, Del Sarto, Ghizzo, Gürcan, Morel and Novikau2017), labelled here as $F_{EP}$ . In this article, $v_{\text{bump}}=4v_{ti}$ is chosen. In order to initialize an EP distribution function which is function of the constants of motion only, we neglect the radial dependence of the magnetic field in $v_{\Vert }(\unicode[STIX]{x1D707},E,\boldsymbol{R})=\sqrt{2(E-\unicode[STIX]{x1D707}B)/m/v_{ti}}\simeq \tilde{v}_{\Vert }(\unicode[STIX]{x1D707},E)$ in the Vlasov equation (details are given in Zarzoso et al. Reference Zarzoso, Biancalani, Bottino, Lauber, Poli, Girardo, Garbet and Dumont2014). Neumann and Dirichlet boundary conditions are imposed to the scalar potential, respectively at the inner and outer boundaries, $s=0$ and $s=1$ .

Figure 1.

Linear frequency (a) and growth rate (b) versus EP concentration.

The scan of linear simulations with different EP concentration, performed in Biancalani et al. (Reference Biancalani, Chavdarovski, Qiu, Bottino, Del Sarto, Ghizzo, Gürcan, Morel and Novikau2017), is reported here for completeness. This defines the linear dynamics of the system. The dependence of the linear frequency and growth rate on the EP concentration is shown in figure 1. For comparison, the GAM frequency for these parameters is $\unicode[STIX]{x1D714}_{\text{GAM}}=1.8\unicode[STIX]{x1D714}_{s}$ . The list of the parameters described above defines a regime where the EGAM comes continuously from a GAM, when increasing the EP concentration from zero. On the other hand, the EGAM would come from a Landau pole, for sufficiently larger values of $q$ (see Zarzoso et al. Reference Zarzoso, Biancalani, Bottino, Lauber, Poli, Girardo, Garbet and Dumont2014).

4 Nonlinear EGAM evolution

In this section, we describe the evolution in time of the nonlinear simulations of EGAMs performed with ORB5. Here, as in the rest of this article, the wave–particle nonlinearity only is considered. As an example, we consider a case with $n_{EP}/n_{i}=0.10$ . A zonal (i.e. axisymmetric) radial electric field is initialized at $t=0$ , with an amplitude of the order of $10^{3}~\text{V}~\text{m}^{-1}$ , and let evolve in time in a nonlinear simulation with ORB5. A typical simulation has a spatial resolution set by a grid of (ns,nchi,nphi) $=$ (256, 64, 16) number of points respectively in the radial, poloidal and toroidal direction, a time step of $\text{d}t=20\unicode[STIX]{x1D6FA}_{i}^{-1}$ , and a number of markers of $(\text{ntot}_{i},\text{ntot}_{EP})=(10^{7},10^{7})$ respectively for the thermal and fast ions. An initial linear phase is observed, where the radial electric field grows exponentially in time. In this phase the linear frequency and growth rate are measured and checked to match with the ones of the linear simulation: $\unicode[STIX]{x1D714}_{L}=1.24\unicode[STIX]{x1D714}_{s}$ , $\unicode[STIX]{x1D6FE}_{L}=0.06\unicode[STIX]{x1D714}_{s}$ . Then, a nonlinear phase is entered, the growth rate gradually decreases to zero, and the radial electric field saturates at $t\simeq 2.5\times 10^{4}\unicode[STIX]{x1D6FA}_{i}^{-1}$ (see figure 2 a), when the electric field reaches a value of $\unicode[STIX]{x1D6FF}E_{r}\simeq 3.5\times 10^{4}~\text{V}~\text{m}^{-1}$ . This value of the saturated electric field can be compared with the prediction of Biancalani et al. (Reference Biancalani, Chavdarovski, Qiu, Bottino, Del Sarto, Ghizzo, Gürcan, Morel and Novikau2017):

(4.1) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FF}E_{r,th}=\frac{2RB\unicode[STIX]{x1D6FD}_{0}^{2}}{\unicode[STIX]{x1D714}_{\text{GAM}}}\unicode[STIX]{x1D6FE}_{L}^{2}=3.5\times 10^{4}\frac{\text{V}}{\text{m}}, & & \displaystyle\end{eqnarray}$$

where the constant $\unicode[STIX]{x1D6FD}_{0}=2.66$ is estimated in Biancalani et al. (Reference Biancalani, Chavdarovski, Qiu, Bottino, Del Sarto, Ghizzo, Gürcan, Morel and Novikau2017) for this regime. For the present equilibrium plasma profiles, the regime is defined by the value of $q=2$ , for which we have an EGAM which comes from a GAM (see also § 3). We emphasize here that the quadratic scaling of the electric field with the linear growth rate shown in (4.1) has been found to be valid for the whole considered range of EP concentrations (the same range used in figure 1). After the saturation, the EGAM enters a deep nonlinear phase ( $t>2.5\times 10^{4}\unicode[STIX]{x1D6FA}_{i}^{-1}$ ), when the electric field starts decreasing in amplitude. In this article we are interested in the first nonlinear phase only, up to the saturation, and we leave the study of the deep nonlinear phase to another dedicated article. In particular, we focus here on the nonlinear modification of the EP distribution function at the time of the saturation, and on the corresponding nonlinear modification of the EGAM frequency.

Figure 2.

Radial electric field in time, at different radial positions (a), and EP distribution function at different times, versus parallel velocity (b) for an EGAM simulation with ORB5 with $n_{EP}/n_{i}=0.10$ .

The EP distribution function redistributes in $v_{\Vert }$ during the first nonlinear phase, causing a relaxation of the drive due to the inverse Landau damping. The EP distribution function of this simulation is shown in figure 2 b). The redistribution of the EPs is observed to occur in a range of velocities between $2.5v_{ti}$ and $4.5v_{ti}$ . The EP distribution function does not change during the linear phase, and when entering the nonlinear phase, the redistribution occurs with higher-velocity EP moving towards lower values of $v_{\Vert }$ , as time increases. Therefore, negative values of the perturbed distribution function are measured at high velocities, and positive at low velocities. The resonance velocity can be calculated as $v_{\Vert \text{res}}=qR\unicode[STIX]{x1D714}_{\text{EGAM}}=3.5v_{ti}$ , and can be measured in figure 2(b) as the velocity where the perturbed distribution function changes sign. We note that this velocity measured in figure 2(b), for this value of $n_{EP}/n_{i}=0.10$ , does not sensibly change in the time range of interest.

Before moving further, we want to consider another case for comparison, with a stronger drive, namely with $n_{EP}/n_{i}=0.176$ . The evolution in time of the radial electric field is shown in figure 3(a). The linear frequency and growth rate is measured and checked to match with the ones of the linear simulation: $\unicode[STIX]{x1D714}_{L}=1.14\unicode[STIX]{x1D714}_{s}$ , $\unicode[STIX]{x1D6FE}_{L}=0.094\unicode[STIX]{x1D714}_{s}$ . The saturated level of the radial electric field is measured at $\unicode[STIX]{x1D6FF}E_{r}\simeq 0.8\times 10^{5}~\text{V}~\text{m}^{-1}$ (in agreement with the prediction of Biancalani et al. (Reference Biancalani, Chavdarovski, Qiu, Bottino, Del Sarto, Ghizzo, Gürcan, Morel and Novikau2017)).

Figure 3.

Radial electric field in time, at different radial positions (a), and EP distribution function at different times, versus parallel velocity (b) for an EGAM simulation with ORB5 with $n_{EP}/n_{i}=0.176$ .

Figure 4.

(a) Perturbed distribution function for a case with $n_{EP}/n_{i}=0.176$ . The times are expressed here in units of $\unicode[STIX]{x1D6FA}_{i}^{-1}$ . (b) Frequencies corresponding to the positive (red) and negative (blue) peaks of the EP perturbed distribution function. The measured EGAM frequency is also shown in green.

The EP distribution function of the case with $n_{EP}/n_{i}=0.176$ is shown in figure 3(b), for different times up to the saturation. Firstly, we note that the range of velocities affected by the nonlinear modification is broader than that for the weaker drive. In fact, the redistribution of the EPs is observed to occur in a range of velocities between $2v_{ti}$ and $5v_{ti}$ . Secondly, we note that the resonance frequency, which is calculated in this case from the linear frequency as $v_{\Vert \text{res}}=qR\unicode[STIX]{x1D714}_{\text{EGAM}}=3.2v_{ti}$ , does not perfectly describe the velocity of the change of sign of the perturbed distribution function at all times. In fact, the resonance velocity is observed to grow in time, from $3.2v_{ti}$ to $3.5v_{ti}$ .

The evolution of the resonance velocity in time is in relation to the EGAM nonlinear frequency modification, i.e. the EGAM chirping. The perturbed EP distribution function can be plotted explicitly (see figure 4 a). The positive peak (clump) and the negative peak (hole) can be seen to form and evolve in time, becoming bigger and centred at higher and higher distances from the linear resonance velocity. The location in velocity space of the peaks can be measured and translated into resonance frequencies as $\unicode[STIX]{x1D714}_{1,2}=v_{\Vert 1,2}/qR$ (see figure 4 b). When compared with the measured EGAM frequency, we note that the nonlinear EGAM frequency modification at the time of the saturation is described with a good approximation by the resonance frequency of the negative peak (hole). This relation offers the possibility to predict the nonlinear frequency by approximating it with the frequency obtained by the velocity of the negative peak of the EP perturbed distribution function (see also Wang et al. Reference Wang, Todo and Kim2013).

In the next section, § 5, we introduce the beam-plasma system (BPS) and the map linking the BPS with the EGAM. This map shows how we can predict the EGAM EP redistribution in velocity space.

5 Dynamics of the energetic particles: analogy with the beam-plasma system

In the EGAM system, the equilibrium magnetic field is not uniform but has a toroidal shape. In general, particles moving in a toroidal magnetic field can perform passing orbits, i.e. follow the magnetic field on both low-field side and high-field side, or perform banana orbits in the space restricted to the low-field side of the tokamak. By construction, in the EGAM problem considered here, the energetic ions are initialized with a bump-on-tail distribution function (beam distribution), with a relatively high parallel mean velocity along the equilibrium magnetic field, plus a smaller isotropic thermal distribution around the mean velocity. Due to the relatively high parallel velocity, the time derivative of their toroidal angle never vanishes (and therefore banana orbits of the EP on the low-field side are not considered in the present treatment). During their motion which is, to leading order, directed along the equilibrium magnetic field, they perform small drifts towards higher values of the minor radius, and then towards lower values of the minor radius, known as the curvature and grad-B drifts. These drifts have zero time average: this defines orbits with an average radial position plus a radial orbit width.

The radial electric field of the EGAM can exchange energy with the energetic particles, due to their radial component of the trajectories, and in particular of the curvature drift (Qiu et al. Reference Qiu, Zonca and Chen2011): $v_{dc}=(v_{\Vert }^{2}/\unicode[STIX]{x1D6FA}_{EP})\boldsymbol{B}\times \unicode[STIX]{x1D735}\boldsymbol{B}/B^{2}$ . Due to the structure of the perturbed distribution function, the effective parallel wavenumber of the EGAM is $k_{\Vert }=1/qR$ , the phase angle is $\unicode[STIX]{x1D6E9}=\unicode[STIX]{x1D703}-\unicode[STIX]{x1D714}_{\text{EGAM}}t$ and its normalized time derivative is $\dot{\unicode[STIX]{x1D6E9}}/\unicode[STIX]{x1D714}_{\text{EGAM}}=(v_{\Vert }-v_{\Vert r})/qR\unicode[STIX]{x1D714}_{\text{EGAM}}$ . In terms of the phase angle, the energetic particles experience a periodic electric field, and their harmonic motion can be expressed as (Qiu et al. Reference Qiu, Zonca and Chen2011):

(5.1) $$\begin{eqnarray}\displaystyle \frac{\ddot{\unicode[STIX]{x1D6E9}}}{\unicode[STIX]{x1D714}_{\text{EGAM}}^{2}}=-\frac{\unicode[STIX]{x1D714}_{b}^{2}}{\unicode[STIX]{x1D714}_{\text{EGAM}}^{2}}\sin \unicode[STIX]{x1D6E9}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D714}_{b}$ is the bounce frequency of the energetic particles in the potential created by the wave. Using these considerations, the analogy with respect to the 1-D beam-plasma system (BPS) turns out to be evident. In fact, the BPS is described as given by a Langmuir wave, excited by a beam of energetic electrons along a given direction $x$ . The Langmuir wave has a perturbed electric field directed along $x$ , oscillating at the plasma frequency $\unicode[STIX]{x1D714}_{p}=\sqrt{n_{e}e^{2}/m_{e}\unicode[STIX]{x1D716}_{0}}$ . In general, the Langmuir wave can be decomposed in Fourier in terms of the wavenumbers $k_{\ell }=\ell (2\unicode[STIX]{x03C0}/L)$ , where $L$ is the periodicity length of the 1-D space domain of the BPS, and $\ell$ is a positive integer (whereas the EGAM has only one possible wavenumber set by the equilibrium, as mentioned above). Considering a single monochromatic wave with a chosen value of $\ell$ , we can focus on that, and we denote the wavenumber as $k$ . The phase angle $\unicode[STIX]{x1D6E9}$ experienced by the electrons in the field of the Langmuir wave is thus $\unicode[STIX]{x1D6E9}=kx-\unicode[STIX]{x1D714}_{p}t$ and the normalized variation in time of the phase angle is $\dot{\unicode[STIX]{x1D6E9}}/\unicode[STIX]{x1D714}_{p}=k(v-v_{r})/\unicode[STIX]{x1D714}_{p}=kv/\unicode[STIX]{x1D714}_{p}-1$ , where the resonant velocity is defined as $v_{r}=\unicode[STIX]{x1D714}_{p}/k$ . As a choice of nomenclature, we refer here to the velocities of the EGAM as $v_{\Vert }$ , and to the velocities of the BPS as $v$ . Similarly, we refer to the wavenumbers of the EGAM as $k_{\Vert }=1/qR$ , and to the wavenumbers of the BPS as $k$ .

In the following, we describe the detailed mapping procedure which links the EGAM framework with the BPS. As already mentioned, the dynamics of the EGAM model can be reduced in the parallel velocity direction and we start from the generic resonance condition written using two suitable normalization constants $\unicode[STIX]{x1D708}_{1,2}$ , i.e.

(5.2) $$\begin{eqnarray}\displaystyle \frac{v_{\Vert }-v_{\Vert 0}}{\unicode[STIX]{x1D708}_{1}}=\frac{v-v_{r}}{\unicode[STIX]{x1D708}_{2}}, & & \displaystyle\end{eqnarray}$$

where we recall that the transit resonance velocity reads $v_{\Vert \text{res}}=qR\unicode[STIX]{x1D714}_{\text{EGAM}}(n_{EP})$ . Using the introduced above standard normalization $\unicode[STIX]{x1D708}_{1}=v_{ti}$ and, for the calculation of this section, $\bar{v}_{\Vert }=v_{\Vert }/v_{ti}$ , in the following, we denote with $\bar{v}_{\Vert \text{min}}\leqslant \bar{v}_{\Vert }\leqslant \bar{v}_{\Vert \text{max}}$ the domain of the positive bump of energetic particles. Imposing the boundary $\bar{v}_{\Vert \text{min}}=0\mapsto v_{\text{min}}=0$ , in order to map one single bump of energetic particles with $\bar{v}_{\Vert }\geqslant 0$ , we get $\unicode[STIX]{x1D708}_{2}=v_{r}/\bar{v}_{\Vert r}$ and the map finally writes

(5.3) $$\begin{eqnarray}\displaystyle v=\frac{v_{r}}{\bar{v}_{\Vert r}}\bar{v}_{\Vert }. & & \displaystyle\end{eqnarray}$$

Let us now introduce the following normalization: $v=\unicode[STIX]{x1D714}_{p}(2\unicode[STIX]{x03C0}/L)^{-1}u$ . In order to fix the dimensionless resonant wavenumber $\ell _{r}$ , we use the condition $k_{1}v_{\text{max}}=\unicode[STIX]{x1D714}_{p}$ , with $v_{\text{max}}=v_{r}\bar{v}_{\Vert \text{max}}/\bar{v}_{\Vert r}$ , which characterizes the spectral features (wavenumbers and periodicity length). This yields $\ell _{r}=\ell _{1}\bar{v}_{\Vert \text{max}}/\bar{v}_{\Vert r}$ , and $\ell _{r}$ is determined arbitrarily fixing $\ell _{1}$ since $\bar{v}_{\Vert \text{max}}$ and $\bar{v}_{\Vert r}$ are given quantities from the EGAM system. We stress how the resonance condition can be rewritten as $\ell _{r}u_{r}=1$ .

The map between the velocities of the two systems is now closed. The bump (positive part) of the EP is described by the shifted Maxwellian distribution function $F_{EP}(v_{\Vert })$ in velocity space (Biancalani et al. Reference Biancalani, Chavdarovski, Qiu, Bottino, Del Sarto, Ghizzo, Gürcan, Morel and Novikau2017). For modelling the EP distribution function of the EGAM in the BPS, let us now discretize the positive bump of $F_{EP}(v_{\Vert })$ in $n$ delta-like beams, equispaced in velocity space and located in $\bar{v}_{\Vert j}$ (with $j=1,\ldots ,n$ ), and assign the numbers of particles $N_{j}$ for each beam distributed according to $F_{EP}$ . The initial conditions on the distribution for BPS simulations are now given by $N_{j}$ particles located at

(5.4) $$\begin{eqnarray}\displaystyle u_{j}=\bar{v}_{\Vert j}/(\ell _{r}\bar{v}_{\Vert r}). & & \displaystyle\end{eqnarray}$$

For the sake of completeness, we mention that, for the simulation of the BPS, we have set $n=600$ , $\ell _{1}=400$ and we have used $N=10^{6}$ total particles. The complete derivation of the BPS dynamical equation used here, is described in Carlevaro et al. (Reference Carlevaro, Falessi, Montani and Zonca2015, Reference Carlevaro, Milovanov, Falessi, Montani, Terzani and Zonca2016) (and references therein), and it can be specified for one single resonance as:

(5.5a-c ) $$\begin{eqnarray}\displaystyle \bar{x}_{i}^{\prime }=u_{i},\quad u_{i}^{\prime }=\text{i}\ell _{r}\bar{\unicode[STIX]{x1D719}}_{r}\text{e}^{\text{i}\ell _{r}\bar{x}_{i}}+\text{c.c.},\quad \bar{\unicode[STIX]{x1D719}}_{r}^{\prime }=-\text{i}\bar{\unicode[STIX]{x1D719}}_{r}+\frac{\text{i}\unicode[STIX]{x1D702}}{2\ell _{r}^{2}N}\mathop{\sum }_{i=1}^{N}\text{e}^{-\text{i}\ell _{r}\bar{x}_{i}}, & & \displaystyle\end{eqnarray}$$

where the particle position along the $x$ direction is labelled by $x_{i}$ , with $i=1,\ldots ,N$ ( $N$ being the total particle number) normalized as $\bar{x}_{i}=x_{i}(2\unicode[STIX]{x03C0}/L)$ . The Langmuir electrostatic scalar potential $\unicode[STIX]{x1D719}(x,t)$ is expressed in terms of the Fourier component $\unicode[STIX]{x1D719}_{r}(k,t)$ and we have used: $\unicode[STIX]{x1D702}=n_{B}/n_{p}$ (for the plasma density $n_{p}$ assumed much greater than the beam one $n_{B}$ ), $\unicode[STIX]{x1D70F}=t\unicode[STIX]{x1D714}_{p}$ (the prime indicates derivative with respect to this variable), $\tilde{\unicode[STIX]{x1D719}}_{r}=(2\unicode[STIX]{x03C0}/L)^{2}e\unicode[STIX]{x1D719}_{r}/m\unicode[STIX]{x1D714}_{p}^{2}$ , $\bar{\unicode[STIX]{x1D719}}_{r}=\tilde{\unicode[STIX]{x1D719}}_{r}\text{e}^{-\text{i}\unicode[STIX]{x1D70F}}$ . Equations (5.5) are solved using a Runge–Kutta (fourth order) algorithm. For the considered time scales and for an integration step $h=0.1$ , both the total energy and momentum (for the explicit expressions, see Carlevaro et al. (Reference Carlevaro, Falessi, Montani and Zonca2015)) are conserved with relative fluctuations of approximately $1.4\times 10^{-5}$ .

The BPS is closed once the density of the beam (drive) is fixed. In order to quantitatively compare the nonlinear features of two systems, we now fix the bounce (trapping) frequency $\unicode[STIX]{x1D714}_{B}$ normalized to the mode frequency equal for the two schemes. For the BPS, the bounce frequency results are proportional to the linear growth rate of the mode $\unicode[STIX]{x1D6FE}_{L,\text{BPS}}$ . The same occurs for the EGAM system, but with a proportionality factor depending on the EP density (Biancalani et al. Reference Biancalani, Chavdarovski, Qiu, Bottino, Del Sarto, Ghizzo, Gürcan, Morel and Novikau2017). In particular, we get:

(5.6a,b ) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x1D714}_{B,\text{BPS}}}{\unicode[STIX]{x1D714}_{p}}=\unicode[STIX]{x1D6FC}\frac{\unicode[STIX]{x1D6FE}_{L,\text{BPS}}}{\unicode[STIX]{x1D714}_{p}},\quad \frac{\unicode[STIX]{x1D714}_{B,\text{EGAM}}}{\unicode[STIX]{x1D714}_{L,\text{EGAM}}}=\unicode[STIX]{x1D6FD}(\bar{n}_{EP})\frac{\unicode[STIX]{x1D6FE}_{L,\text{EGAM}}}{\unicode[STIX]{x1D714}_{L,\text{EGAM}}}. & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D6FC}\simeq 3.3$ (see well-known literature results in Chen & Zonca (Reference Chen and Zonca2016), Wu et al. (Reference Wu, White, Chen and Rosenbluth1995) and also Carlevaro et al. (Reference Carlevaro, Montani and Zonca2018)) while for the EGAM we have $\unicode[STIX]{x1D6FD}=\unicode[STIX]{x1D6FD}_{0}\sqrt{\unicode[STIX]{x1D714}_{L,\text{EGAM}}/\unicode[STIX]{x1D714}_{\text{GAM}}}$ , with $\unicode[STIX]{x1D6FD}_{0}=2.66$ in this regime (see § 3), and $\unicode[STIX]{x1D714}_{L,\text{EGAM}}$ depending on the EP density (Biancalani et al. Reference Biancalani, Chavdarovski, Qiu, Bottino, Del Sarto, Ghizzo, Gürcan, Morel and Novikau2017). For the four selected EGAM simulations, we get $\unicode[STIX]{x1D6FD}(0.07,0.10,0.176,0.30)\simeq [2.21,2.17,2.07,1.98]$ . Using standard normalization for frequencies, i.e. $\bar{\unicode[STIX]{x1D6FE}}_{L,\text{BPS}}=\unicode[STIX]{x1D6FE}_{L,\text{BPS}}/\unicode[STIX]{x1D714}_{p}$ and $\bar{\unicode[STIX]{x1D6FE}}_{L,\text{EGAM}}=\unicode[STIX]{x1D6FE}_{L,\text{EGAM}}/\unicode[STIX]{x1D714}_{\text{GAM}}$ , $\bar{\unicode[STIX]{x1D714}}_{L,\text{EGAM}}=\unicode[STIX]{x1D714}_{L,\text{EGAM}}/\unicode[STIX]{x1D714}_{\text{GAM}}$ and equalling the bounce frequencies, we finally get

(5.7) $$\begin{eqnarray}\displaystyle \bar{\unicode[STIX]{x1D6FE}}_{L,\text{BPS}}=\frac{\unicode[STIX]{x1D6FD}}{\unicode[STIX]{x1D6FC}}\frac{\bar{\unicode[STIX]{x1D6FE}}_{L,\text{EGAM}}}{\bar{\unicode[STIX]{x1D714}}_{L,\text{EGAM}}}. & & \displaystyle\end{eqnarray}$$

This condition preserves the linear and nonlinear features of the two systems and it is used in the evaluation for the drive of BPS simulations from the linear dispersion relation which formally reads as

(5.8) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D716}=1-\frac{\unicode[STIX]{x1D714}_{p}^{2}}{\unicode[STIX]{x1D714}^{2}}=\frac{\unicode[STIX]{x1D702}\unicode[STIX]{x1D714}_{p}^{2}}{k^{2}}\int _{-\infty }^{+\infty }\text{d}v\frac{k\unicode[STIX]{x2202}_{v}\hat{F}_{B}(v)}{kv-\unicode[STIX]{x1D714}}, & & \displaystyle\end{eqnarray}$$

where $\hat{F}_{B}(v)$ is the initial beam distribution function. Here, the dielectric function $\unicode[STIX]{x1D716}$ can be expanded near $\unicode[STIX]{x1D714}\simeq \unicode[STIX]{x1D714}_{p}$ to deal with Langmuir modes as in (5.5), i.e.  $\unicode[STIX]{x1D716}\simeq 2(\bar{\unicode[STIX]{x1D714}}-1)$ (where $\bar{\unicode[STIX]{x1D714}}=\unicode[STIX]{x1D714}/\unicode[STIX]{x1D714}_{p}$ ). Let us now use the expansion $\bar{\unicode[STIX]{x1D714}}=\bar{\unicode[STIX]{x1D714}}_{0}+\text{i}\bar{\unicode[STIX]{x1D6FE}}_{L,\text{BPS}}$ , where $\bar{\unicode[STIX]{x1D714}}_{0}$ is the real part of the normalized Langmuir frequency $\bar{\unicode[STIX]{x1D714}}$ . Using the linear character of the mapping which yields the normalization $\hat{F}_{B}(v)=\unicode[STIX]{x1D705}F_{EP}(\bar{v}_{\Vert })$ (with $\unicode[STIX]{x1D705}=\text{const}.$ ), (5.8) can be written in terms of the EGAM system variables as

(5.9) $$\begin{eqnarray}\displaystyle 2(\bar{\unicode[STIX]{x1D714}}_{0}+\text{i}\bar{\unicode[STIX]{x1D6FE}}_{L,\text{BPS}}-1)-\frac{\unicode[STIX]{x1D702}\bar{v}_{\Vert r}}{M}\int _{-\infty }^{+\infty }\text{d}\bar{v}_{\Vert }\frac{\unicode[STIX]{x2202}_{\bar{v}_{\Vert }}F_{EP}}{\bar{v}_{\Vert }/\bar{v}_{\Vert r}-\bar{\unicode[STIX]{x1D714}}_{0}-\text{i}\bar{\unicode[STIX]{x1D6FE}}_{L,\text{BPS}}}=0, & & \displaystyle\end{eqnarray}$$

where $M=\int _{-\infty }^{+\infty }\text{d}\bar{v}_{\Vert }F_{EP}$ . This equation is numerically integrated assuming (5.7), which guarantees the requested features described above, and provides the drive parameter $\unicode[STIX]{x1D702}$ closing the map procedure.

6 Nonlinear EP redistribution of the EGAM, and comparison with the beam-plasma model

In this section, we compare the results of the EP redistribution due to the EGAM, with the redistribution due to the beam-plasma system. In particular, we aim at predicting, from BPS information, the nonlinear parallel velocity spread in the positive bump of the distribution function.

In the BPS, the single mode dynamics proceeds in an initial exponential mode growth followed by nonlinear saturation. Here the particles get trapped and begin to bounce back and forth in the potential well generating clumps. A measure of the clumps width $\unicode[STIX]{x0394}u_{NL}^{c}$ for a generic initial half-Gaussian velocity distribution, which can be directly extrapolated to the analysis of the present work, has been evaluated in Carlevaro et al. (Reference Carlevaro, Montani and Zonca2018) for several cases outlining the following scaling rule as function of the linear drive:

(6.1) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x0394}u_{NL}^{c}/u_{r}=(6.64\pm 0.12)\bar{\unicode[STIX]{x1D6FE}}_{L}. & & \displaystyle\end{eqnarray}$$

In order to include the dynamic role also of passing but nearly resonant particles (Escande & Elskens Reference Escande and Elskens2002), i.e. the region involved in the effective wave–particle power exchange, in the following analysis we consider as the proper nonlinear particle velocity spread the scaled quantity $\unicode[STIX]{x0394}u_{NL}\simeq \unicode[STIX]{x1D712}\unicode[STIX]{x0394}u_{NL}^{c}$ with $\unicode[STIX]{x1D712}\simeq 1.28$ . This estimate is derived (Carlevaro et al. Reference Carlevaro, Montani and Zonca2018) characterizing the active overlap of different nonlinear fluctuations (Escande & Doveil Reference Escande and Doveil1981; Lichtenberg & Lieberman Reference Lichtenberg and Lieberman2010) and corresponds to the finite distortion of the distribution function, including effects at the edges of the plateau (defined as the flattened region of the distribution function, mainly coinciding with the clump size). We finally obtain the desired formula for the prediction of the nonlinear velocity spread due to the EGAM:

(6.2) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x0394}v_{\Vert NL}}{v_{\Vert L,\text{res}}}=8.5\bar{\unicode[STIX]{x1D6FE}}_{L}^{\text{BPS}} & & \displaystyle\end{eqnarray}$$

and by substituting the value of $\bar{\unicode[STIX]{x1D6FE}}_{L}^{\text{BPS}}$ we obtain:

(6.3) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x0394}v_{\Vert NL}}{v_{\Vert L,\text{res}}}=2.57\frac{\unicode[STIX]{x1D6FD}_{0}}{\sqrt{\unicode[STIX]{x1D714}_{\text{GAM}}}}\frac{\unicode[STIX]{x1D6FE}_{L}^{\text{EGAM}}}{\sqrt{\unicode[STIX]{x1D714}_{L,\text{EGAM}}}}. & & \displaystyle\end{eqnarray}$$

Equation (6.3) has been derived for the EGAM system, using the scaling derived in (5.7) and the normalized mapping of (5.4). For the regime of interest in this article, we have $\unicode[STIX]{x1D6FD}_{0}=2.66$  (Biancalani et al. Reference Biancalani, Chavdarovski, Qiu, Bottino, Del Sarto, Ghizzo, Gürcan, Morel and Novikau2017), and $\unicode[STIX]{x1D714}_{\text{GAM}}=1.8\unicode[STIX]{x1D714}_{s}$ , therefore we obtain:

(6.4) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x0394}v_{\Vert NL}}{v_{\Vert L,\text{res}}}=\frac{5.1}{\sqrt{\unicode[STIX]{x1D714}_{s}}}\frac{\unicode[STIX]{x1D6FE}_{L}^{\text{EGAM}}}{\sqrt{\unicode[STIX]{x1D714}_{L,\text{EGAM}}}}. & & \displaystyle\end{eqnarray}$$

Figure 5.

Energetic-particle distribution function averaged in space, and measured at $\unicode[STIX]{x1D707}\simeq 0$ , versus parallel velocity, for $n_{EP}/n_{i}=0.07$ (a) and $n_{EP}/n_{i}=0.10$ (b). The vertical dashed and continuous lines are the resonance velocity (centre), with the borders of the nonlinear velocity predicted by $v_{\Vert \text{res}}\pm \unicode[STIX]{x0394}v_{\Vert NL}$ .

Figure 6.

Same as figure 5 but for $n_{EP}/n_{i}=0.176$ (a) and $n_{EP}/n_{i}=0.30$ (b).

Let us now analyse the predictivity of (6.3), simplified as in (6.4) for the regime of interest. Four different simulations are considered, with different values of the energetic-particle concentration: $n_{EP}/n_{i}\in [0.07,0.10,0.176,0.30]$ . The corresponding linear frequencies and growth rates are $\unicode[STIX]{x1D714}_{L,\text{EGAM}}=[1.30,1.24,1.14,1.04]\unicode[STIX]{x1D714}_{s}$ and $\unicode[STIX]{x1D6FE}_{L,\text{EGAM}}=[0.04,0.06,0.09,0.11]\unicode[STIX]{x1D714}_{s}$ . In all simulations, the distribution function is averaged in space, and measured at $\unicode[STIX]{x1D707}\simeq 0$ . Snapshots of the distribution function at different times are selected, going from the linear phase to the instant when the first peak of the electric field is reached, i.e. at saturation. In figures 5 and 6, the distribution functions measured with ORB5 are shown to be modified by the nonlinear interaction with the EGAM, with a certain width in the velocity space. The position of the linear resonance velocity and the nonlinearly modified resonance velocity are shown for each simulation respectively as a dashed vertical black line, and a continuous vertical black line. We can now compare with the predicted values of (6.4): $\unicode[STIX]{x0394}v_{\Vert NL}/v_{\Vert L,\text{res}}=[0.17,0.28,0.45,0.55]$ . The corresponding predicted range is delimited by vertical red lines in figures 5 and 6. The dashed lines correspond to the range calculated with respect to the linear resonance velocity, and the continuous lines with respect to the nonlinearly modified resonance velocity. Note that the predicted width of the velocity domain of EP redistribution, centred at the nonlinearly modified resonance velocity, fits very well with the results of ORB5. Note also that, for the two cases with lowest drive, the nonlinear modification of the resonance velocity is negligible, and therefore the predicted width of the velocity domain of EP redistribution fits very well the results of ORB5, even when centred at the linear resonance velocity.

Figure 7.

Plot of the BPS distribution function taken at saturation time (dashed-black line) mapped back in the $v_{\Vert }$ space, over the evolution of the EP profile (coloured lines in standard colour scheme from blue ( $t=2000\unicode[STIX]{x1D6FA}_{i}^{-1}$ ) to red ( $t=20\,000\unicode[STIX]{x1D6FA}_{i}^{-1}$ )) for $n_{H}/n_{i}=0.17$ (a) and $n_{H}/n_{i}=0.3$ (b).

In summary, two distinct regimes can be identified. When the instabilities are weakly driven (figure 5), a very good match between the estimated deviation $v_{\Vert \text{res}}\pm \unicode[STIX]{x0394}v_{\Vert NL}$ and the nonlinear EP redistribution is observed, where $v_{\Vert \text{res}}$ is the linear resonance velocity. Otherwise, in the strongly driven regime (figure 6) the importance of the frequency chirping comes out. In particular, nonetheless the value $v_{\Vert \text{res}}\pm \unicode[STIX]{x0394}v_{\Vert NL}$ remains very predictive, as long as the nonlinearly modified resonance velocity is chosen as centre of the EP redistribution due the EGAM.

The nonlinear frequency shift which is characteristic of the EGAM, can not be intrinsically implemented in the BPS model. In fact, differences remain in the physics of the EGAM and the BPS. This clearly emerge in figure 7, where the distribution function of the BPS at saturation (dashed-black line), mapped back in the $v_{\Vert }$ space, is overplotted on the evolution of the EP profile for $n_{EP}/n_{i}=0.176$ and $n_{EP}/n_{i}=0.3$ . In particular, it is evident how the discrepancy due the fixed character (at ${\sim}\unicode[STIX]{x1D714}_{p}$ ) of the Langmuir resonance gives rise to a very different morphology of the distribution function, although it predicts well the effective nonlinear velocity spread. Also the inclusion of additional modes with artificial ad hoc damping rates results in a drastically non-comparable nonlinear dynamics, underlining the intrinsic differences of the physical systems.

7 Conclusions and discussion

Geodesic acoustic modes (GAMs), i.e. finite-frequency zonal (i.e. axisymmetric) flows with mainly radial electric field polarization, are known to be important in tokamaks due to their interaction with turbulence. GAMs can also be excited by energetic particles (EPs), taking the name of EGAMs. In the view of understanding and predicting the EP redistribution in tokamaks, the linear and nonlinear dynamics of EGAMs should be properly theoretically understood.

In this article, we have investigated the nonlinear dynamics of EGAMs with particular interest in the EP redistribution in velocity space. As a tool of investigation of the nonlinear inverse Landau damping, which is responsible of the EGAM saturation and EP redistribution, a comparison of the EGAM dynamics with the beam-plasma system (BPS) is done. Although the BPS describes the interaction of EP with a different mode, i.e. the Langmuir wave in a 1-D geometry, nevertheless, analogies with the EGAM had been suggested in previous articles. These analogies are investigated here and used to build a mapping of the two systems. The mapping is used to understand and predict the EP redistribution of the EGAM.

The EGAM is investigated here with collisionless electrostatic simulations with the gyrokinetic (GK) particle-in-cell code ORB5. Only the wave–particle nonlinearity is retained in the simulations presented here, meaning that the markers for thermal ions follow unperturbed trajectories in phase space. The GK model allows us to study the EGAM dynamics, retaining the crucial physics of the resonances with thermal and fast ions. The BPS is investigated with the a 1-D code treating the thermal plasma as a cold dielectric medium and describing the dynamics of the fast particles (fast electrons in the case of the BPS) as an N-body problem solved with an Hamiltonian formulation.

The GK simulations of the EGAM show that, after a first linear growth, the EGAM enters a first nonlinear phase where the EP distribution function suffers a modification due to the EGAM field. A saturation of the EGAM field occurs, and then a deep nonlinear phase is entered, where nonlinear oscillations of the fields are observed. In this article, we are interested in the first nonlinear phase only, up to the first saturation. The EP population is observed to redistribute, with EP going from higher to lower values of the parallel component of the velocity, as the EGAM grows in amplitude. The resonant velocity is also measured with ORB5 in the different cases considered.

The mapping of the EGAM and the BPS is then described, and the comparison with the EP redistribution in velocity space is shown. In particular, the main result is the prediction of the width of the velocity space $\unicode[STIX]{x0394}u_{NL}$ which is affected by the EGAM, around the resonance. The implication of this result is evident, as reduced models are needed for predicting the nonlinear dynamics of instabilities in tokamaks, instead of using numerically expansive GK simulations. In fact, BPS simulations are numerically much cheaper, and therefore the mapping described here offers a tool for predicting the EP redistribution in regimes where EGAMs experience wave–particle nonlinear saturation.

We have also found a transition among two regimes: for weakly driven EGAMs, the resonant velocity does not evolve in time during the first nonlinear phase, exactly like in the BPS problem; on the other hand, above a certain threshold in drive, the resonant velocity slightly increases in time, and therefore a difference of the EGAM and BPS is found. The increase of the resonant velocity in the high-drive regime is consistent with the nonlinear frequency modification (i.e. the frequency chirping) which is present for the EGAM, and is absent for the BPS considered here. The difference can be cured by calculating the nonlinear width $\unicode[STIX]{x0394}u_{NL}$ around the new value of the resonant velocity. This proves that the one-to-one correspondence of the EP redistribution around the resonance is completely described by the nonlinear inverse Landau damping which is included in the model of the BPS. Therefore, we can state that the EP are redistributed around the resonance by the EGAM for a purely 1-D problem which is the same as the BPS, namely the nonlinear inverse Landau damping. Having tried to include many modes in the BPS, we have observed that this can modify the EP redistribution, but leads to a overestimation of the EP redistribution due to a lack of damping in the different modes.

In the view of the continuation of this work, several next steps can be taken in the direction of getting closer to more and more realistic scenarios. For example, the inclusion of wave–wave coupling of the EGAM with itself is under investigation. Kinetic electron effects should also be considered, as the electron resonances might affect the Landau damping. The inclusion of turbulence is also in progress, as a mean of modifying the EGAM saturation level and EP redistribution.