Gyrokinetic investigation of the damping channels of Alfvén modes in ASDEX Upgrade

Vannini, F.; Biancalani, A.; Bottino, A.; Hayward-Schneider, T.; Lauber, P.; Mishchenko, A.; Novikau, I.; Poli, E.; ASDEX Upgrade Team

doi:10.1063/1.5134979

The linear destabilization and nonlinear saturation of energetic-particle driven Alfvénic instabilities in tokamaks strongly depend on the damping channels. In this work, the collisionless damping mechanisms of Alfvénic modes are investigated within a gyrokinetic framework by means of global simulations with the particle-in-cell code ORB5 and compared with the eigenvalue code LIGKA and reduced models. In particular, the continuum damping and the Landau damping (of ions and electrons) are considered. The electron Landau damping is found to be dominant compared to the ion Landau damping for experimentally relevant cases. As an application, the linear and nonlinear dynamics of toroidicity induced Alfvén eigenmodes and energetic-particle driven modes in ASDEX Upgrade is investigated theoretically and compared with experimental measurements.

I. INTRODUCTION

In burning plasmas relevant for magnetic fusion energy research, an important role is played by energetic particles (EPs). With the term EPs, we refer to fusion reaction products (alpha particles) or superthermal ions or electrons resulting from plasma heating. Such particles possess higher velocities compared to those typical of the background plasma. In a typical tokamak, the EP velocity is of the order of the Alfvén speed $v_{A} = \sqrt{B / (4 π ρ_{m 0})}$ (with $ρ_{m, 0}$ the background plasma mass density and B the background magnetic field strength). Therefore, they can excite, through resonant wave-particle interactions, waves whose frequency falls in the magnetohydrodynamic (MHD) domain [ $O (10^{- 2} Ω_{c i})$ ]. Among them, the most relevant are the nearly incompressible, anisotropic, shear Alfvén waves (SAWs). These are instabilities with a group velocity parallel to the background magnetic field and localized at the surface where $ω = k_{| |} v_{A}$ ⁠, with $k_{∥}$ the component of the wave vector parallel to the background magnetic field (⁠ $k_{∥} = k \cdot B / B$ ⁠) and ω the frequency of the wave. This situation, strictly valid in slab geometry, is modified in tokamak devices where discrete Alfvén eigenmodes (AEs) can exist in the frequency gaps of the SAW continuum spectrum, creating a zoology of modes: global AEs (GAEs¹), toroidal AEs (TAEs²), beta-induced AEs (BAEs³), and others. In addition, tokamak plasmas may be characterized by forced oscillations called energetic particle continuum modes (EPMs). These are non-normal modes of the SAW continuum spectra, merging as discrete fluctuations at the frequency that maximizes the wave-EP power exchange, above the threshold condition of the continuum damping.⁴ The excitation of these modes creates a transport channel for the EPs which can lead to losses of EPs before their thermalization, causing a less effective heating and also possibly damaging the vessel of the machine. They are also believed to be responsible for abrupt large events (ALE) observed in the Japanese tokamak (JT-60U).⁵

The present paper will detail the studies carried out on the Landau damping and on the continuum damping.^6,7 In order to retain the description of all the nonlinear effects (such as wave-wave and wave-particle interaction, as well as finite-Larmor-radius and finite-orbit-width effects), the simulations have been mainly performed with the gyrokinetic code ORB5. It is an electromagnetic, global, nonlinear, particle in cell (PIC) code^8,9 (whose model, if properly set, contains the MHD equations as a subset). When possible, the simulation results will be compared with the analytical theory. The analytical prediction of the continuum damping has been obtained studying the ideal MHD vorticity equation. For the Landau damping, an analytical estimation has been obtained by adding a kinetic term to the MHD vorticity equation and then following a perturbative approach. These two analytical derivations will be detailed in the dedicated sections. The simulations always consider global modes, fulfilling the condition $k_{⊥} ρ_{i} ≪ 1$ ⁠, where $k_{⊥}$ is the perturbation wave vector perpendicular to the background magnetic field and ρ_i is the bulk ion Larmour radius. Therefore, non-ideal effects of electron inertia and finite ion Larmor radius¹⁰ are considered negligible and no kinetic modification of the Alfvén modes is observed; thus, radiative damping is negligible.¹¹

The paper is structured as follows. In Sec. II, a description of the model implemented in the code ORB5 is given. Section III will be dedicated to the documentation of the phase mixing for perturbations of the continuum spectrum. Section IV will be dedicated to the documentation of the Landau damping in the absence of continuum damping, acting on TAEs. The numerical results shown in both Secs. III and IV will be obtained in the absence of EPs. In Sec. III, the numerical simulations will be performed in the cylinder limit using simplified profiles. In the simulations in Sec. IV, a small but finite inverse aspect ratio will be considered, using the equilibrium profiles of the International Tokamak Physics Activity (ITPA¹²). In Sec. V, the studies on the linear and nonlinear growth rate and frequency spectra conducted considering experimental profiles from the Non-Linear Energetic-particle Dynamics (NLED)- ASDEX-Upgrade (AUG) case¹³ will be presented. In this section, EPs have also been retained and have been observed to drive, depending on the shape of their density profile, an EPM or a TAE. TAEs driven by energetic ions have been treated in Ref. 14 (where all the particle species have been treated kinetically) and in Ref. 15 (where the bulk and fast ions have been treated as gyro-kinetic and, for the electron, a fluid model has been used). In the present paper, all the species are considered as drift-kinetic. In Sec. V, results of a comparison with the linear gyrokinetic code LIGKA¹⁶ are shown. Finally, a summary and future outlook will be given in Sec. VI.

II. MODEL

Since the Alfvén waves have a frequency much smaller than the typical ion cyclotron frequency $(Ω_{c i})$ and their amplitude in the core is small compared to the background quantities, a good description of their propagation and interaction with the bulk plasma can be given through gyrokinetic theory. This allows retaining a kinetic description of the events under consideration, reducing the 6D problem to a 5D one, by averaging the fast gyromotion. In this way, the numerical costs are sensibly reduced.

ORB5 is a global, nonlinear, gyrokinetic, electromagnetic, PIC code, which can take into account collisions and sources.^8,17 The gyrokinetic model of ORB5¹⁸ contains the reduced MHD equations as a subset.¹⁹ Additional discussion about the link between reduced MHD and gyrokinetics can be found in Refs. 20 and 21. In this section, we give a brief description of the gyrokinetic model implemented in ORB5 and briefly show how the implemented equations are solved. We refer for more exhaustive explanations to Ref. 17, which also gives a more complete description of the recent updates in ORB5. The magnetic equilibria used by ORB5 in this work are ideal-MHD equilibria (solution of the Grad-Shafranov equation) from the CHEASE code²² for the NLED-AUG case, or ad hoc equilibria constituted by circular, concentric magnetic surfaces. It deals with a straight-field line set of coordinates. The magnetic surfaces are labeled by $r = \sqrt{ψ / ψ_{edge}}$ ⁠, which plays the role of radial coordinates. Here, ψ is the poloidal magnetic flux function. The angular dependence is given by the toroidal coordinate $φ$ and by the poloidal magnetic angle,

χ = \frac{1}{q (r)} \int_{0}^{θ} d θ^{'} \frac{B \cdot \nabla φ}{B \cdot \nabla θ^{'}},

(1)

where $θ^{'}$ is the geometrical poloidal angle and q(s) is the safety factor, defined as

q (r) = \frac{1}{2 π} \int_{0}^{2 π} d θ^{'} \frac{B \cdot \nabla φ}{B \cdot \nabla θ^{'}} .

(2)

