2. The model

The numerical tool used here is the gyrokinetic particle-in-cell code ORB5, originally written for electrostatic turbulence studies (Jolliet et al. Reference Jolliet, Bottino, Angelino, Hatzky, Tran, Mcmillan, Sauter, Appert, Idomura and Villard2007) and then extended to its electromagnetic multispecies version (Bottino et al. Reference Bottino, Vernay, Scott, Brunner, Hatzky, Jolliet, McMillan, Tran and Villard2011; Mishchenko et al. Reference Mishchenko, Bottino, Biancalani, Hatzky, Hayward-Schneider, Ohana, Lanti, Brunner, Villard and Borchardt2019; Lanti et al. Reference Lanti, Ohana, Tronko, Hayward-Schneider, Bottino, McMillan, Mishchenko, Scheinberg, Biancalani and Angelino2020). ORB5 is global, i.e. it resolves modes with structure comparable with the minor radius. Thus, it is appropriate for studying AMs with low toroidal mode number.

In this paper, we use two versions of the code. In the first part of the numerical experiment, we run self-consistent nonlinear simulations of AMs driven by EPs in the presence of turbulence, therefore the electromagnetic version of the code is used, as done by Biancalani et al. (Reference Biancalani, Bottino, Brunner, Di Siena, Görler, Hatzky, Jenko, Könies, Lauber and Novikau2019, Reference Biancalani, Bottino, Di Siena, Gürcan, Hayward-Schneider, Jenko, Lauber, Mishchenko, Morel and Novikau2021). In the second part of the numerical experiment, we want to study the isolated effect of the nonlinearly modified profiles on ITG turbulence, therefore we use the electrostatic version of the code to avoid the development of AMs.

In the electromagnetic version of the code, the magnetic potential is split into the Hamiltonian and symplectic parts $A_{\|} = A_{\|}^{{(h)}} + A_{\|}^{{(s)}}$ (see Mishchenko et al. Reference Mishchenko, Könies, Kleiber and Cole2014 for details). The perturbed equations of motion in mixed-variable formulation are (Mishchenko et al. Reference Mishchenko, Bottino, Biancalani, Hatzky, Hayward-Schneider, Ohana, Lanti, Brunner, Villard and Borchardt2019)

(2.1)\begin{gather} \dot{{\boldsymbol{R}}}^{(1)} = \frac{{\boldsymbol{b}}}{B_{\|}^*} \times \boldsymbol{\nabla} \left\langle \phi - v_{\|} A_{\|}^{(s)} - v_{\|} A_{\|}^{(h)} \right\rangle - \frac{q}{m} \langle A^{(h)}_{\|} \rangle {\boldsymbol{b}}^*, \end{gather}

(2.2)\begin{gather}\dot{v}_{\|}^{(1)} ={-} \frac{q}{m} \left[ {\boldsymbol{b}}^* \boldsymbol{\cdot} \boldsymbol{\nabla} \left\langle \phi - v_{\|} A_{\|}^{(h)} \right\rangle + \frac{\partial{}}{\partial{t}} \left\langle A_{\|}^{(s)} \right\rangle \right] - \frac{\mu}{m} \frac{{\boldsymbol{b}} \times \boldsymbol{\nabla} B}{B_{\|}^*} \boldsymbol{\cdot} \boldsymbol{\nabla} \left\langle A_{\|}^{(s)} \right\rangle, \end{gather}

where

(2.3)\begin{gather} {\boldsymbol{b}}^* = {\boldsymbol{b}}^*_0 + \frac{\boldsymbol{\nabla}\langle A_{\|}^{{(s)}} \rangle\times {\boldsymbol{b}}}{B_{\|}^*} , \quad {\boldsymbol{b}}^*_0 = {\boldsymbol{b}} + \frac{m v_{\|}}{q B_{\|}^*} \boldsymbol{\nabla}\times{\boldsymbol{b}}, \end{gather}

(2.4)\begin{gather}B_{\|}^* = B + \frac{m v_{\|}}{q} {\boldsymbol{b}}\boldsymbol{\cdot}\boldsymbol{\nabla}\times{\boldsymbol{b}}. \end{gather}