All the quantities in the code are normalized through four reference parameters: the ion mass (m_i), the ion charge (⁠ $q_{i} = e Z_{i}$ ⁠, with e the electric charge and Z_i the atomic number), the value of the magnetic field strength on axis (⁠ $B_{0} = | B (r = 0) |$ ⁠), and the value of the electron temperature at a specified reference position r₀, $T_{e} (r_{0})$ ⁠. All other normalized quantities are obtained through these: the time units are provided in the inverse of the ion-cyclotron frequency, $Ω_{c i} = q_{i} B_{0} / (m_{i} c)$ ⁠, the velocity units are normalized through the ions sound velocity (⁠ $c_{s} = \sqrt{q_{s} T_{e} (r_{0}) / m_{i}}$ ⁠, with the temperature measured in keV), the length units through the ion sound Larmor radius (⁠ $ρ_{s} = c_{s} / Ω_{c i}$ ⁠), and the densities are normalized by means of their average in space. The Vlasov–Maxwell gyrokinetic equations are derived through variational principles from a gyrokinetic Lagrangian. The Lagrangian is discretized through finite elements and PIC methods. An Action principle is then applied to the discrete Lagrangian itself leading to the discrete gyrokinetic Vlasov and field equations. Those equations reflect the conservation properties of the discrete Lagrangian in the limit $Δ t \to 0$ ⁠. An immediate consequence is that it is possible to consistently derive discrete conserved quantities (such as the energy²³) that are also used in ORB5 to test the quality of the simulations performed. In the ordering is present, separating the effects given from the geometry of the non-uniform magnetic field, from those related to the fluctuations of the electromagnetic perturbation. This means that (as can be derived^18,24) the small parameter related to the variation of the background magnetic field $ε_{B} = ρ_{t h} L_{B}$ (with ρ_th the thermal Larmor radius and L_B the typical variation length of the magnetic field) and the small parameter related to the fluctuating electromagnetic field (⁠ $ε_{δ}$ ⁠) are related through

ε_{B} = O (ε_{δ}^{2}) .

(3)

In this way, the action functional, written in “p_z-formulation,” appears as the following:

\begin{matrix} A = \int_{t_{0}}^{t_{1}} L d t = \sum_{s} \int d t d Ω (\frac{q_{s}}{c} A^{*} \cdot \dot{X} + \frac{m_{s} c}{q_{s}} μ \dot{θ} - H_{0}) f_{s} \\ - ε_{δ} \sum_{s = e} \int d t d Ω H_{1} f_{s} - ε_{δ} \int d t d Ω H_{1}^{d k} f_{e} \\ - ε_{δ}^{2} \sum_{s = e} \int d t d Ω H_{2} f_{e q, s} - α ε_{δ}^{2} \int d t d Ω H_{2}^{d k} f_{e q, e} - α ε_{δ}^{2} \int d t d V \frac{| \nabla_{⊥} A_{1, | |} |^{2}}{8 π}, \end{matrix}

(4)

where α = 0 gives the electrostatic model and α = 1 the electromagnetic one. In Eq. (4), $d Ω = d V d W$ ⁠, with $d W = B_{| |}^{*} d μ d p_{z}$ ⁠. A sum over the species “s” also appears. The symplectic magnetic field is defined through the symplectic magnetic potential $A^{*} = A + (c / q_{s}) p_{z} \hat{b}$ ⁠, with $\hat{b}$ the unit vector parallel to the background magnetic field. The canonical gyrocenter momentum is $p_{z} = m_{s} v_{∥} + α ε_{δ} (q_{s} / c) A_{1 ∥}$ ⁠. In the action functional, some approximations have been made. The quasi-neutrality allows considering in Eq. (4) only the contribution given from the magnetic potential, neglecting the one given from the perturbed electric field. Also the incompressibility of the parallel perturbed magnetic field is assumed $B_{1, | |} = o (B_{1, ⊥})$ and only the perpendicular component of the perturbed magnetic potential is retained: $B_{1, | |} = \nabla \times (A_{1, | |} b) \sim \nabla A_{1, | |} \times b$ ⁠. In Eq. (4), it must be noted that while H₀, H₁ multiply the total distribution functions f_s, f_e, H₂ is related only to the equilibrium distribution function. Thanks to this choice, nonlinear second order terms do not appear in the gyrocenter dynamics and the field equations are linear. The gyrocenter Hamiltonians appearing are

\begin{matrix} H_{0} = \frac{p_{z}^{2}}{2 m_{s}} + μ B H_{1} = q_{s} ⟨ ϕ_{1} - α A_{1, ∥} \frac{p_{z}}{m_{s} c} ⟩, \\ H_{2} = - \frac{m_{s} c^{2}}{2 B^{2}} | \nabla_{⊥} ϕ_{1} |^{2} + α \frac{q_{s}^{2}}{2 m_{s} c^{2}} {⟨ A_{1, ∥} ⟩}^{2}, \end{matrix}

(5)

where the gyroaveraging operator has been introduced $⟨ f ⟩ = \frac{1}{2 π} \int_{0}^{2 π} d θ f$ ⁠. The gyroaveraging is removed for the electrons that are treated as drift-kinetic,

H_{1}^{d k} = - e (ϕ_{1} - α A_{1, ∥} \frac{p_{z}}{m_{s} c}) H_{2}^{d k} = α \frac{e^{2}}{2 m_{e} c^{2}} A_{1, ∥}^{2} .

(6)

For the distribution of the species s, the linear gyrokinetic Vlasov equation is

\frac{d f_{s}}{d t} = \frac{\partial f_{s}}{\partial t} + \dot{X} \cdot \nabla f_{s} + {\dot{p}}_{z} \frac{\partial f_{s}}{\partial p_{z}} = 0,

(7)

where the gyrokinetic characteristics can be derived from Eq. (4) and are

{\begin{matrix} \dot{X} = \frac{c \hat{b}}{q_{s} B_{∥}^{⋆}} \times \nabla H + \frac{\partial H}{\partial p_{z}} \frac{B^{⋆}}{B_{∥}^{⋆}}, \\ {\dot{p}}_{z} = - \frac{B^{⋆}}{B_{∥}^{⋆}} \cdot \nabla H . \end{matrix}

(8)

The field equations, quasineutrality and Ampère, are both derived from Eq. (4) via functional derivatives on the perturbed field. ORB5 splits the total distribution function in a background distribution function f₀ and a time dependent δf and discretize, this latter through numerical particles (markers) used to sample the phase space. Through an operator splitting approach, the code solves first the conlisionless dynamics (using a 4th-order Runge–Kutta method). The quasineutrality and Ampère equations are solved using the Galerkin methods and the perturbed fields are discretized through cubic B-splines finite elements defined on a grid $(N_{s}, N_{χ}, N_{ϕ})$ ⁠. Finally, it is important to mention that recently the mixed-representation (“pullback” scheme²⁵) has solved the so-called “cancelation problem” for electromagnetic simulations.

III. PHASE MIXING/CONTINUUM DAMPING

In the present section, low temperatures are chosen to minimize the Landau damping. Here, we study the process of phase mixing of an initial perturbation oscillating at the frequency of the continuum, which leads to a decay of the perturbation, scaling $\sim 1 / t$ ⁠. Following, e.g., Ref. 26, we will call this process (local) continuum damping. In the literature, an estimation of the global continuum damping can be found, referring to global eigenmodes with frequency within the gap, touching the shear Alfvén continuum and having an exponential decay rate.^27,28 This is not what is considered in this section.

In order to do so, it is first important to understand what are the equations governing the Alfvén waves. These will be obtained under the validity of the ideal magnetohydrodynamic (MHD) theory by treating MHD equations with a perturbative approach. The Alfvén wave's dynamics can be expressed starting from the quasi-neutrality condition $\nabla \cdot δ J = 0$ (⁠ $δ J$ being the perturbed current) that is rewritten in terms of its components parallel and perpendicular to the background magnetic field (⁠ $\hat{b} = B / B$ ⁠), yielding

\nabla \cdot δ J_{⊥} + B \cdot \nabla \frac{δ J_{∥}}{B} = 0.

(9)

An expression for the vorticity equation can be elegantly obtained following Ref. 29, where the Hain and Lust equation is derived. In the present paper, the authors will follow Ref. 30. In order to obtain a simplified but relevant set of equations, modes with $k_{⊥} ≫ k_{∥}$ are considered, so that the timescale between incompressible shear Alfvén waves and compressional waves can be separated. To further simplify the problem, we consider a pressureless plasma (P = 0) obtaining the following vorticity equation:

B \cdot \nabla [\frac{1}{B} \nabla_{⊥}^{2} (\frac{1}{B} B \nabla δ ϕ)] - \nabla \cdot (\frac{1}{v_{A}^{2}} \frac{\partial^{2}}{\partial t^{2}} \nabla_{⊥} δ ϕ) = 0.

(10)

A differential equation for the perturbed scalar potential $δ ϕ$ is thus obtained. It is linked to the perturbed magnetic potential (⁠ $δ A \approx δ A \hat{b}$ ⁠) through the condition $δ E_{∥} = 0$ ⁠, derived from the ideal Ohm's law. In this section, a non-uniform plasma equilibrium with a cylindrical limit will be considered. a will denote the typical length scale perpendicular to the equilibrium magnetic field while R₀ will represent the typical length scale parallel to it. The equilibrium magnetic field, in a coordinate system $(r, θ, z)$ ⁠, will be assumed to be $B = (0, B_{0, θ} (r), B_{0, z} (r))$ ⁠. By assuming a shear Alfvén oscillation of the scalar potential $δ ϕ (r, θ, ϕ, t)$ of the form

δ ϕ = \sum_{m, n} δ ϕ_{m, n} (r) e^{i (m θ - \frac{n z}{R_{0}} - ω t)},

(11)

where m is the poloidal mode number, we can now write Eq. (10) in cylindrical coordinates as follows:

\frac{1}{r} \frac{\partial}{\partial r} r [{(\frac{m}{q (r)} - n)}^{2} + \frac{R_{0}^{2}}{v_{A}^{2}} \frac{\partial^{2}}{\partial t^{2}}] \frac{\partial}{\partial r} (\frac{δ ϕ}{r}) = \frac{m^{2}}{r^{2}} [{(\frac{m}{q (r)} - n)}^{2} + \frac{R_{0}^{2}}{v_{A}^{2}} \frac{\partial^{2}}{\partial t^{2}}] δ ϕ,

(12)

where the local safety-factor has been defined,

q (r) = \frac{r B_{0, z}}{R_{0} B_{0, θ}} .

(13)

The shear Alfvén wave dispersion relation is then found to be

ω_{A}^{2} = v_{A}^{2} k_{m, n}^{2} = \frac{v_{A}^{2}}{R_{0}^{2}} {(\frac{m}{q (r)} - n)}^{2} .

(14)

Equation (14) shows that the shear Alfvén waves are local plasma oscillations, having a frequency spectrum that varies continuously throughout the plasma radial direction. Exploiting the local nature of the continuum plasma oscillation, Eq. (12) can be reduced to

\frac{1}{r} \frac{\partial}{\partial r} r [{(\frac{m}{q (r)} - n)}^{2} - ω^{2} \frac{R_{0}^{2}}{v_{A}^{2}}] \frac{\partial δ ϕ}{\partial r} = 0,

(15)

which integrated in the radial domain, becomes a differential equation for the radial electric field E_r,

(ω_{A}^{2} + \frac{\partial^{2}}{\partial t^{2}}) E_{r} = 0 \Rightarrow E_{r} = E_{0} e^{- i ω_{A} (r) t} .

(16)

Assuming now a dispersion relation of the form $ω_{A} (r) = ω_{A 0} + ω_{A}^{'} (r - r_{0})$ and by Fourier transforming the radial electric field in the radial coordinate, the following relation is obtained:

(f E_{r}) (k_{r}) = \sqrt{2 π} E_{0} e^{- i (ω_{A 0} - ω^{'} r_{0}) t} δ (k_{r} + ω_{A}^{'} t) k_{r} \propto - ω_{A}^{'} t .

(17)

The obtained linear dependence in time of the radial wave number is a proof of the phase mixing together with the fact that, as $E_{r} (r, t) = - i k_{r} (t) ϕ (r, t)$ ⁠, the scalar potential exhibits the characteristic decay called continuum damping,

δ ϕ \propto {| ω_{A}^{'} t |}^{- 1}

(18)

as was proven in Ref. 30 (see also Refs. 26, 31, and 32 for the application to Geodesic Acoustic Modes, GAMs).

The simulations presented in this section have considered simplified geometry and profiles. EPs are not involved in the simulations here presented and the bulk ions are hydrogen ions. The inverse aspect ratio is $ε = 0.01$ ⁠. Flat density (⁠ $n_{e} = n_{H} = 2.22 \times 10^{20} m^{- 3}$ ⁠) and temperature (⁠ $T_{e} = T_{H} = 0.01 keV$ ⁠) profiles have been taken into account. This choice leaves all the radial dependence of the dispersion relation in the safety factor profile, which exhibits a linear radial dependence: $q = q_{0} + q_{1} \cdot r$ ⁠, so that: $ω_{A} = \frac{v_{A}}{R_{0}} \frac{1}{q_{0} + q_{1} r}$ ⁠. Moreover, in this temperature regime the Landau damping can be neglected.³³ The involved equilibrium corresponds to an Alfvén continuum (see Fig. 1) for a mode perturbation with poloidal and toroidal mode numbers (m, n) = (1, 0). In Table I are listed: minor and major radii, value on axis of the equilibrium magnetic field, ion cyclotron frequency $(Ω_{c i})$ ⁠, Alfvén frequency on-axis (⁠ $ω_{A 0} = ω_{A} (r = 0)$ ⁠), and ratio of the last two and order of magnitude of $k_{⊥} ρ_{i}$ ⁠. In Table II are listed the main parameters of the simulation in use, that is: spatial grid, time step, and number of chosen ions and electrons.

FIG. 1.

View large Download slide

Frequency spectra. The results of ORB5 are given by the blue line. The theoretical prediction for ω_A, corresponding to Eq. (14), corresponds to the dotted line.

TABLE I.

Important constants used in the simulations.

$a_{0} (m)$	$R_{0} (m)$	$B_{0} (T)$	$Ω_{c i} (rad / s)$	$ω_{A 0} (rad / s)$	$Ω_{c i} / ω_{A 0}$	$k_{⊥} ρ_{i}$
0.1	10	3	$2.87 \times 10^{8}$	$4.38 \times 10^{5}$	655	$10^{- 3}$

$a_{0} (m)$	$R_{0} (m)$	$B_{0} (T)$	$Ω_{c i} (rad / s)$	$ω_{A 0} (rad / s)$	$Ω_{c i} / ω_{A 0}$	$k_{⊥} ρ_{i}$
0.1	10	3	$2.87 \times 10^{8}$	$4.38 \times 10^{5}$	655	$10^{- 3}$

TABLE II.

Main simulation parameters for the study of the continuum damping.

$(n_{r}, n_{χ}, n_{φ})$	$Δ t (Ω_{c i}^{- 1})$	Number of ions	Number of electrons
(256, 256, 4)	20	10⁷	10⁷

In Fig. 2, the measured values for the wave numbers k_r are shown for a simulation having $q_{0} = 1.75$ and $q_{1} = 0.5$ ⁠. They have been measured by fitting the mode structure with a sinusoidal function at times where a maximum has been reached at the radial position r = 0.6. By linearly fitting the measured wave numbers, it is possible to extract the coefficient $k_{r, 1}$ [representing the slope of the line $k_{r} (t) = k_{r, 0} + k_{r, 1} t$ ]. $k_{r, 1}$ is found to be in reasonable agreement with the radial derivative of the frequency (⁠ $ω'_{A}$ ⁠), according to the theoretical expectations [Eq. (17)]. In Fig. 3, the dynamics of the scalar potential at some radial positions are shown, together with the predicted decay, Eq. (18).

FIG. 2.

View large Download slide

Radial wave number dependence on time. The results of ORB5 are given by dots. The theoretical prediction for this simulation (⁠ $q_{0} = 1.75$ and $q_{1} = 0.5$ ⁠) is that $k_{r, 1} = 0.123 ω_{A 0}$ ⁠, while the measured value is $k_{r, 1} = 0.107 ω_{A 0}$ ⁠.

FIG. 3.

View large Download slide

Perturbation amplitude dependence on time. Analytical estimations are given by the dashed lines (curves decaying in time as $ϕ \sim {| ω_{A}^{'} t |}^{- 1}$ ⁠). The scalar potential measured at different radial positions is given by continuous lines.

Finally, in Fig. 4 the obtained values of the coefficients $k_{r, 1}$ have been plotted against different values of the slope of the safety factor profiles (q₁) in use in the different simulations and compared with Eq. (17). Given the reasonable agreement found between the results of the numerical simulations and the theory, we can say to have verified the relevance of the continuum damping as the main damping mechanism for this specific case and to have observed the presence of phase mixing.

FIG. 4.

View large Download slide

Dependence of $k_{r, 1}$ on the slope of the safety factor profile. Analytical estimations are given by the dashed lines and the results of ORB5 are given by dots. No EPs are present here.

IV. LANDAU DAMPING

In the present section, the Landau damping will be studied. The attention will be focused also on a particular Alfvén eigenmode, the toroidal Alfvén eignemode (TAE). Its characteristic frequency lies in the forbidden frequency window (gap) of the shear Alfvén continuum created by the coupling between two close poloidal harmonics $(m, m + 1)$ ⁠, because of the finite tokamak toroidicity.³⁴ TAEs (like discrete AEs in general), can exist essentially free of continumm damping.^2,35 A TAE is located at a radial position r₀ satisfying: $q (r_{0}) = \frac{2 m + 1}{2 n}$ ⁠. The theoretical derivation in Refs. 34 and 36 will be now followed. Here, a kinetic transverse part of the wave-induced current $δ J_{⊥}^{k}$ is added to the ideal MHD current, so that Eq. (9) becomes

\nabla \cdot (δ J^{MHD} + δ J_{⊥}^{k}) = 0.

(19)

Equation (19) is then multiplied by $δ ϕ$ and integrated in the overall plasma volume obtaining

\int d x δ J^{MHD} \cdot \nabla δ ϕ + \int d x δ J_{⊥}^{k} \cdot \nabla_{⊥} δ ϕ = 0,

(20)

where it was assumed as boundary conditions $\int d x \cdot δ J δ ϕ = 0$ ⁠. Calling ω₀ the frequency of the wave solution of the ideal MHD vorticity equation, we can consider $ω = ω_{0} + δ ω$ the solution of the new vorticity equation Eq. (19), being $δ ω ≪ ω_{0}$ ⁠. Following a perturbative approach, an expression for $γ = ℑ ω$ is obtained from Eq. (20),

γ = \frac{2 π}{c^{2}} \frac{\sum_{m, n} \int_{V} d^{3} x δ J_{⊥ m, n}^{k} \cdot \nabla δ ϕ_{m, n}^{*}}{\sum_{m, n} \int_{V} d^{3} x \frac{1}{v_{A}^{2}} [| δ ϕ'_{m, n} |^{2} + {(\frac{m}{r})}^{2} | δ ϕ_{m, n} |^{2}]} .

(21)

Subscripts m and n denote the poloidal and toroidal mode numbers, respectively. In order to obtain a simplified equation for γ, some further calculations have been done and will be now described. Assuming a Maxwellian distribution function F₀ and focusing our attention on TAE (that is assuming to have a perturbation $δ ϕ$ strongly peaked at the radial position where we expect to have a TAE), we obtain in cylindrical coordinates,

γ = \sum_{j} γ_{j} γ_{j} = - β_{j} q_{0}^{2} \frac{v_{A}}{2 q_{0} R_{0}} [G_{m j} + n q_{0} r_{L θ, j} \frac{1}{n_{0, j}} \frac{\partial n_{0, j}}{\partial r} (H_{m j} + η J_{m j})],

(22)

with

Ω_{θ, j} = \frac{e B_{p}}{m_{j} c}, r_{L θ} = \frac{v_{t h, j}}{Ω_{θ, j}}, β_{j} = 8 π \frac{n_{0, j} T_{j}}{B_{0}^{2}}, η_{j} = \frac{\partial log (T_{j})}{\partial log (n_{0, j})}, λ_{j} = v_{A} / v_{t h, j}

(23)

and

{\begin{matrix} g_{m, j} (λ_{j}) = \frac{π}{2} λ_{j} (1 + 2 λ_{j}^{2} + 2 λ_{j}^{4}) e^{- λ_{j}^{2}}, G_{m j} = g_{m, j} (λ_{j}) + g_{m, j} (λ_{j} / 3), \\ h_{m, j} (λ_{j}) = \frac{π}{2} (1 + 2 λ_{j}^{2} + 2 λ_{j}^{4}) e^{- λ_{j}^{2}}, H_{m j} = h_{m, j} (λ_{j}) + \frac{1}{3} h_{m, j} (λ_{j} / 3), \\ j_{m, j} (λ_{j}) = \frac{π}{2} (\frac{3}{2} + 2 λ_{j}^{2} + λ_{j}^{4} + 2 λ_{j}^{6}) e^{- λ_{j}^{2}}, J_{m j} = j_{m, j} (λ_{j}) + \frac{1}{3} j_{m, j} (λ_{j} / 3) . \end{matrix}

(24)

In Eq. (22), γ has been decomposed in the species contributions (the sum over j). It is formally identical to the one derived in Ref. 37. The difference lies in the fact that in Ref. 37 the authors have obtained the estimation for γ starting from energy principles, while here everything has been done by adding a correction to the MHD quasi-neutrality equation and thus to the Alfvén dynamics.³⁸ Since we are interested in the study of the Landau damping, we will not consider the EP contribution, which actually drives the mode unstable. Equation (22) depends on the ratio between the Alfvén speed and the thermal velocity of the considered species. In this section, an equilibrium with small, but finite value of inverse aspect ratio will be considered, $ε = 0.1$ ⁠. The temperature and density profiles are flat. The magnetic equilibrium and profiles are those of the ITPA-TAE international benchmark case¹² and the safety factor profile is shown in Fig. 5. In Fig. 6, the dependence of the damping against the electron mass is shown, together with the analytical theory. In the present section, the bulk ions are hydrogens ions, so the realistic mass ratio among the bulk species is $m_{H} / m_{e} \sim 2000$ ⁠. As can be noticed, the electron Landau damping depends sensitively on its mass, but it shows a weaker dependence than expected compared to the analytical theory. One should note that the analytical theory presented here considers only well passing particles, whereas the contribution of the species changes depending on weather the resonant particles are passing, trapped or in between. It was proved in Ref. 39 that barely trapped electrons are the relevant ones for EP driven modes. We suppose this to be valid also for the Alfvén modes considered here and to be the main justification for the observed discrepancy between the measured damping and the analytical theory. It would be of great interest to develop a more accurate theory capable of better fittings the results, but this will represent a task for a future work. For computational reasons, the electron mass has been then taken 200 times lighter than the ions mass (⁠ $m_{e} = m_{H} / 200$ ⁠). Nevertheless, we can notice that the chosen mass affects the final result of approximately 30% of its value (see Fig. 6) and, for the purpose of this paper, this is the required precision we have chosen.

FIG. 5.

View large Download slide

Safety factor of the ITPA-TAE case (⁠ $q (r) ≃ 1.71 + 0.15 r^{2}$ ⁠).

FIG. 6.

View large Download slide

Damping rate dependence on the electron mass for ORB5 simulations.

Scans against temperature have been performed. The relative importance of electron and ion Landau damping depends on electron and ion temperatures. In this regime, the species temperatures are kept flat and, at the reference case, they are equal (T_H = T_e). In Fig. 7, the dependence of the Landau damping against the electron temperature is shown and the results are compared with the analytical theory. Scans against the ion temperature have also been performed (keeping the reference value of the electron temperature fixed) and no substantial variation has been observed. The observed discrepancy in Fig. 7 between the analytical theory and the simulations results is motivated using the arguments discussed above.

FIG. 7.

View large Download slide

Landau damping dependence vs the electron temperature.

In Tables III and IV, other important details of the simulations are presented.

TABLE III.

Important constants in use in the ITPA-TAE case.

$a_{0} (m)$	$R_{0} (m)$	$B_{0} (T)$	$Ω_{c i} (rad / s)$	$ω_{A 0} (rad / s)$	$Ω_{c i} / ω_{A 0}$	$k_{⊥} ρ_{i}$
1	10	3	$2.87 \times 10^{8}$	$1.46 \times 10^{6}$	196	$10^{- 2}$

$a_{0} (m)$	$R_{0} (m)$	$B_{0} (T)$	$Ω_{c i} (rad / s)$	$ω_{A 0} (rad / s)$	$Ω_{c i} / ω_{A 0}$	$k_{⊥} ρ_{i}$
1	10	3	$2.87 \times 10^{8}$	$1.46 \times 10^{6}$	196	$10^{- 2}$

TABLE IV.

Main simulation parameters of the ITPA-TAE case. Number of grid points used, time step, and number of markers involved.

$(n_{r}, n_{χ}, n_{φ})$	$Δ t (Ω_{c i}^{- 1})$	Number of ions	Number of electrons
(256, 288, 48)	5	$4 \times 10^{6}$	$32 \times 10^{6}$

The chosen initial potential perturbation is peaked around r = 0.5 and is constituted by one single toroidal mode number n = 6 while the poloidal mode numbers $9 \leq m \leq 12$ are considered (see Fig. 8). A TAE is located at r = 0.5 in the gap of the continuum spectra created by the coupling of the poloidal modes m = 10 and m = 11, as it is shown in Fig. 9. In Fig. 7, the dependence of the damping rate against the value of the electron temperature is shown. The damping rate is found to increase with the increasing electron temperature. This is evidence that the dominant damping is the electron Landau damping. The errorbar of the measured points corresponds to 20% of their value. This because, as it is shown in Fig. 10, the damping rate value has a dependence on the chosen width of the perturbation. For completeness in Fig. 7, the approximated analytical electron Landau damping formula is also shown (dashed line). A reasonable qualitative agreement is found between the predicted decay and the simulation results. Finally in Fig. 6, the dependence of the measured damping rate of ORB5 simulations against the electron mass has been shown. For decreasing electron masses, the absolute value of the damping rate is shown to decrease, consistently with theory of the electron Landau damping. In summary, it has been shown that the bulk electrons provide the main damping mechanism of the observed Alfvén modes in this particular regime. Several approximations have been done in the analytical theory, “inter alia” only passing particles are considered thus neglecting the contributions of barely trapped electrons, which are thought to be important, and which are included in our numerical simulations.

FIG. 8.

View large Download slide

Mode structure without EPs.

FIG. 9.

View large Download slide

Frequency spectra without EPs.

FIG. 10.

View large Download slide

Damping rate dependence on the width of the initial Gaussian beam.

V. NLED-AUG CASE

In the present section, the results of numerical simulations involving a realistic scenario will be presented.

The shot number $# 31213$ of ASDEX-Upgrade (AUG) has been selected as a benchmark within the Non-Linear Energetic-particle Dynamics (NLED) Eurofusion enabling research project.¹³ Here, an early off-axis neutral beam injection (NBI), with injection energy of ∼93 keV occurs with an injection angle (angle between the horizontal axis and the beam-line) of $7.13 °$ ⁠. The magnetic equilibrium measured at the time $t = 0.84 s$ is considered in the present simulations (see Fig. 11). This case is referred to as “NLED-AUG case.” Further description of this case can be found in Ref. 13. The NLED-AUG case is found to be of great interest because of its rich linear and nonlinear dynamics arising from the interaction of the modes with the EPs.

FIG. 11.

View large Download slide

Poloidal view of the magnetic equilibrium in use.

Tables V and VI contain the details of the main parameters considered in the simulations. Table VII shows the values of the bulk species profiles on axis in the absence of EPs. The bulk ions, as well as the EPs (when considered), are constituted by deuterium ions. The EP temperature will be considered to be radially flat and equal to $T_{E P} = 113 keV$ ⁠. For the EP density profiles, an off-axis density profile (see Fig. 12) fitting the experimental profiles is considered with a Maxwellian distribution function. For comparison, we also run simulations with an on-axis EPs density profiles (see Fig. 13). Note that when EPs are included the electron density profiles are changed in order to match quasi-neutrality $n_{e} = Z_{D} \cdot n_{D} + Z_{E P} \cdot n_{E P}$ ⁠.

TABLE V.

Important constants in use in the NLED-AUG case.

$a_{0} (m)$	$R_{0} (m)$	$B_{0} (T)$	$Ω_{c i} (rad / s)$	$ω_{A 0} (rad / s)$	$Ω_{c i} / ω_{A 0}$	$k_{⊥} ρ_{i}$
0.482	1.666	2.202	$1.0539 \times 10^{8}$	$4.98 \times 10^{6}$	21.15	$10^{- 2}$

$a_{0} (m)$	$R_{0} (m)$	$B_{0} (T)$	$Ω_{c i} (rad / s)$	$ω_{A 0} (rad / s)$	$Ω_{c i} / ω_{A 0}$	$k_{⊥} ρ_{i}$
0.482	1.666	2.202	$1.0539 \times 10^{8}$	$4.98 \times 10^{6}$	21.15	$10^{- 2}$

TABLE VI.

Main simulation parameters of the NLED-AUG case. Number of grid points used, time step, and number of markers involved.

$(n_{r}, n_{χ}, n_{φ})$	$Δ t (Ω_{c i}^{- 1})$	Number of bulk ions	Number of electrons	Number of EPs
(256, 288, 48)	1	$16 \times 10^{6}$	$64 \times 10^{6}$	$16 \times 10^{6}$

$(n_{r}, n_{χ}, n_{φ})$	$Δ t (Ω_{c i}^{- 1})$	Number of bulk ions	Number of electrons	Number of EPs
(256, 288, 48)	1	$16 \times 10^{6}$	$64 \times 10^{6}$	$16 \times 10^{6}$

TABLE VII.

Profile parameters of the NLED-AUG case.

$T_{elelctrons} (s = 0) (keV)$	$T_{bulk ions} (s = 0) (keV)$	$n_{elelctrons} (s = 0) (m^{- 3})$	$n_{bulk ions} (s = 0) (m^{- 3})$	$n_{E P} (s = 0) (m^{- 3})$
0.709	2.48	$1.672 \times 10^{19}$	$1.6018 \times 10^{19}$	$6.98 \times 10^{17}$

$T_{elelctrons} (s = 0) (keV)$	$T_{bulk ions} (s = 0) (keV)$	$n_{elelctrons} (s = 0) (m^{- 3})$	$n_{bulk ions} (s = 0) (m^{- 3})$	$n_{E P} (s = 0) (m^{- 3})$
0.709	2.48	$1.672 \times 10^{19}$	$1.6018 \times 10^{19}$	$6.98 \times 10^{17}$

FIG. 12.

View large Download slide

Off-axis EPs density profile. EP concentration of 3%.

FIG. 13.

View large Download slide

On-axis EP density profile. EP concentration of 3%.

In Fig. 14, the safety factor profile is shown. In Fig. 15, the temperature profiles of the bulk species are shown. The safety factor profile has a reversed shear, with $q_{min} (r = 0.5) ≃ 2.28$ ⁠. In Fig. 16, the dependence of the growth rate against the electron mass is shown. Since in this case ions are constituted by deuterium ions, the realistic mass ratio among the bulk species is $m_{D} / m_{e} \sim 3600$ ⁠. For numerical reasons, the electron mass is chosen to be $m_{e} = m_{D} / 500$ ⁠, with m_D the deuterium mass. As can be seen from Fig. 16, the chosen fictitious mass increases the real value of the growth rate by approximately 20%. This represents the chosen precision for the studies that we are carrying.

FIG. 14.

View large Download slide

Safety factor profile.

FIG. 15.

View large Download slide

Bulk species temperature profiles.

FIG. 16.

View large Download slide

Growth rate dependence against the electron mass.

The modes observed in the present scenario, strongly depend on the shape of the EP density profile. In particular, as will be shown, a TAE is observed in the presence of an on-axis density profile. An EPM is as well observed when an off-axis density profile is taken into account. The present section is divided into three subsections. The first and the second will retain only the linear dynamics. In the first, the dependence of γ against the electron temperature will be studied, with and without EP contribution. When the EPs are here considered, they will assume an on-axis density profile and TAEs will be driven unstable. The results of a benchmark with LIGKA will also be presented. In Subsection V B, the presence of an EPM, driven unstable by EPs assuming an off-axis density profile, will be shown. Finally in Subsection V C, the results of simulations involving the nonlinear dynamics will be presented.

It is important to remind that, unless specified, the bulk and energetic ions will be treated as drift-kinetic. The initial perturbation considered will take into account just one toroidal mode number (n = 1) and the poloidal mode number $0 \leq m \leq 7$ ⁠.

A. Linear TAEs

Figures 17, 18, and 19 show the frequency spectra, mode structure, and poloidal view of the scalar potential $ϕ$ ⁠, obtained considering an on-axis density profile for the EPs (see Fig. 13) and a concentration equal to $⟨ n_{E P} ⟩ / ⟨ n_{e} ⟩ = 3 %$ (where $⟨ \dots ⟩$ indicates the volume average). In Fig. 17, the continuum spectra obtained with the linear gyrokinetic code LIGKA¹⁶ is shown (red crosses), together with the analytical curve for the continuum spectra calculated in cylindrical coordinates and including the toroidicity effects³⁴ (green dotted line).

FIG. 17.

View large Download slide

Frequency spectra measured with an on-axis EP density profile.

FIG. 18.

View large Download slide

Mode structure observed with an on-axis EP density profile.

FIG. 19.

View large Download slide

Observed scalar potential in the poloidal plane with an on-axis EP density profile.

The radial dependence of the EP density profiles is expressed by the formula: $n_{E P} ≃ {(1 - r^{α})}^{β}$ ⁠. The coefficients α, β have been chosen in order to have the maximum gradient of n_EP at the position where an Alfvén mode is expected. The numerical analysis shows that a mode lying in the gap of the continuum spectra created by the poloidal modes m = 2 and m = 3 is observed. It appears to be peaked at the radial position $r ≃ 0.738$ ⁠. Due to the radial localization and frequency, this is identified as a TAE.

In Fig. 20, the dependence of the growth rate against the electron temperature, keeping the bulk ion temperature constant (⁠ $T_{D} (r = 0) = 3.5 keV$ ⁠), is shown. In Fig. 21, the dependence of the growth rate against the bulk ion temperature, keeping the electron temperature constant (⁠ $T_{e} (r = 0) = 0.707 keV$ ⁠), is shown. The EP temperature is flat $T_{E P} = 113 keV$ and the EPs have a concentration of 3%. In Fig. 22, the dependence of the damping rate (simulations without fast particles) is shown, against the value of the electron temperature. This study of the dependence of the growth or damping rate against the electron temperature shows that the electrons are mainly responsible for the damping of the Alfvén modes even in this realistic scenario, which is identified here as electron Landau damping. Since a realistic scenario is considered here, the approximate theoretical predictions for the Landau damping described in Sec. IV is outside its validity regime and therefore is not shown.

FIG. 20.

View large Download slide

Scan in the electron temperature for a growing mode. EP concentration 3%.

FIG. 21.

View large Download slide

Scan in the bulk ions temperature for a growing mode. EP concentration 3%.

FIG. 22.

View large Download slide

Scan in the electron temperature for a damped mode.

In Figs. 23 and 24, the results of a first benchmark between ORB5 and the code LIGKA are shown. Here, a scan in the EPs temperature is depicted. The EPs have an on-axis profiles. The dependence of growth rate and frequency against T_EP shown by the two codes is the same. The damping rate of the modes (measured at $T_{E P} = 30 keV$ ⁠), is $γ_{ORB 5} = - 0.003 ω_{A 0}$ and $γ_{LIGKA} = - 0.00215 ω_{A 0}$ ⁠. The observed differences between the results of LIGKA and ORB5 are mainly attributed to differences present in the equilibrium reconstruction underlying both codes. Moreover, because the mode is truly global (touching both the inner axis and the edge), the details at the boundary conditions are relevant. The observed discrepancies are another cause of the differences of the measurements between the two codes.

FIG. 23.

View large Download slide

Scan in T_EP. TAE growth rate, calculated with LIGKA and ORB5 for an EP concentration equal to 3% (same density, temperature profiles in use).

FIG. 24.

View large Download slide

Scan in T_EP. TAE frequency, calculated with LIGKA and ORB5 for an EP concentration equal to 3% (same density, temperature profiles in use).

B. Linear EPM

Figures 25, 26, and 27 show the frequency spectra, mode structure, and poloidal view of the scalar potential $ϕ$ ⁠, obtained considering an off-axis density profile for the EPs (see Fig. 12). A mode sitting at the radial position $r ≃ 0.22$ is observed. The dominant poloidal component of the scalar potential appears to be that having m = 2. Since the frequency observed in this simulation lies on the lower branch of the continuum, this mode is identified as an EPM. With $T_{E P} = 113 keV$ ⁠, the measured frequency is observed to change from 135 to $150 kHz$ when finite Larmor radius effects are taken into account. Considering a closer value to the experimental one for the EP temperature (⁠ $T_{E P} \sim 30 keV$ ⁠), obtained by fitting the pressure of the injected beam, the measured frequency turns out to be of $110 kHz$ ⁠. In Fig. 28, the spectrogram obtained with Mirnov coils is shown. A big variety of EP driven modes can be found. At $t = 0.84 s$ ⁠, the modes with frequencies around $50 kHz$ have been identified as EGAMs (see Refs. 39–41). We focus here on the Alfvén modes with frequency lying in the domain between 100 and $150 kHz$ ⁠. In Fig. 28, the values of the measured frequencies obtained in simulations keeping finite orbit effects are shown. The withe cross has been obtained in a simulation where $T_{E P} \sim 113 keV$ ⁠. The light-blue cross instead has been obtained in a simulation where $T_{E P} \sim 30 keV$ ⁠. Strong approximations have been done here for the EPs (flat temperature profile, Maxwellian distribution function). Despite that, the estimate of the frequency can be observed to lie in the range of the frequencies measured in the experiments. In order to have a more precise comparison, a more accurate EP distribution function will be implemented in ORB5 and results shown in a dedicated paper.

FIG. 25.

View large Download slide

Frequency spectra measured with an off-axis EPs density profile.

FIG. 26.

View large Download slide

Mode structure observed with an off-axis EP density profile.

FIG. 27.

View large Download slide

Observed scalar potential in the poloidal plane with an off-axis EP density profile.

FIG. 28.

View large Download slide

Experimental spectrogram obtained with Mirnov Coil compared with theoretical predictions at one selected time. The theoretical predictions are obtained treating the fast ions as gyrokinetic and having temperature of $\sim 30 keV$ (light-blue cross) and $\sim 113 keV$ (white cross).

C. Nonlinear simulations

In this subsection, results involving the nonlinear dynamics of the Alfvén waves are presented, when both on-axis (Figs. 29–32) and off-axis (Figs. 33–36) density profiles for the EPs are considered. With an on-axis density profile of the EPs, a mode sitting in the frequency gap is observed (TAE), Figs. 29 and 30. Its mode structure and frequency spectra are not observed to change passing from the linear to the nonlinear phase, confirming its nature of an eigenmode of this system, which is only weakly perturbed by EPs.

FIG. 29.

View large Download slide

Frequency spectra in the linear phase. EPs have an on-axis profile.

FIG. 30.

View large Download slide

Frequency spectra in the saturation phase. EPs have an on-axis profile.

FIG. 31.

View large Download slide

Mode structure in the linear phase. EPs have an on-axis profile.

FIG. 32.

View large Download slide

Mode structure in the saturation phase. EPs have an on-axis profile.

FIG. 33.

View large Download slide

Frequency spectra in the linear phase. EPs have an off-axis profile.

FIG. 34.

View large Download slide

Frequency spectra in the saturation phase. EPs have an off-axis profile.

FIG. 35.

View large Download slide

Mode structure in the linear phase. EPs have off-axis profile.

FIG. 36.

View large Download slide

Mode structure in the saturation phase. EPs have an off-axis profile.

When an off-axis density profile for the EPs is considered, a mode with dominant poloidal mode number m = 2 and peaked around $r ≃ 0.22$ is observed (see Fig. 35). This is consistent with what was observed in Sec. V B, when just the linear effects in the simulations were involved. Passing to the nonlinear phase a secondary mode with m = 2 and m = 3 is observed to grow around the radial position $r ≃ 0.738$ ⁠. This second mode is identified as the previously described TAE. This happens, because in the first linear phase the EPs drive the EPM unstable, which appears in fact to be dominant. In the nonlinear phase, the coexistence of the EPM and TAE is observed, due to an earlier saturation of the EPM (see Fig. 37).

FIG. 37.

View large Download slide

Time evolution of the dominant poloidal modes of the scalar potential (m = 2, 3) at the radial positions where the TAE and EPM are located. EPs have an off-axis profile.

VI. CONCLUSION

The presence of Alfvén modes in burning plasma can negatively affect the energy confinement and can also cause damage in the confining machine. Because of their importance, the present paper has dealt with the main damping mechanisms affecting the Alfvén modes, trying to outline them and to understand in which regime they are acting. The attention has been focused on toroidal Alfvén eigenmodes and energetic particle modes. These studies have been mainly carried by means of numerical simulations conducted with the code ORB5. The obtained results have been compared, when possible, with the presented analytical theory developed in a simplified geometry and with the results of the linear code LIGKA.

In Sec. III, simulations with a very small inverse aspect ratio (⁠ $ε = 0.01$ ⁠) have been considered in order to lead the analysis in the cylinder limit. Simplified profiles have been taken into account and very low electron temperature has been considered in order to have the continuum damping dominant over the Landau damping. The developed theory has been used to analyze the results of simulations without energetic particles. The dependence of the radial wave number of the mode against the time has been observed (phase mixing). Also, the scalar potential has been found to decay as $δ ϕ \propto {| ω_{A}^{'} t |}^{- 1}$ (continuum damping).

In Sec. IV, higher bulk ion and electron temperatures have been considered in order to observe the Landau damping to be dominant over the continuum damping. The numerical simulations have been conducted using plasma equilibrium and profiles from the ITPA-TAE international benchmark case. In order to separate the ion and electron contribution to the damping, the dependence of the damping rate against the bulk species temperature has been studied. The results obtained analyzing the slope of the scalar potential from numerical simulations have been compared with the analytical theory in use. This has been achieved adding to the ideal MHD vorticity equation a kinetic term and then, through a perturbative approach, a simplified estimation for the Landau damping has been obtained using cylindrical coordinates. A reasonable agreement has been found between the analytical theory in use and the numerical results. This has also proved that the electrons are mainly responsible for the damping.

In Sec. V, a realistic plasma equilibrium taken from a shot in ASDEX Upgrade has been considered. The results of the linear numerical simulations have shown the dependence of the damping rate against the bulk electron temperature describing, also in this case, the action of the Landau damping. A benchmark with the code LIGKA has shown good agreement for the frequency and growth rate dependence on the EP temperature. Finally, the nonlinear simulations have shown the interaction of an EPM and a TAE in the scenario with an off-axis EP density profile.

The future works will extend the developed theory in order to find a better agreement between the predicted estimation of the decaying mode and the numerical simulations. In Ref. 4, it was suggested that all the Alfvén fluctuations can be explained within the framework of a single general fishbone-like dispersion relation (GFLDR). This could represent the starting point to improving the analytical prediction of the damping rate which would be a very interesting analytical and numerical task. In this respect, comparison with the GFLDR would be enhanced by evaluation of the generalized inertia for general geometry, following the recent analysis of Ref. 42. Future and deeper benchmark with the code LIGKA and the Hybrid MHD-Gyrokinetic code HYMAGYC⁴³ would be of great interest in order to better understand the linear and nonlinear dynamics contained in the NLED-AUG case.

ACKNOWLEDGMENTS

Simulations presented in this work were performed on the CINECA Marconi supercomputer within the ORBFAST and OrbZONE projects.

One of the authors, F. Vannini, would like to thank Xin Wang, Zhixin Lu, and Gregorio Vlad for useful, interesting discussions and for great help provided in understanding Alfvén dynamics, Gyrokinetics, and MHD theory. The authors wish to acknowledge stimulating discussions with F. Zonca, G. Fogaccia, A. Könies, J. Gonzalez-Martin, and A. Di Siena. This work was partly performed in the frame of the “Multi-scale Energetic particle Transport in fusion devices” ER project.

This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training program 2014–2018 and 2019–2020 under Grant Agreement No. 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission.

References

1.

K.

Appert

,

R.

Gruber

,

F.

Troyuon

, and

J.

Vaclavik

,

Plasma Phys.

24

,

1147

(

1982

).

https://doi.org/10.1088/0032-1028/24/9/010

Google Scholar

Crossref

2.

C. Z.

Cheng

and

M. S.

Chance

,

Phys. Fluids

29

,

3695

(

1986

).

https://doi.org/10.1063/1.865801

Google Scholar

Crossref

3.

H. S.

Zhang

,

Z.

Lin

,

I.

Holod

,

X.

Wang

,

Y.

Xiao

, and

W. L.

Zhang

,

Phys. Plasmas

17

,

112505

(

2010

).

https://doi.org/10.1063/1.3498761

Google Scholar

Crossref

4.

L.

Chen

and

F.

Zonca

,

Rev. Mod. Phys.

88

,

015008

(

2016

).

https://doi.org/10.1103/RevModPhys.88.015008

Google Scholar

Crossref

5.

A.

Bierwage

,

K.

Shinohara

,

Y.

Todo

,

N.

Aiba

,

M.

Ishikawa

,

G.

Matsunaga

,

M.

Takechi

, and

M.

Yagi

,

Nat. Commun.

9

,

3282

(

2018

).

https://doi.org/10.1038/s41467-018-05779-0

Google Scholar

Crossref

PubMed

6.

E.

Barston

,

Ann. Phys.

29

,

282

(

1964

).

https://doi.org/10.1016/0003-4916(64)90222-2

Google Scholar

Crossref

7.

Z.

Sedláček

,

J. Plasma Phys.

5

,

239

–

263

(

1971

).

https://doi.org/10.1017/S0022377800005754

Google Scholar

Crossref

8.

S.

Jolliet

,

A.

Bottino

,

P.

Angelino

,

R.

Hatzky

,

T.

Tran

,

B.

Mcmillan

,

O.

Sauter

,

K.

Appert

,

Y.

Idomura

, and

L.

Villard

,

Comput. Phys. Commun.

177

,

409

(

2007

).

https://doi.org/10.1016/j.cpc.2007.04.006

Google Scholar

Crossref

9.

A.

Bottino

,

T.

Vernay

,

B.

Scott

,

S.

Brunner

,

R.

Hatzky

,

S.

Jolliet

,

B. F.

McMillan

,

T. M.

Tran

, and

L.

Villard

,

Plasma Phys. Controlled Fusion

53

,

124027

(

2011

).

https://doi.org/10.1088/0741-3335/53/12/124027

Google Scholar

Crossref

10.

M. N.

Rosenbluth

and

P. H.

Rutherford

,

Phys. Rev. Lett.

34

,

1428

(

1975

).

https://doi.org/10.1103/PhysRevLett.34.1428

Google Scholar

Crossref

11.

R. R.

Mett

and

S. M.

Mahajan

,

Phys. Fluids B

4

,

2885

(

1992

).

https://doi.org/10.1063/1.860459

Google Scholar

Crossref

12.

A.

Könies

,

S.

Briguglio

,

N.

Gorelenkov

,

T.

Fehér

,

M.

Isaev

,

P.

Lauber

,

A.

Mishchenko

,

D.

Spong

,

Y.

Todo

,

W.

Cooper

,

R.

Hatzky

,

R.

Kleiber

,

M.

Borchardt

,

G.

Vlad

,

A.

Biancalani

, and

A.

Bottino

,

Nucl. Fusion

58

,

126027

(

2018

).

https://doi.org/10.1088/1741-4326/aae4e6

Google Scholar

Crossref

13.

P.

Lauber

, www2.ipp.mpg.de/pwl/nled_aug/data.html for further details on the NLED-AUG case.

14.

A.

Mishchenko

,

R.

Hatzky

, and

A.

Könies

,

Phys. Plasmas

15

,

112106

(

2008

).

https://doi.org/10.1063/1.3021453

Google Scholar

Crossref

15.

Y.

Nishimura

,

Phys. Plasmas

16

,

030702

(

2009

).

https://doi.org/10.1063/1.3088028

Google Scholar

Crossref

16.

P.

Lauber

,

S.

Günter

,

A.

Könies

, and

S.

Pinches

,

J. Comput. Phys.

226

,

447

(

2007

).

https://doi.org/10.1016/j.jcp.2007.04.019

Google Scholar

Crossref

17.

E.

Lanti

,

N.

Ohana

,

N.

Tronko

,

T.

Hayward-Schneider

,

A.

Bottino

,

B.

McMillan

,

A.

Mishchenko

,

A.

Scheinberg

,

A.

Biancalani

,

P.

Angelino

,

S.

Brunner

,

J.

Dominski

,

P.

Donnel

,

C.

Gheller

,

R.

Hatzky

,

A.

Jocksch

,

S.

Jolliet

,

Z. X.

Lu

,

J. P.

Martin Collar

,

I.

Novikau

,

E.

Sonnendrücker

,

T.

Vernay

, and

L.

Villard

,

Comput. Phys. Commun.

(published online,

2019

).

https://doi.org/10.1016/j.cpc.2019.107072

Google Scholar

18.

N.

Tronko

,

A.

Bottino

, and

E.

Sonnendruecker

,

Phys. Plasmas

23

,

082505

(

2016

).

https://doi.org/10.1063/1.4960039

Google Scholar

Crossref

19.

N.

Miyato

,

B. D.

Scott

,

D.

Strintzi

, and

S.

Tokuda

,

J. Phys. Soc. Jpn.

78

,

104501

(

2009

).

https://doi.org/10.1143/JPSJ.78.104501

Google Scholar

Crossref

20.

B.

Scott

and

J.

Smirnov

,

Phys. Plasmas

17

,

112302

(

2010

).

https://doi.org/10.1063/1.3507920

Google Scholar

Crossref

21.

J. A.

Krommes

,

Phys. Plasmas

20

,

124501

(

2013

).

https://doi.org/10.1063/1.4851996

Google Scholar

Crossref

22.

H.

Lütjens

,

A.

Bondeson

, and

O.

Sauter

,

Comput. Phys. Commun.

97

,

219

(

1996

).

https://doi.org/10.1016/0010-4655(96)00046-X

Google Scholar

Crossref

23.

A.

Bottino

and

E.

Sonnendrücker

,

J. Plasma Phys.

81

,

435810501

(

2015

).

https://doi.org/10.1017/S0022377815000574

Google Scholar

Crossref

24.

N.

Tronko

,

A.

Bottino

,

C.

Chandre

, and

E.

Sonnendruecker

,

Plasma Phys. Controlled Fusion

59

,

064008

(

2017

).

https://doi.org/10.1088/1361-6587/aa68af

Google Scholar

Crossref

25.

A.

Mishchenko

,

A.

Bottino

,

A.

Biancalani

,

R.

Hatzky

,

T.

Hayward-Schneider

,

N.

Ohana

,

E.

Lanti

,

S.

Brunner

,

L.

Villard

,

M.

Borchardt

et al.,

Comput. Phys. Commun.

238

,

194

–

202

(

2019

).

https://doi.org/10.1016/j.cpc.2018.12.002

Google Scholar

Crossref

26.

A.

Biancalani

,

F.

Palermo

,

C.

Angioni

,

A.

Bottino

, and

F.

Zonca

,

Phys. Plasmas

23

,

112115

(

2016

).

https://doi.org/10.1063/1.4967703

Google Scholar

Crossref

27.

F.

Zonca

and

L.

Chen

,

Phys. Rev. Lett.

68

,

592

(

1992

).

https://doi.org/10.1103/PhysRevLett.68.592

Google Scholar

Crossref

PubMed

28.

H. L.

Berk

,

J. W.

Van Dam

,

Z.

Guo

, and

D. M.

Lindberg

,

Phys. Fluids B

4

,

1806

(

1992

).

https://doi.org/10.1063/1.860455

Google Scholar

Crossref

29.

K.

Appert

,

R.

Gruber

, and

J.

Vaclavik

,

Phys. Fluids

17

,

1471

(

1974

).

https://doi.org/10.1063/1.1694918

Google Scholar

Crossref

30.

G.

Vlad

,

F.

Zonca

, and

S.

Briguglio

,

La Rivista del Nuovo Cimento

22

,

1

(

1999

).

https://doi.org/10.1007/BF02874568

Google Scholar

Crossref

31.

F.

Zonca

and

L.

Chen

,

Europhys. Lett.

83

,

35001

(

2008

).

https://doi.org/10.1209/0295-5075/83/35001

Google Scholar

Crossref

32.

F.

Palermo

,

A.

Biancalani

,

C.

Angioni

,

F.

Zonca

, and

A.

Bottino

,

Europhys. Lett.

115

,

15001

(

2016

).

https://doi.org/10.1209/0295-5075/115/15001

Google Scholar

Crossref

33.

L.

Chen

and

A.

Hasegawa

,

Phys. Fluids

17

,

1399

(

1974

).

https://doi.org/10.1063/1.1694904

Google Scholar

Crossref

34.

G. Y.

Fu

and

J. W.

Van Dam

,

Phys. Fluids B

1

,

1949

(

1989

).

https://doi.org/10.1063/1.859057

Google Scholar

Crossref

35.

C. Z.

Cheng

,

L.

Chen

, and

M. S.

Chance

,

Ann. Phys.

161

,

21

(

1985

).

https://doi.org/10.1016/0003-4916(85)90335-5

Google Scholar

Crossref

36.

Y. I.

Kolesnichenko

,

V. V.

Lutsenko

,

H.

Wobig

, and

V.

Yakovenko

,

Phys. Plasmas

9

,

517

(

2002

).

https://doi.org/10.1063/1.1432993

Google Scholar

Crossref

37.

R.

Betti

and

J. P.

Freidberg

,

Phys. Fluids B

4

,

1465

(

1992

).

https://doi.org/10.1063/1.860057

Google Scholar

Crossref

38.

Y.

Kolesnichenko

,

A.

Könies

,

V.

Lutsenko

,

M.

Drevlak

,

Y.

Turkin

, and

P.

Helander

,

Nucl. Fusion

56

,

066004

(

2016

).

https://doi.org/10.1088/0029-5515/56/6/066004

Google Scholar

Crossref

39.

I.

Novikau

,

A.

Biancalani

,

A.

Bottino

,

A.

Di Siena

,

P.

Lauber

,

E.

Poli

,

E.

Lanti

,

L.

Villard

,

N.

Ohana

, and

S.

Briguglio

, “

Implementation of energy transfer technique in ORB5 to study collisionless wave-particle interactions in phase-space

,”

Comput. Phys. Commun.

(to be published). arXiv:1903.05024 [physics.plasm-ph].

40.

A.

Di Siena

,

A.

Biancalani

,

T.

Görler

,

H.

Doerk

,

I.

Novikau

,

P.

Lauber

,

A.

Bottino

, and

E.

Poli

,

Nucl. Fusion

58

,

106014

(

2018

).

https://doi.org/10.1088/1741-4326/aad51d

Google Scholar

Crossref

41.

P.

Lauber

, “

Off-axis NBI heated discharges at ASDEX Upgrade: EGAMs, RSAEs, TAE bursts

,” in

13th Energetic Particle Physics TG Meeting

(

2014

).

Google Scholar

42.

M. V.

Falessi

,

N.

Carlevaro

,

V.

Fusco

,

G.

Vlad

, and

F.

Zonca

,

Phys. Plasmas

26

,

082502

(

2019

).

https://doi.org/10.1063/1.5098982

Google Scholar

Crossref

43.

G.

Fogaccia

,

G.

Vlad

, and

S.

Briguglio

,

Nucl. Fusion

56

,

112004

(

2016

).

https://doi.org/10.1088/0029-5515/56/11/112004

Google Scholar

Crossref

2020

Author(s)

Gyrokinetic investigation of the damping channels of Alfvén modes in ASDEX Upgrade

I. INTRODUCTION

II. MODEL

III. PHASE MIXING/CONTINUUM DAMPING

IV. LANDAU DAMPING

V. NLED-AUG CASE

A. Linear TAEs

B. Linear EPM

C. Nonlinear simulations

VI. CONCLUSION

ACKNOWLEDGMENTS

References

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

Gyrokinetic investigation of the damping channels of Alfvén modes in ASDEX Upgrade

I. INTRODUCTION

II. MODEL

III. PHASE MIXING/CONTINUUM DAMPING

IV. LANDAU DAMPING

V. NLED-AUG CASE

A. Linear TAEs

B. Linear EPM

C. Nonlinear simulations

VI. CONCLUSION

ACKNOWLEDGMENTS

References

Citing articles via

Submit your article

Sign up for alerts

Related Content

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only