We also have an equation for $\partial A_{\|}^{(s)} / \partial t$ (Ohm's law, see Mishchenko et al. Reference Mishchenko, Könies, Kleiber and Cole2014):

(2.5)\begin{equation} \frac{\partial{}}{\partial{t}}A_{\|}^{(s)} + {\boldsymbol{b}} \boldsymbol{\cdot} \boldsymbol{\nabla} \phi = 0 \end{equation}

and the field equations (here we use the notation as in Mishchenko et al. Reference Mishchenko, Könies, Kleiber and Cole2014):

(2.6)\begin{gather} -\boldsymbol{\nabla}\boldsymbol{\cdot}\left(\frac{n_0}{B \omega_{ci}}\boldsymbol{\nabla}_{{\perp}}\phi\right) = \bar{n}_{1i} - \bar{n}_{1e}, \end{gather}

(2.7)\begin{gather}\sum_{s = i,e}\frac{\beta_s}{\rho_{s}^{2}} A_{\|}^{(h)} - \boldsymbol{\nabla}_{{\perp}}^{2} A_{\|}^{(h)} = \mu_{0} \sum_{s = i,e} \bar{j}_{\|1s} + \boldsymbol{\nabla}_{{\perp}}^2 A_{\|}^{(s)}. \end{gather}

Note that this splitting of $A_\parallel$ into $A_\parallel ^{(h)}$ and $A_\parallel ^{(s)}$ is arbitrary and does not impose additional constraints on the dynamics, but it leads to improved numerical stability when $A_\parallel ^{(h)} \ll A_\parallel ^{(s)}$, which is typically true for Alfvénic physics (see Mishchenko et al. Reference Mishchenko, Könies, Kleiber and Cole2014). It has to be noted that (2.5) represents a choice used to separate the symplectic from the Hamiltonian part of $A_\|$. This choice is in principle arbitrary, but we have found that, in practice, this provides a numerically sane scheme in most situations. Regarding the energy conservation, this is guaranteed by the derivation of the pullback mitigation scheme using the Lie-transform Lagrangian approach (see Mishchenko et al. Reference Mishchenko, Könies, Kleiber and Cole2014).

For the electrostatic simulations, only the thermal ions and energetic ions (the EP species) are treated kinetically, whereas the electrons are treated adiabatically. For these simulations, only the scalar potential is needed, so the gyrokinetic Poisson law, (2.6), is solved.

For noise control purposes, as well as for maintaining some of the plasma profiles close to their initial state, a modified Krook operator is used for the thermal ions and electrons (not for the EPs):

(2.8)\begin{equation} \frac{\mathrm{d} f_s}{\mathrm{d} t} = S(f_s), \end{equation}

with $S(f_s)=-\gamma _K(f_s-f_{0,s}) + S_\mathrm {corr}(f_s)$. The coefficient $\gamma _K$ is chosen empirically such that the signal-to-noise ratio is maintained at a sufficiently high level, while only weakly affecting the physics of interest. It is typically chosen as 5 % or 10 % of the maximal growth rate. Here, we use the same value as used by Biancalani et al. (Reference Biancalani, Bottino, Di Siena, Gürcan, Hayward-Schneider, Jenko, Lauber, Mishchenko, Morel and Novikau2021), because the dynamics is very similar. The term $S_\mathrm {corr}(f_s)$ is such that a number of moments $M(v)$ are conserved by the source term: $\langle \int S(f_s)M(v)d^3v\rangle =0$, where $\langle Q \rangle$ stands for the flux-surface average of a quantity $Q$. Here, $S_\mathrm {corr}$ is used to conserve the undamped ExB Zonal Flow residual. For more details, see McMillan et al. (Reference McMillan, Jolliet, Tran, Villard, Bottino and Angelino2008) and Lanti et al. (Reference Lanti, Ohana, Tronko, Hayward-Schneider, Bottino, McMillan, Mishchenko, Scheinberg, Biancalani and Angelino2020).

3. Magnetic equilibrium and plasma profiles

For continuity with the previous work, we choose the magnetic equilibrium and plasma profiles very similar to those used in the case labelled here as the ‘EPS-2019 case’, published by Biancalani et al. (Reference Biancalani, Bottino, Brunner, Di Siena, Görler, Hatzky, Jenko, Könies, Lauber and Novikau2019, Reference Biancalani, Bottino, Lauber, Mishchenko and Vannini2020, Reference Biancalani, Bottino, Di Siena, Gürcan, Hayward-Schneider, Jenko, Lauber, Mishchenko, Morel and Novikau2021). The only two differences are the shape of the q-profile and the localisation of the density and temperature gradients, as described below.

Like in the EPS-2019 case, a magnetic equilibrium with inverse aspect ratio $\epsilon =0.1$ is considered (the major radius is $R_0 = 10$ m, the minor radius is $a=1.0$ m), with circular concentric flux surfaces. The magnetic field on axis is $B_0 = 3.0$ T. Differently from the EPS-2019 case, the q-profile is nearly monotonic here (see figure 1), with a rational surface at mid-radius, allowing a better localisation of the modes of interest. In particular, we have a value of $q(0)=1.79$ at the axis, a minimum of $q_{{\rm min}}=1.787$ at $s=0.33$, a value of $q=1.8$ at $s = 0.525$ and a value at the edge of $q(1) = 2.53$. Here, the flux radial coordinate $s$ is defined as $s=\sqrt {\psi _{{\rm pol}}/\psi _{{\rm pol}}({\rm edge})}$. The rational surface at $s=0.525$ corresponds to a normalised radius chosen as reference position, with value $\rho _r=0.5$. Here $\rho$ is a normalised radial coordinate defined as $\rho =r/a$.

Figure 1. Safety factor profile in $s$ coordinate.

In tokamaks, the type of turbulence under investigation strongly depends on the dimensionless parameter $\rho ^* = \rho _s/a$ (with $\rho _s = \sqrt {T_e/m_i}/\varOmega _i$ being the sound Larmor radius). Here, like in the EPS-2019 case, we choose a value of $\rho ^*$ similar to the CYCLONE base case (originally chosen as an international benchmark case for ITG turbulence in a DIII-D configuration): $\rho ^* = \rho _s/a = 0.00571$ (therefore, $Lx=2/\rho ^* = 350$). For comparison, note that $\rho ^* = 1/100$ in Ishizawa et al. (Reference Ishizawa, Imadera, Nakamura and Kishimoto2021). The electron thermal to magnetic pressure ratio of $\beta _e = 8{\rm \pi} \langle n_e \rangle T_e(\rho _r)/B_0^2 = 5\times 10^{-4}$ (with $\langle n_e \rangle$ being the volume averaged electron density). The equilibrium density and temperature profiles are different in the two parts of this study, and they are defined in the following way.

1. (i) In the first part, we use very similar equilibrium density and temperature profiles as in the EPS-2019 case. The initial profiles are described by the following equations:

   (3.1)\begin{gather} n(\rho)/n(\rho_r) = \exp [-\varDelta \kappa_n \tanh ((\rho-\rho_r)/\varDelta)], \end{gather}
   (3.2)\begin{gather}T(\rho)/T(\rho_r) = \exp [-\varDelta \kappa_t \tanh ((\rho-\rho_r)/\varDelta)], \end{gather}
   where both thermal ions and electrons have $\kappa _n = 0.3$ and $\kappa _t = 1.0$ respectively for the density and the temperature profiles. The only difference here, with respect to the EPS-2019 case, is that the density and temperature gradients are slightly more localised around mid-radius. This is done by selecting, for both thermal ions and electrons, a slightly smaller value of $\varDelta$: $\varDelta = 0.15$. Regarding the EP species, we select $\kappa _n = 20.0$ and $\kappa _t = 0.0$, $\varDelta = 0.10$, and $T_{{\rm EP}}/T_e=10$ at the reference position. In this first part of the study, we run nonlinear electromagnetic simulations of turbulence, zonal flows, AMs and EPs, similarly to the simulations shown by Biancalani et al. (Reference Biancalani, Bottino, Brunner, Di Siena, Görler, Hatzky, Jenko, Könies, Lauber and Novikau2019, Reference Biancalani, Bottino, Di Siena, Gürcan, Hayward-Schneider, Jenko, Lauber, Mishchenko, Morel and Novikau2021) for the EPS-2019 case.

2. (ii) In the second part of this study, the initial density and temperature profiles are not given by (3.1). However, we take the temperature and density profiles given in output from ORB5 in the electromagnetic simulations performed in the first part. These profiles are initialised in electrostatic ITG-turbulence simulations (without AMs and EPs).

4. First part of the study: EM simulation

In this section, we show the result of the nonlinear self-consistent EM simulation. The dynamics is very similar to that shown by Biancalani et al. (Reference Biancalani, Bottino, Brunner, Di Siena, Görler, Hatzky, Jenko, Könies, Lauber and Novikau2019, Reference Biancalani, Bottino, Di Siena, Gürcan, Hayward-Schneider, Jenko, Lauber, Mishchenko, Morel and Novikau2021). The slightly different profiles are found to affect the radial localisation of the AM. The evolution of the fields in time can be observed in figure 2(a). EPs are switched on at $t=6\times 10^4 \varOmega _i^{-1}$. The AM is a beta-induced Alfvén eigenmode (BAE) with $n=5$, $m=9$ (see figure 2b).

Figure 2. (a) Evolution of the radial electric field in time for a simulation with EPs (blue) and without EPs (red). In the simulation with EPs, an AM is driven unstable. (b) Structure of the AM in the poloidal plane at $t=73\,000$ $\varOmega _i^{-1}$. The AM is identified as a BAE.

Like in the case of Biancalani et al. (Reference Biancalani, Bottino, Brunner, Di Siena, Görler, Hatzky, Jenko, Könies, Lauber and Novikau2019, Reference Biancalani, Bottino, Di Siena, Gürcan, Hayward-Schneider, Jenko, Lauber, Mishchenko, Morel and Novikau2021), the BAE carries heat fluxes. This can be seen in figure 3(a). These heat fluxes modify the equilibrium profiles. As an example, the temperature profile of the thermal ions is shown at $t=0$ and $t=75\,000$ $\varOmega _i^{-1}$ in figure 3(b). Note that the temperature is flattened at the radial position of the BAE, i.e. around $\rho =0.4$. Note also that, like shown by Biancalani et al. (Reference Biancalani, Bottino, Di Siena, Gürcan, Hayward-Schneider, Jenko, Lauber, Mishchenko, Morel and Novikau2021), the direct interaction of EPs and turbulence is negligible in this regime. For completeness, note that the EP transport due to turbulence has been proven to be consistent with quasilinear theory, and to become negligible for large values of the EP temperature (White & Mynick Reference White and Mynick1989; Zhang et al. Reference Zhang, Decyk, Holod, Xiao, Lin and Chen2010; Chen & Zonca Reference Chen and Zonca2016).

Figure 3. (a) Evolution of the heat fluxes in time for a nonlinear electro-magnetic simulation with EPs. (b) Result of the effect of the heat fluxes on the modification of the temperature profiles of deuterium (similar effects are found on the electrons). The same profile measured at $t=75\,000$ $\varOmega _i^{-1}$ is used for the electrostatic simulations in the second part of this study.

We now take the temperature and density profiles measured at the end of this nonlinear simulation, and we use them as starting conditions for the electrostatic turbulence simulations discussed in the next section.

5. Second part of the study: (A) linear dynamics of ITG modes

In this section, we show the results of linear electrostatic simulations (therefore without AMs). We also keep only two species in the dynamics: thermal ion and electrons. No EPs are initialised. Two cases are compared: in one case, we load the original profiles (i.e. as in § 4); in the other case, we take the profiles of the EM simulation, modified by the presence of the AM, and we use these for the ES turbulence simulation. The goal is to see how the dynamics of the ITG modes differs in the two cases.

These linear simulations show that both the structure and the growth rate of the ITGs are different in the two kinds of profiles. The structure is shown in figure 4. The maximum linear growth rate for the case with the original (i.e unperturbed) profiles is measured as $\gamma _{{\rm lin}}\simeq 1.2\times 10^{-4} \varOmega _i$. The maximum linear growth rate for the case with the profiles modified by the AM mode is measured as $\gamma _{{\rm lin}}\simeq 1.0\times 10^{-4} \varOmega _i$. Therefore, ITGs are found to be linearly mitigated by the presence of the AM.

Figure 4. Structure of the ITG dominant mode in linear simulations (a) with original unperturbed profiles and (b) with the profiles modified by the AM.

6. Second part of the study: (B) Nonlinear dynamics of ITG turbulence

In this section, we show the results of nonlinear electrostatic simulations. Exactly the same case as shown in § 5 is considered. We want to compare the dynamics with unperturbed profiles and with profiles modified by the heat flux carried by the AM. A possible way to give an estimation of the turbulence intensity is by measuring the heat flux. In figure 5(a), the time evolution of the ion heat flux of the nonlinear ITG simulations is shown for both cases. Note that the heat flux in the simulation with modified profiles is approximately a factor of two lower than the heat flux of the simulation with unperturbed profiles. In figure 5(b), we can see that the time averaged radial heat flux is nearly suppressed at the radial position of the AM, where the profiles are flattened. In summary, we can observe an indirect mechanism of turbulence reduction of the EPs (by means of the AM).

Figure 5. Turbulence intensity measured with the ion heat fluxes. The two cases with original (blue) and with the modified (red) profiles are depicted. (a) Time evolution of the heat flux. Note that the heat flux is lower when the modified profiles are used. (b) Radial structure of the heat flux, averaged in time.

7. Application to an experimental case: the effect of AMs on the equilibrium profiles in the NLED-AUG case

In the previous sections, we have selected a tokamak configuration and we have performed self-consistent electromagnetic simulations, measured the nonlinearly modified temperature profile, and used these profiles to study the effect on ITG linear and nonlinear dynamics. For continuity with previous work, the chosen case is a tokamak configuration with low inverse aspect ratio, circular concentric flux surfaces and relatively low beta. This case allows self-consistent electromagnetic simulation including multiple scales, at a relatively low numerical cost.

In this section, we want to draw some conclusions on more experimentally relevant scenarios. We take the NLED-AUG case (Lauber Reference Lauber2014) as an example, because the linear and nonlinear physics of AMs has been studied in detail for this case with ORB5 (Vannini et al. Reference Vannini, Biancalani, Bottino, Hayward-Schneider, Lauber, Mishchenko, Novikau and Poli2020, Reference Vannini, Biancalani, Bottino, Hayward-Schneider, Lauber, Mishchenko, Poli and Vlad2021, Reference Vannini, Biancalani, Bottino, Hayward-Schneider, Lauber, Mishchenko, Poli, Rettino, Vlad and Wang2022b; Rettino et al. Reference Rettino, Hayward-Schneider, Biancalani, Bottino, Lauber, Chavdarovski, Weiland, Vannini and Jenko2022) and benchmarked with other codes (Vlad et al. Reference Vlad, Wang, Vannini, Briguglio, Carlevaro, Falessi, Fogaccia, Fusco, Zonca and Biancalani2021). Therefore, we can claim that our simulations correctly include the main nonlinear dynamics of these AMs. This preparatory phase is crucial to be able to make statements on the physics of Alfvén modes in experimentally relevant predictions.

The NLED-AUG simulations presented here use an experimental (shaped) magnetic equilibrium, and experimental density and temperature profiles of all the species. The equilibrium distribution function of the thermal species is a Maxwellian, and that of the EPs is an isotropic slowing-down. The isotropic slowing-down is a good approximation of the distribution function of particles entering the plasma with a certain energy (like alpha particles) and interacting with the other particles due to Coulomb collisions. In this case, the injection energy is $\mathcal {E} = 93$ keV and the EPs have a density concentration of $\langle n_{{\rm EP}} \rangle / \langle n_e \rangle = 0.0949$, where $\langle \cdots \!\rangle$ indicates volume average. The reader should refer to Vannini et al. (Reference Vannini, Biancalani, Bottino, Hayward-Schneider, Lauber, Mishchenko, Poli, Rettino, Vlad and Wang2022a) for more details on this case. Note that for the physics of interest in this paper, there is no qualitative dependence of the results on the initial conditions of the energetic particles (as long as they drive the AM unstable). In fact, the energetic particles are used here to drive the AM unstable and then the AM modifies the plasma profiles. It is only due to the plasma profile modification that the turbulence dynamics is modified in our ES simulations (no EPs are present there). In this simulation, modes among $n=0$ and $n=6$ are kept. In the linear phase, the most unstable mode is the mode with $n=2$. The simulation is fully nonlinear, meaning that the markers of all species are pushed along their full (equilibrium + perturbation) trajectories.

The radial structure of the mode, as given by Vannini et al. (Reference Vannini, Biancalani, Bottino, Hayward-Schneider, Lauber, Mishchenko, Poli, Rettino, Vlad and Wang2022a), shows a mode peaked in the core, at $s=\sqrt {\psi /\psi _{{\rm edge}}}\simeq 0.25$ (see Figure 6). The corresponding modification of the temperature profiles of thermal ions and electrons can be studied. In figure 7, the temperature profiles normalised with $T_e({\rm ref})$, at the beginning and at the end of the simulations (with $T_e({\rm ref})$ being the electron temperature measured at the axis, $\rho =0$, at the beginning of the simulation, $t=0$). Note that a sensible flattening of the temperature profiles is caused by the AM for both ions and electrons. This is especially visible near the AM radial localisation, i.e. $s \simeq 0.25$. As a consequence, we indicate that the mechanism of indirect interaction of EPs and ITG turbulence shown in this paper can be important not only in simplified configurations, but also in experimentally relevant regimes. We leave the analysis of this AUG shot and the comparison with experimental measurements to a dedicated work.

Figure 6. Radial structure of the AM in the NLED-AUG case.

Figure 7. Radial profiles of the temperature of (a) thermal ions and (b) electrons at the initial state of the nonlinear simulations of AMs and after the nonlinear saturation.