# Calculation of toroidal Alfvén eigenmode mode structure in general
axisymmetric toroidal geometry

Guangyu Wei  Institute for Fusion Theory and Simulation and Department of
Physics, Zhejiang University, Hangzhou 310027, China  Center for Nonlinear
Plasma Science and C.R. ENEA Frascati, Via E. Fermi 45, Frascati, Italy
Matteo Valerio Falessi  Center for Nonlinear Plasma Science and C.R. ENEA
Frascati, Via E. Fermi 45, Frascati, Italy  Istituto Nazionale di Fisica
Nucleare (INFN), Sezione di Roma, Piazzale Aldo Moro 2, Roma, Italy  Tao Wang
Institute for Fusion Theory and Simulation and Department of Physics, Zhejiang
University, Hangzhou 310027, China  Center for Nonlinear Plasma Science and
C.R. ENEA Frascati, Via E. Fermi 45, Frascati, Italy  Fulvio Zonca
<EMAIL_ADDRESS>Center for Nonlinear Plasma Science and C.R. ENEA
Frascati, Via E. Fermi 45, Frascati, Italy  Institute for Fusion Theory and
Simulation and Department of Physics, Zhejiang University, Hangzhou 310027,
China  Zhiyong Qiu  Key Laboratory of Frontier Physics in Controlled Nuclear
Fusion and Institute of Plasma Physics, Chinese Academy of Sciences, Hefei
230031, China  Center for Nonlinear Plasma Science and C.R. ENEA Frascati, Via
E. Fermi 45, Frascati, Italy

###### Abstract

A workflow is developed based on the ideal MHD model to investigate the linear
physics of various Alfvén eigenmodes in general axisymmetric toroidal
geometry, by solving the coupled shear Alfvén wave (SAW) and ion sound wave
(ISW) equations in ballooning space. The model equations are solved by the
FALCON code in the singular layer, and the corresponding solutions are then
taken as the boundary conditions for calculating parallel mode structures in
the whole ballooning space. As an application of the code, the frequencies and
mode structures of toroidal Alfvén eigenmode (TAE) are calculated in the
reference equilibria of the Divertor Tokamak Test facility (DTT) with positive
and negative triangularities, respectively. By properly handling the boundary
conditions, we demonstrate finite TAE damping due to coupling with the local
acoustic continuum, and find that the damping rate is small for typical plasma
parameters.

††preprint: AIP/123-QED

##  I Introduction

In magnetically confined fusion plasmas, excitation of shear Alfvén wave (SAW)
instabilities by energetic particles (EPs) via wave-particle resonance has
been one of the major concerns due to their potential role in EP transport
Chen and Zonca (2016). Therefore, a comprehensive understanding of the
dynamics of SAW fluctuations in fusion plasmas is of crucial importance. It is
well known that, in toroidal plasma devices such as tokamaks, SAWs and ion
sound waves (ISWs) are characterized by continuously varying continuous
spectra. The periodic variation of equilibrium quantities along the magnetic
field lines leads to the breakup of the continuous spectrum and formation of
SAW-ISW frequency gaps Cheng, Chen, and Chance (1985); Cheng and Chance
(1986), where discrete eigenmodes may exist. The discrete Alfvén eigenmodes
(AEs), such as toroidal Alfvén eigenmode (TAE) Cheng, Chen, and Chance (1985);
Cheng and Chance (1986); Heidbrink  _et al._ (1991) and beta-induced Alfvén
eigenmode (BAE) Heidbrink  _et al._ (1993), which are inside the frequency
gaps and are free of significant continuum damping Rosenbluth  _et al._
(1992); Zonca and Chen (1992, 1993), can be easily driven unstable by EPs with
a low threshold condition and, thus, have attracted a lot of attention in
fusion research. If the EP drive is sufficiently strong to overcome the
continuum damping, it can give rise to another type of discrete fluctuations
known as energetic particle modes (EPMs) Chen (1994) originated from the SAW
continuous spectrum. Therefore, it is particularly important to precisely
calculate the structures of SAW continuous spectrum and the frequencies of
various AEs in realistic magnetic configurations with plasma non-uniformities
included, which can serve as a valuable reference for analyzing and
interpreting simulation results and experimental observations.

The coupled SAW-ISW continuous spectrum can be obtained by solving the
equations that describe the fluctuations propagation along the equilibrium
magnetic field lines Chu  _et al._ (1992); Falessi  _et al._ (2019, 2020).
This is equivalent to the typical approach that takes a spectral decomposition
of the fluctuations in poloidal and toroidal angles and computes the null
space (kernel) of the matrix of the highest order radial derivative Cheng and
Chance (1986), corresponding to radially singular SAW-ISW structures. It is
worthy noting that, the former method is convenient when dealing with high
toroidal mode numbers, and can be straightforwardly extended to kinetic
descriptionFalessi  _et al._ (2019). Besides, by solving the differential
equations of SAW-ISW continuous spectrum in ballooning space, the obtained
solutions also determine the asymptotic behaviors of the discrete AEs in the
frequency gaps Chen and Zonca (2017); Falessi  _et al._ (2019).

In previous worksFalessi  _et al._ (2019, 2020); Pucella  _et al._ (2022),
some of the authors developed the Floquet Alfvén continuum code (FALCON) to
calculate the coupled SAW-ISW continuous spectrum in general axisymmetric
toroidal geometry. FALCON, along with its associated library called
EQUIlibrium Postrpocessing codE (EQUIPE), are open source tools avaible to the
fusion community. They have been collaboratively developed111The access to the
code repository can be obtained by writing to the code maintainer at
<EMAIL_ADDRESS>validatedPucella  _et al._ (2022); Piron  _et al._
(2023) and verified in the last years. In this work, we are dedicated to
develop a workflow for the calculation of AEs within this framework and adopt
the ideal MHD limit to illustrate the general methodology in a simple case.
The extension of present calculation to kinetic description is straightforward
thanks to the adoption of ballooning mode representation Connor, Hastie, and
Taylor (1978); Lu, Zonca, and Cardinali (2012), and will be pursued in future
work. The local dispersion curves calculated by FALCONFalessi  _et al._
(2019); Pucella  _et al._ (2022) can be used to construct the global
continuous spectrum. Meanwhile, they can be taken as the boundary condition
for calculating the parallel mode structures of discrete gap modes. Therefore,
we have developed a general workflow, based on the FALCON code, to calculate
the frequencies and mode structures of AEs. As an application of the code, in
this work, we calculate the TAE Cheng, Chen, and Chance (1985); Cheng and
Chance (1986) spectrum in one of the reference equilibria of the Divertor
Tokamak Test facility (DTT)Albanese and Pizzuto (2017), and demonstrate the
existence of a small but finite TAE damping rate in the ideal MHD limit. This
damping results from the coupling between TAE and the acoustic continuum, so
it is generally much smaller than the conventional SAW continuum damping for
typical plasma parameters, due to the different polarizations of the eigenmode
and the local ISW continuum. The dependence of this effect on plasma shaping
is illustrated here as we calculate TAE in a DTT equilibrium with negative
triangularity (NT), demonstrating the general applicability of the workflow
across different geometries. We find minimal deviation from previous results
with positive triangularity (PT), which is consistent with the experimental
observations on DIII-D tokamakZeeland  _et al._ (2019).

The paper is organized as follows. In Sec. II, we introduce the model
equations and numerical approach adopted in the code. In Sec. III, we apply
the code to analyze the frequencies and parallel mode structure of TAE in DTT
equilibria with positive and negative triangularities and discuss the weak
effect of triangularity on TAE. In Sec. IV, we give a brief summary of current
work and outline the future prospects for code development.

##  II Theoretical Framework

###  II.1 Model equations

Following Ref. 13, we introduce our model equations that describe the SAW and
ISW fluctuations coupled by the magnetic curvature, which consist of the
vorticity equation

$\displaystyle
B_{0}\nabla_{\parallel}\left(\frac{1}{B_{0}}\nabla_{\perp}^{2}\nabla_{\parallel}\Phi_{s}\right)-\bm{\nabla}_{\perp}\cdot\left(\frac{\partial_{t}^{2}}{v_{A}^{2}}\bm{\nabla}_{\perp}\Phi_{s}\right)$
(1)
$\displaystyle+\frac{8\pi}{B_{0}}(\bm{b}_{0}\times\bm{\kappa})\times\bm{\nabla}_{\perp}\left[\frac{\bm{b}_{0}}{B_{0}}\cdot(\bm{\nabla}P_{0}\times\bm{\nabla}_{\perp}\Phi_{s})\right]$
$\displaystyle+\frac{4\pi}{cB_{0}}\left[\bm{b}_{0}\times\bm{\nabla}_{\perp}(\nabla_{\parallel}\Phi_{s})\right]\cdot\bm{\nabla}J_{0\parallel}$
$\displaystyle=-\frac{8\pi}{cB_{0}}(\bm{b}_{0}\times\bm{\kappa})\cdot\bm{\nabla}_{\perp}\delta
P_{\mathrm{comp}},$

and the parallel force balance equation

$\displaystyle\partial_{t}^{2}\delta
P_{\mathrm{comp}}-c_{S}^{2}B_{0}\nabla_{\parallel}\left(\frac{1}{B_{0}}\nabla_{\parallel}\delta
P_{\mathrm{comp}}\right)$ (2) $\displaystyle=-\frac{2\Gamma
P_{0}c}{B_{0}}(\bm{b}_{0}\times\bm{\kappa})\cdot\bm{\nabla}_{\perp}\partial_{t}^{2}\Phi_{s},$

where $\Phi_{s}$ is the perturbed stream function, related with the
perpendicular plasma displacement
$\delta\bm{\xi}_{\perp}=(c/B_{0})\bm{b}_{0}\times{\bm{\nabla}\Phi_{s}}$,
$\delta P_{\mathrm{comp}}=-\Gamma P_{0}(\bm{\nabla}\cdot\bm{\xi})$ is the
compressional component of the perturbed pressure,
$\bm{\kappa}=\bm{b}_{0}\cdot\bm{\nabla}\bm{b}_{0}$ is the magnetic curvature,
$v_{A}^{2}=B_{0}^{2}/(4\pi\rho_{0})$ and $c_{S}^{2}=\Gamma P_{0}/\rho_{0}$ are
Alfvén and sound speed, respectively, and other notations are standard. In Eq.
(1), the first three terms on the left-hand side represent field line bending,
plasma inertia and curvature-pressure coupling, respectively, while the fourth
term, representing the well-known kink drive, will be neglected in the
following study, since it is typically small for high mode number AEs
characterizing DTT plasmas as expected for reactor relevant conditionsWang
_et al._ (2018, 2019). The right-hand side is the contribution of curvature
coupling with plasma compressibility. In Eq. (2), the left-hand side can be
regarded as the linear response of ISW in a straight equilibrium magnetic
field, and the right-hand side is the correction of finite magnetic curvature.

In this work, we consider a general tokamak geometry with the equilibrium
magnetic field expressed by

$\bm{B}_{0}=F(\psi)\bm{\nabla}\varphi+\bm{\nabla}\varphi\times\bm{\nabla}\psi,$
(3)

where $\psi$ is the poloidal magnetic flux function, and $\varphi$ is the
physical toroidal angle. Following our previous worksFalessi  _et al._ (2019),
we make use of Boozer coordinates $(r,\theta,\zeta)$, where $r(\psi)$ is a
radial-like flux coordinate, poloidal angle $\theta$ is chosen such that
$\mathcal{J}B_{0}^{2}$ depends only on $\psi$, with
$\mathcal{J}=(\bm{\nabla}\psi\times\bm{\nabla}\theta\cdot\bm{\nabla}\zeta)^{-1}$
being the Jacobian of the coordinates, and the generalized toroidal angle
$\zeta$ is defined such that the magnetic field lines are straight on a flux
surface, i.e.,

$q=\frac{\bm{B}_{0}\cdot\bm{\nabla}\zeta}{\bm{B}_{0}\cdot\bm{\nabla}\theta}$
(4)

is a flux function. Focusing on modes with high toroidal mode number and
finite magnetic shear, we adopt the ballooning mode representation Connor,
Hastie, and Taylor (1978); Lu, Zonca, and Cardinali (2012), where an arbitrary
fluctuation field can be decomposed as

$\displaystyle f(r,\theta,\zeta)$
$\displaystyle=\sum_{m\in\mathbb{Z}}A_{n}(r)e^{i(n\zeta-m\theta)}\int
d\vartheta e^{i(m-nq)\vartheta}\hat{f}_{n}(\vartheta,r)$ (5)
$\displaystyle=2\pi A_{n}(r)\sum_{\ell\in\mathbb{Z}}e^{in\zeta-
inq(\theta-2\pi\ell)}\hat{f}_{n}(\theta-2\pi\ell,r),$

and the time dependence of the fluctuation field ($\propto e^{-i\omega t}$) is
left implicit for the sake of simplicity in notation. On the right hand side
of Eq. (5), $A(r)$ is the radial envelope of the mode, and
$\hat{f}_{n}(\vartheta,r)$ is the corresponding parallel mode structure in the
ballooning space, with $\vartheta$ being the "extended poloidal angle" along
the magnetic field line. By writing Eq. (5), we have assumed the scale
separation between radial envelope and single poloidal harmonics, which
enables us to reduce the original two-dimensional problem in the poloidal
plane to two one-dimensional problems Zonca and Chen (1993), i.e., solving for
the parallel mode structure $\hat{f}_{n}(\vartheta,r)$ in ballooning space at
the lowest order and for the radial envelope $A(r)$ at the next step. In this
work, we focus on the solution of the parallel mode structure.

Applying the ballooning mode representation and ignoring the effect of radial
envelope variation, the coupled SAW-ISW equations Eq. (1) and Eq. (2) become
Falessi  _et al._ (2019)

$\displaystyle\left[\partial_{\vartheta}^{2}-\frac{\partial_{\vartheta}^{2}\hat{\kappa}_{\perp}}{\hat{\kappa}_{\perp}}+\frac{\omega^{2}\mathcal{J}^{2}B_{0}^{2}}{v_{A}^{2}}-8\pi\mathcal{J}^{2}\frac{rB_{0}P_{0}^{\prime}}{q\hat{\kappa}_{\perp}\psi^{\prime}}\left(\kappa_{g}\frac{\bm{\nabla}\psi\cdot\hat{\bm{\kappa}}_{\perp}}{\hat{\kappa}_{\perp}|\bm{\nabla}\psi|}-\kappa_{n}\frac{rB_{0}}{q\hat{\kappa}_{\perp}|\bm{\nabla}\psi|}\right)\right]g_{1}(\vartheta)$
(6a)
$\displaystyle=-\left(2\Gamma\beta\right)^{1/2}\frac{\mathcal{J}^{2}B_{0}^{2}}{qR_{0}}\left(\kappa_{g}\frac{\bm{\nabla}\psi\cdot\hat{\bm{\kappa}}_{\perp}}{\hat{\kappa}_{\perp}|\bm{\nabla}\psi|}-\kappa_{n}\frac{rB_{0}}{q\hat{\kappa}_{\perp}|\bm{\nabla}\psi|}\right)g_{2}(\vartheta),$
$\left(\frac{\mathcal{J}^{2}B_{0}^{2}}{q^{2}R_{0}^{2}}+\frac{\omega_{S}^{2}}{\omega^{2}}\partial_{\vartheta}^{2}\right)g_{2}(\vartheta)=-\left(2\Gamma\beta\right)^{1/2}\frac{\mathcal{J}^{2}B_{0}^{2}}{qR_{0}}\left(\kappa_{g}\frac{\bm{\nabla}\psi\cdot\hat{\bm{\kappa}}_{\perp}}{\hat{\kappa}_{\perp}|\bm{\nabla}\psi|}-\kappa_{n}\frac{rB_{0}}{q\hat{\kappa}_{\perp}|\bm{\nabla}\psi|}\right)g_{1}(\vartheta),$
(6b)

where $g_{1}(\vartheta)$ and $g_{2}(\vartheta)$ are the projections of
$\Phi_{s}$ and $\delta P_{\mathrm{comp}}$ in ballooning space with convenient
normalization, $\beta$ is the ratio of kinetic and magnetic pressure,
$\omega_{S}$ is the sound frequency, $\kappa_{g}$ and $\kappa_{n}$ are
geodesic and normal curvature, and $\hat{\bm{\kappa}}_{\perp}$ is the
normalized perpendicular wave vector. The definitions of these quantities are
as follows:

$\displaystyle
g_{1}(\vartheta)=\frac{\hat{\kappa}_{\perp}\hat{\Phi}_{s}(\vartheta)}{\left(\beta
q^{2}\right)^{1/2}}\frac{ck_{\vartheta}}{B_{0}R_{0}},\quad
g_{2}(\vartheta)=\frac{i\delta\hat{P}_{\mathrm{comp}}(\vartheta)}{\left(2\Gamma\right)^{1/2}P_{0}},$
$\displaystyle\beta=\frac{8\pi
P_{0}}{B_{0}^{2}},\quad\omega_{S}^{2}=\frac{c_{S}^{2}}{q^{2}R_{0}^{2}},\quad\kappa_{g}=\frac{F(\psi)}{|\bm{\nabla}\psi|}\frac{\partial_{\theta}B_{0}}{\mathcal{J}B_{0}^{2}},$
$\displaystyle\kappa_{n}=\frac{|\bm{\nabla}\psi|}{B_{0}}\left(\partial_{\psi}B_{0}+\frac{4\pi}{B_{0}}\partial_{\psi}P_{0}\right)+\frac{\bm{\nabla}\psi\cdot\bm{\nabla}\theta}{|\bm{\nabla}\psi|}\frac{\partial_{\theta}B_{0}}{B_{0}},$
$\displaystyle\hat{\bm{\kappa}}_{\perp}=\frac{\bm{k}_{\perp}}{k_{\vartheta}}=s\vartheta\bm{\nabla}r+r\bm{\nabla}\theta-\frac{r}{q}\bm{\nabla}\zeta,$

where $k_{\vartheta}=-nq/r$, and $s=rq^{\prime}/q$ is the magnetic shear. For
simplicity, here, we have omitted the slow equilibrium radial scale dependence
in the functions $g_{1}$ and g2, noting that it is parameterized by the radial
variation of equilibrium quantities. Eqs. (6), describing the SAW-ISW mode
structures propagating along the magnetic field lines, are the starting point
of the present work.

Eqs. (6) with proper boundary conditions define a well posed eigenvalue
problem that generally doesn’t allow for analytical solutions except in some
particular cases with limiting parameters and simplified geometry Zonca and
Chen (1993); Chen and Zonca (2017). In this work, we developed a workflow to
solve Eqs. (6) numerically with proper boundary conditions. Compared with
other codes that directly solve the ideal MHD equations in real space,
adopting the ballooning mode representation in our model equations allows us
not only to compute the high-n modes with higher resolution, but also to study
the underlying physics in more detail, as will be shown below. In the
following, we will briefly introduce our methodology for solving this problem,
which was firstly proposed in Ref. 11.

###  II.2 Solving method

Taking the $|\vartheta|\rightarrow\infty$ limit, we can reduce Eqs. (6) to
Falessi  _et al._ (2019)

$\displaystyle\left(\partial_{\vartheta}^{2}-\frac{\partial_{\vartheta}^{2}|\bm{\nabla}r|}{|\bm{\nabla}r|}+\frac{\omega^{2}\mathcal{J}^{2}B_{0}^{2}}{v_{A}^{2}}\right)g_{1}$
(7a)
$\displaystyle=-\left(2\Gamma\beta\right)^{1/2}\frac{\mathcal{J}^{2}B_{0}^{2}}{qR_{0}}\frac{s\vartheta}{|s\vartheta|}\kappa_{g}g_{2},$
$\left(\frac{\mathcal{J}^{2}B_{0}^{2}}{q^{2}R_{0}^{2}}+\frac{\omega_{S}^{2}}{\omega^{2}}\partial_{\vartheta}^{2}\right)g_{2}=-\left(2\Gamma\beta\right)^{1/2}\frac{\mathcal{J}^{2}B_{0}^{2}}{qR_{0}}\frac{s\vartheta}{|s\vartheta|}\kappa_{g}g_{1}.$
(7b)

Notice that in Eq. (5), due to finite magnetic shear, the extended poloidal
angle $\vartheta$ also plays the role of a dimensionless radial wave vector.
So Eqs. (7), obtained from Eqs. (6) in $|\vartheta|\rightarrow\infty$ limit,
can actually describe the radial singular structures of plasma response, i.e.,
the continuous spectrum Chen and Zonca (2017); Falessi  _et al._ (2019, 2020).
More generally, by solving Eqs. (7), we obtain the asymptotic behaviors of the
parallel mode structures in the radial singular layer for both continuous
spectrum as well as AEs Zonca and Chen (2014a, b). These solutions will then
be taken as the boundary condition for Eqs. (6) to calculate parallel mode
structure in the whole ballooning space and the corresponding dispersion
relation. Here, we will consider the $\theta_{k}\equiv
k_{r}/(nq^{\prime})\rightarrow 0$ limit corresponding to the most unstable
mode radial envelope structure, with $k_{r}$ the radial envelope wave number.
At the same time, the radial envelope for the ground state of a radially
localized mode structure can be computed analytically. Here, this will provide
a simple ( "ad-hoc" ) example, which we will use to illustrate the two-
dimensional radial mode structure that can be reconstructed following the
present approach.

Following Ref. 11, we note that all the coefficients in Eqs. (7), related with
the equilibrium quantities, are $2\pi$-periodic functions of $\vartheta$; so,
Floquet theory Floquet (1883) can be employed to investigate the solutions of
this system. In this paper, we review only the most important results of
Floquet theory, and interested readers are referred to Refs. 25 and 11 for
more detailed discussions. Rewriting Eqs. (7) as four first order coupled
differential equations, we know from Floquet theory that, given the value of
$\omega$, they must have solutions in the form of

$\mathbf{x}_{i}(\omega;\vartheta)=\mathbf{P}_{i}(\vartheta)e^{i\nu_{i}\vartheta},\quad
i=1,2,3,4,$ (8)

where $\mathbf{x}_{i}$ is a column vector with four elements corresponding to
the solutions of $g_{1}$, $\partial_{\vartheta}g_{1}$, $g_{2}$ and
$\partial_{\vartheta}g_{2}$, respectively. $\mathbf{P}_{i}(\vartheta)$ is a
$2\pi$-periodic function, and $\nu_{i}$ is generally a complex number called
characteristic Floquet exponent. One of the major advantages of Floquet
solutions is that their asymptotic behaviors are totally determined by the
characteristic Floquet exponent $\nu_{i}$, so we can choose the solutions that
satisfy the physical boundary condition (outgoing/decaying wave Zonca and Chen
(2014a, b) in $\vartheta$ space), with only the information of
$\nu_{i}(\omega,r)$. Specifically, solutions in the continuum are traveling
waves that propagate toward $|\vartheta|\rightarrow\infty$, so the Floquet
exponent must satisfy $\mathrm{Im}\nu=0$ and
$\mathrm{sgn}(\vartheta)\partial\mathrm{Re}\nu/\partial\omega>0$. Solutions in
the frequency gap are standing waves which exponentially decay as
$|\vartheta|\rightarrow\infty$, corresponding to AEs, so the Floquet exponent
must satisfy $\mathrm{sgn}(\vartheta)\mathrm{Im}\nu>0$. The FALCON code,
developed by Falessi et al, calculates the local dispersion curves
$\nu(\omega,r)$ at each radial position and uses
$\nu^{2}(\omega,r)=(nq(r)-m)^{2}$ to compute the global continuous spectrum
Falessi  _et al._ (2019). The advantage and accuracy of this approach is due
to the fact that the dispersion curves $\nu(\omega,r)$ are independent of $n$
and of any fast radial variation. In this work, we use the dispersion curves
calculated by FALCON to determine the frequencies and parallel mode structures
of AEs.

Since for generally up-down asymmetric equilibrium, the solutions of Eqs. (6)
are asymmetric for parity transformation in $\vartheta$ space, we should
calculate the boundary condition on the left ($\vartheta\rightarrow-\infty$)
and right ($\vartheta\rightarrow+\infty$) separately. At each boundary, among
the four solutions in Eq. (8), generally only two of them satisfy the physical
boundary conditions. So we have the general solutions of Eqs. (7), for
negative and positive $\vartheta$, expressed by

$\displaystyle\mathbf{x}^{-}(\omega;\vartheta)=w_{1}^{-}\mathbf{x}_{1}^{-}+w_{2}^{-}\mathbf{x}_{2}^{-},$
(9)
$\displaystyle\mathbf{x}^{+}(\omega;\vartheta)=w_{1}^{+}\mathbf{x}_{1}^{+}+w_{2}^{+}\mathbf{x}_{2}^{+},$

respectively, where $\mathbf{x}_{1,2}^{-}$ and $\mathbf{x}_{1,2}^{+}$ are
selected from $\mathbf{x}_{i}$ in Eq. (8) taking into account the boundary
conditions, while $w_{1,2}^{-}$ and $w_{1,2}^{+}$ are weight coefficients of
the linear combination. Since Eqs. (7) is the limiting form of Eqs. (6) as
$|\vartheta|\rightarrow\infty$, it is possible to connect the solutions of
Eqs. (6), denoted by $\mathbf{X}^{-}(\vartheta)$ (for negative $\vartheta$)
and $\mathbf{X}^{+}(\vartheta)$ (for positive $\vartheta$), with the solutions
of Eqs. (7): $\mathbf{X}^{-}(-2\pi\ell)=\mathbf{x}^{-}(\omega;0)$ and
$\mathbf{X}^{+}(2\pi\ell)=\mathbf{x}^{+}(\omega;0)$, where $\ell$ is a
sufficiently large positive integer. The solution in the whole ballooning
space can be formally expressed by

$\displaystyle\mathbf{X}^{-}(\vartheta)=w_{1}^{-}\mathbf{X}_{1}^{-}+w_{2}^{-}\mathbf{X}_{2}^{-},$
(10)
$\displaystyle\mathbf{X}^{+}(\vartheta)=w_{1}^{+}\mathbf{X}_{1}^{+}+w_{2}^{+}\mathbf{X}_{2}^{+},$

where $\mathbf{X}_{1,2}^{-}$ and $\mathbf{X}_{1,2}^{+}$ are obtained by
numerically integrating Eqs. (6) from $-2\pi\ell$ and $2\pi\ell$ with boundary
condition $\mathbf{X}_{1,2}^{-}(-2\pi\ell)=\mathbf{x}_{1,2}^{-}(\omega;0)$ and
$\mathbf{X}_{1,2}^{+}(2\pi\ell)=\mathbf{x}_{1,2}^{+}(\omega;0)$, respectively.
The two solutions given in Eqs. (10) must match at $\vartheta=0$, i.e.,
$\mathbf{X}^{-}(0)=\mathbf{X}^{+}(0)$, which yields a linear matrix equation

$\mathbf{M}(\omega)\mathbf{w}=\mathbf{0},$ (11)

where $\mathbf{w}=(w_{1}^{-},w_{2}^{-},w_{1}^{+},w_{2}^{+})$ is a column
vector composed of the four weight coefficients, and $\mathbf{M}$ is a
$4\times 4$ matrix given by

$\begin{pmatrix}\mathbf{X}_{1}^{-}(0)&\mathbf{X}_{2}^{-}(0)&-\mathbf{X}_{1}^{+}(0)&-\mathbf{X}_{2}^{+}(0)\end{pmatrix}.$
(12)

Condition for Eq. (11) having nontrivial solution leads to the local
dispersion relation $D(\omega)\equiv\mathrm{det}(\mathbf{M})=0$, where we have
omitted its explicit dependence on the local parameters for simplicity. In
practice, the eigenvalue $\omega$ is found by scanning the frequency gap and
minimizing the value of $|D(\omega)|$. The corresponding solution of
$\mathbf{w}$ is then substituted into Eqs. (10) to construct the parallel mode
structure.

##  III Numerical results

The model equations and numerical approach introduced in previous section for
solving SAW-ISW eigenmodes is general. In this section, we will illustrate,
taking the well-known TAE as an example, the calculation of mode frequency and
parallel mode structure in general geometry of experimental interest. For
consistency with Ref. 11, we use the equilibrium of a DTT reference scenario
Albanese and Pizzuto (2017) to illustrate the application of our code. The
magnetic equilibrium has been originally calculated by means of the free
boundary equilibrium evolution code CREATE-NL Albanese  _et al._ (2003) and
further refined using the high-resolution equilibrium solver CHEASE Lütjens,
Bondeson, and Sauter (1996). We consider a single null configuration with
magnetic axis and plasma edge excluded to avoid the possible singularity, as
shown in FIG. 1. The profiles of safety factor and magnetic shear are depicted
in FIG. 2 (a), with $\rho_{tor}$ the normalized toroidal radiusFalessi  _et
al._ (2019). Profile of $\bar{\beta}\equiv 8\pi P_{0}/\bar{B}_{0}^{2}$ is
depicted in FIG. 2 (b), with $\bar{B}_{0}$ the magnetic field on axis.

Figure 1: Contour lines of $\psi$ and $\theta$ coordinates.

Figure 2: (a) Profiles of safety factor $q$ and magnetic shear $s$. (b)
Profile of $\bar{\beta}$.

###  III.1 Boundary condition

Figure 3: (a) Real parts and (b) imaginary parts of the characteristic Floquet
exponents $\nu_{i}(\Omega)$.

Following the workflow described in Sec. II.2, we firstly calculate all the
coefficients of Eq. (6) using EQUIPE. Then we call FALCON to compute the
dispersion curves and the boundary condition of Eqs. (6). FIG. 3 (a) and (b)
show the real and imaginary parts of $\nu_{i}(\Omega)$ at $\rho_{tor}=0.54$,
respectively. Here, $\Omega=\omega R_{0}/\bar{v}_{A0}$ with $\bar{v}_{A0}$ the
Alfvén velocity on the magnetic axis, and for convenience, $\Omega$ is chosen
as the vertical coordinates to show different branches and gaps more clearly.
Then we can pick out the SAW-ISW continuum with $\mathrm{Im}{\nu}=0$, as
depicted in FIG. 4 (a), where the shade of the dot’s color is proportional to
the value of Alfvénicity Falessi  _et al._ (2019, 2020); Pucella  _et al._
(2022)

$\mathcal{A}=\frac{\int_{0}^{2\pi}g_{1}^{2}(\vartheta;\nu,r)d\vartheta}{\int_{0}^{2\pi}[g_{1}^{2}(\vartheta;\nu,r)+g_{2}^{2}(\vartheta;\nu,r)]d\vartheta},$
(13)

which characterizes the degree of SAW polarization of the fluctuation. Here,
$g_{1}$ and $g_{2}$ denote the solutions of Eqs. (7). According to FIG. 4 (a),
there is clearly a SAW polarized frequency gap at $\Omega\sim 0.5$, which is
caused by the interaction of two counter-propagating SAWs along the magnetic
field line. We know that TAE is located in this frequency gap, and the
existence of acoustic polarized continuum in the same region implies the
possible damping of TAE due to the coupling with ISW continuum, as will be
shown in the next subsection. Some insights into the coupling of AEs and ISW
continuum can be obtained by checking the properties of the characteristic
Floquet exponents in FIG. 3 (a) and (b). Within the range
$0.38\lesssim\Omega\lesssim 0.68$, $\nu_{2}=-\nu_{3}$ is purely a real number,
which is related to the ISW continuum solution, while $\nu_{1}=-\nu_{4}$ has a
real part of $-0.5$ and a negative imaginary part, which is related to the SAW
gap AE solution. The asymptotic solution of TAE parallel mode structure in
this frequency range is generally given by the mixture of both SAW gap and ISW
continuum solutions

$\displaystyle\bm{x}^{-}(\omega;\vartheta)=w_{1}^{-}P_{1}(\vartheta)e^{i\nu_{1}\vartheta}+w_{2}^{-}P_{2}(\vartheta)e^{i\nu_{2}\vartheta},$
(14)
$\displaystyle\bm{x}^{+}(\omega;\vartheta)=w_{1}^{+}P_{3}(\vartheta)e^{i\nu_{3}\vartheta}+w_{2}^{+}P_{4}(\vartheta)e^{i\nu_{4}\vartheta},$

where outgoing wave boundary conditions have been imposed, and $\omega$ is
assumed to be located at a position where
$\partial\nu_{2}/\partial\omega=-\partial\nu_{3}/\partial\omega<0$, otherwise,
subscripts $2$ and $3$ in Eqs. (14) should be exchanged.

###  III.2 Frequency and parallel mode structure of TAE in DTT

In this subsection, we calculate the frequency and parallel mode structure of
TAE in the DTT reference equilibrium. FIG. 4 (b) depicts the SAW-ISW
continuous spectrum of $n=5$ calculated by FALCON in this equilibrium, with
Alfvénicity denoted by the color-bar. Given the boundary condition in Sec.
III.1, and adopting the workflow introduced in Sec. II.2, we calculate the TAE
frequency and parallel mode structure at the radial position
$\rho_{tor}=0.54$, where $\bar{\beta}=2.28\%$. The real frequency of TAE is
$\Omega_{r}=0.649$, close to the upper SAW continuum, as marked by the red
plus "+" in FIG. 4. Meanwhile, coupling of TAE with the local ISW continuum
will lead to a finite damping rate of the mode, due to the generation of short
scale radial singular structure. Considering the different polarization of the
AE and the resonantly excited ISW continuum (FIG. 4 (a)), the damping rate of
TAE $\Omega_{i}$ is expected to be very small. This justifies the use of a
perturbation method to calculate $\Omega_{i}$. Upon determining the value of
$\Omega_{r}$ by finding the roots of $D_{r}(\Omega_{r})$, we then obtain the
value of $\Omega_{i}$ as

$\Omega_{i}=-\frac{D_{i}}{\partial
D_{r}/\partial\Omega}\Big{|}_{\Omega=\Omega_{r}},$ (15)

where $D_{r}$ and $D_{i}$ are real and imaginary parts of $D(\Omega_{r})$,
respectively. In this way, the obtained value of $\Omega_{i}=-1.61\times
10^{-4}$, which is much smaller than $\Omega_{r}$ as expected.

Figure 4: (a) Local dispersion curves $\nu(\Omega)$ at $\rho_{tor}=0.54$, with
the color-bar showing the Alfvénicity defined by Eq. (13), whose value is
connected with the polarization of the corresponding fluctuationFalessi  _et
al._ (2019, 2020). (b) SAW-ISW continuous spectrum of $n=5$ calculated by
FALCON in DTT, and the color-bar shows the Alfvénicity. The red plus "+" marks
the frequency of TAE calculated at $\rho_{tor}=0.54$.

The small but finite TAE coupling with ISW continuum and the consequent
damping rate is consistent with the parallel mode structures of $g_{1}$ and
$g_{2}$ in ballooning space, as shown in FIG. 5. $g_{1}$, representing the
response of SAW, is characterized by an exponentially decaying behavior
corresponding the formation of a gap mode. Conversely, $g_{2}$, which
represents the response of ISW, exhibits a periodic oscillating behavior
typical of a continuum mode. The oscillation of $g_{2}$ is much faster than
$g_{1}$, which can be explained by the fact that an Alfvén wave with
relatively low $k_{\parallel}$ is interacting with sound waves with higher
$k_{\parallel}$. The parallel mode structure in large-$\vartheta$ limit
implies that there is an energy flux from large scale to small scale. This
physical picture can be seen more clearly using a variational principle Qiu,
Zonca, and Chen (2012). Specifically, we perform the following operation on
Eqs. (6)

$g_{1}^{*}\cdot\mathrm{(6a)}-g_{1}\cdot\mathrm{(6a)}^{*}+\left[g_{2}^{*}\cdot\mathrm{(6b)}-g_{2}\cdot\mathrm{(6b)}^{*}\right],$

and integrate over the ballooning space. Thus, we get

$\displaystyle
4i\omega_{r}\omega_{i}\int_{-\infty}^{\infty}\frac{\mathcal{J}^{2}B_{0}^{2}}{v_{A}^{2}}|g_{1}|^{2}d\vartheta\simeq-\left(g_{1}^{*}\partial_{\vartheta}g_{1}-g_{1}\partial_{\vartheta}g_{1}^{*}\right)\Big{|}_{-\infty}^{\infty}$
(16)
$\displaystyle-\frac{\omega_{S}^{2}}{\omega_{r}^{2}}\left(g_{2}^{*}\partial_{\vartheta}g_{2}-g_{2}\partial_{\vartheta}g_{2}^{*}\right)\Big{|}_{-\infty}^{\infty},$

where "*" is the standard notation for complex conjugate. In Eq. (16), the
left-hand side represents the changing rate of the total energy, and the two
terms on the right-hand side represent the outgoing energy flux, contributed
by $g_{1}$ and $g_{2}$, respectively. With the parallel mode structures of
$g_{1}$ and $g_{2}$ already obtained numerically, we plot the energy fluxes of
$g_{1}$ and $g_{2}$ in FIG. 6. The combination of Eq. (16) and FIG. 6
demonstrates that the energy of the mode primarily originates from $g_{1}$,
whereas the energy flux is attributed to $g_{2}$, which is consistent with the
physical picture that the TAE in the gap is slightly modified by the acoustic
continuum. Furthermore, we calculate the damping rate directly given by Eq.
(16), resulting in a value of $\Omega_{i}=-1.84\times 10^{-4}$, which agrees
well with the previous result from Eq. (15), obtained from a perturbation
expansion.

The damping rate of TAE resulting from the coupling with acoustic continuum is
notably small, about $\mathcal{O}(10^{-4})$ of the real frequency. This
damping rate can easily be covered by numerical errors in most codes. Here,
the local analysis that we have adopted gives us the capability to investigate
such a subtle physical effect. However, we should note that, the value of the
damping rate is not always so small, as it depends on the strength of SAW-ISW
coupling and, ultimately, on the mode polarizations Falessi  _et al._ (2020).
It can be expected that, SAW and ISW may be more strongly coupled for higher
$\bar{\beta}$, and the damping rate could be more relevant. For example, at
$\rho_{tor}=0.385$, $\bar{\beta}=3.53\%$, the obtained TAE real frequency is
$\Omega_{r}=0.770$, with a damping rate of $\Omega_{i}=-1.43\times 10^{-3}$.
It is also worth mentioning that Eqs. (6) generally have multiple solutions
satisfying the outgoing wave boundary conditions in the complex $\Omega$
space. The solution we have reported is the Alfvén polarized one,
corresponding to the well-known TAE. Meanwhile, there also exist acoustic
polarized solutions in the same frequency range, which are highly damped due
to the strong coupling with the continuum and are, therefore, less relevant
and ignored here. We emphasize, again, that mode polarization is the control
physical quantity that determines how strongly damped the corresponding
fluctuation structure is (cf. the color shade in Fig. 4, reflecting the
AlfvénicityFalessi  _et al._ (2019, 2020); Pucella  _et al._ (2022)).

For a more intuitive understanding of the physics mechanism underlying the TAE
damping by resonant ISW continuum excitation, it is helpful to investigate the
mode structures in real space. The $m$-th harmonics of the normalized
$\Phi_{s}$ and $\delta P_{\mathrm{comp}}$ in real space, obtained by direct
Fourier transform of $g_{1}/\hat{\kappa}_{\perp}$ and $g_{2}$ consistent with
Eq. (5), are shown in FIG. 7. There are clearly two dominant radial singular
structures of $\delta P_{\mathrm{comp}}$ around $nq-m=\pm 4.5$, which reflect
result of acoustic continuum resonances where
$\omega^{2}-(nq-m)^{2}\omega_{S}^{2}=0$ is satisfied. Meanwhile, the sidebands
around $nq-m=\pm 3.5,\pm 5.5$ arise from the coupling of different poloidal
harmonics. By analogy with SAW resonant absorptionChen and Hasegawa (1974);
Hasegawa and Chen (1974), one may expect the radial singular structures
illustrated here bear acoustic continuum damping due to the same process. In
the present case, the multiple resonance absorption layers must be attributed
to the rich SAW-ISW coupling due to geometry effectsFalessi  _et al._ (2019,
2020); Pucella  _et al._ (2022). This makes the accurate calculation of the
resonant absorption quite involved in the real space and shows the better
accuracy of our present approach. Meanwhile, it is evident that the damping
rate is relatively small since the resonances occur in high $|k_{\parallel}|$
region, where the dominant TAE mode structure ($\Phi_{s}$) rapidly decays away
from its peak.

###  III.3 Global mode structure of TAE in DTT

In this study, we have omitted the slow radial variation of the parallel mode
structures and focused on the solution of the local problem, leaving the self-
consistent solution of the global problem to future studies. Nonetheless, as
anticipated in Sec. II, the radial envelope for the ground state of a radially
localized mode structure can be computed analytically Zonca and Chen (1993);
Li  _et al._ (2023) and yields $A(r)=\exp[-(r-r_{0})^{2}/\Delta r^{2}]$, where
$r_{0}$ is the reference rational surface, and

$\Delta
r^{2}=\frac{2}{|nq^{\prime}|}\left[\frac{\partial^{2}D/\partial\theta_{k}^{2}}{\partial^{2}D/\partial
r^{2}}\right]^{1/2}_{r=r_{0}}.$

Here, $D(\omega,r,\theta_{k})$ denotes the local TAE dispersion function.
Using this analytical ("ad-hoc") radial envelope, we can schematically
construct the global fluctuation fields in real space with the help of Eq.
(5). FIG. 8 (a) and FIG. 8 (b) show the two- and three-dimensional structures
of the perturbed stream function $\Phi_{s}$ with toroidal mode number $n=20$.
We note that FIG. 8 (a) is similar to the FIG. (7) in Ref. 31, although a more
accurate kinetic description is adopted in that work. Note that the radial
singular structures due to SAW-ISW coupling are not visible in FIG. 8 due to
the fact that the effect is finite but small, as discussed above. Actually,
this radial singular will be completely eliminated in kinetic description,
when the finite Larmor radius effect comes into play and modifies the
potential functions Li  _et al._ (2023) in the model equations. Detailed
studies of these physics effects will be reported elsewhere.

Figure 5: Normalized $g_{1}$ and $g_{2}$ in $\vartheta$ space. Figure 6:
Energy fluxes in $\vartheta$ space. The black and red lines represents the
energy flux contributed by $g_{1}$ and $g_{2}$, respectively. Figure 7: Mode
structures of the $m$-th harmonics of the normalized $\Phi_{s}$ and $\delta
P_{\mathrm{comp}}$ in real space. The horizontal axis $nq-m$ can be regarded
as the radial coordinate.

Figure 8: (a) Two-dimensional and (b) three-dimensional mode structures of a
$n=20$ TAE obtained with an "ad-hoc" radial envelope.

As a final remark, it is worth mentioning that, for the high frequency modes
such as TAE analyzed here, slow sound approximation (SSA) Chu  _et al._
(1992); Falessi  _et al._ (2019) is usually adopted to further simplify Eqs.
(6). The range of validity of this approximation has been systematically
studied in Ref. 12. In SSA, one drops the $\partial_{\vartheta}^{2}$ term in
Eq. (6b) and turn it into an algebraic equation, so Eqs. (6) are reduced to
one second order differential equation. By reapplying the methodology outlined
in Sec. II.2, we calculate TAE in SSA using the same parameters. The frequency
of TAE in SSA is $\Omega=0.645$, very close to the result in the full system.
The corresponding parallel mode structure is shown in FIG. 9, which indicates
that $g_{1}$ has little variation compared to that in the full system, while
$g_{2}$ is significantly modified. This is indeed expected, since ISW
resonances are removed in SSA, with the obvious consequence of removing TAE
damping due to corresponding ISW resonance absorption. The comparison above
between the results of SSA and the full system demonstrates that the former
encapsulates most of the significant characters of TAE, such as the frequency
and the parallel model structure of $g_{1}$, while the latter gives a finite
amplitude of the parallel mode structure in singular layer and the resultant
TAE damping rate. We note that, studying the fully compressible response will
be crucial for the investigation of low-frequency AEs in future works,
although rigorous calculation of the full system only introduces a
perturbative correction compared to that of SSA for the high-frequency mode
such as TAE analyzed here.

Figure 9: Normalized $g_{1}$ and $g_{2}$ in $\vartheta$ space, adopting slow
sound approximation.

###  III.4 Effect of NT on TAE

Recently, researches on using the NT plasma shaping as tokamak reactor design
are revitalized due to its potential advantage in achieving a better plasma
confinement while avoiding the threat of edge localized modes (ELMs) Camenen
_et al._ (2007); Austin  _et al._ (2019); Nelson  _et al._ (2023). However,
present researches on NT have not yet fully explored the spectrum of
$n\sim\mathcal{O}(10)$ AEs Wang  _et al._ (2018), whose impact on alpha
particle transport could play a crucial role in the burning plasma power
balance. In order to study the effect of different triangularities on TAEs, we
calculate TAE in one of the proposed DTT NT quilibria and compare the results
with those obtained previously in PT equilibrium. The cross section of the NT
equilibrium is depicted in FIG. (10), which is obtained by mirror flipping the
last magnetic surface in FIG. (1) and recalculating the equilibrium using
CHEASE code at fixed $q$ and pressure profiles, enabling us to isolate and
analyze the specific impact of triangularity.

FIG. (11) depicts the continuous spectra of PT and NT equilibria calculated by
EQUIPE-FALCON, where the slow sound approximation is adopted to focus on the
variation of the TAE gap. The lower branches of the continuous spectra in both
equilibria are situated within the same frequency range, while the upper
branch is slightly lower in the NT equilibrium. We calculate the TAE frequency
in this NT equilibrium, which is found to be $\Omega=0.616-1.67\times
10^{-4}i$, and also slightly lower than that of PT equilibrium
($\Omega=0.649-1.61\times 10^{-4}i$). The real frequencies of TAEs in these
two equilibria are marked by black and red pluses in FIG. 11, respectively,
and their relative difference is negligibly small. The parallel mode structure
of TAE is depicted in FIG. 12, which, as expected, shows a similar pattern
with before (FIG. 5). Therefore, we conclude that, in the ideal MHD limit,
simply changing the triangularity of the equilibrium does not have a
significant influence on the frequency or mode structure of TAE. This result
is not surprising, considering that TAE spectrum has a predominant feature at
$|\mathcal{J}B_{0}k_{\parallel}|\sim 1/2$. Thus, the TAE dispersion relation,
written in the form of a variational principleZonca and Chen (2014a, b),
readily shows that triangularity effects average out at the leading order.
This is also consistent with toroidal frequency gap formation being dominanted
by the $\cos\theta$-type variation of the magnetic field (i.e., toroidicity),
while triangularity merely introduces a small component of higher order
sidebands. In addition, since TAEs are radially localized deeper in the core
region as compared to ELMs, triangularity may also attenuate its effect as one
moves from plasma edge towards the magnetic axis. Our present results are
consistent with other analyses reported in the literaturesGhai  _et al._
(2021); Oyola et al. (2023) about the relatively weak impact of NT on AE
fluctuation structures. They also confirm that NT could instead be relevant in
the AE excitation mechanism by EP via modification of both resonance
conditions and EP orbitsGhai  _et al._ (2021); Oyola et al. (2023); Zonca  _et
al._ (2015). Adopting the workflow proposed here, a more detailed
investigation on this issue with kinetic effects included will be reported in
our future work.

Figure 10: Contour lines of $\psi$ and $\theta$ coordinates in NT equilibrium.
Figure 11: Continuous spectra of PT and NT equilibria, adopting slow sound
approximation.  Figure 12: Normalized $g_{1}$ and $g_{2}$ in $\vartheta$ space
in NT equilibrium.

##  IV Summary and prospect

In this work, we illustrate the numerical calculation of AEs in the ideal MHD
limit, adopting the general theoretical framework discussed in Refs. 23, 24,
11 and 12. By solving the model equations in the ballooning space, the
asymptotic behaviors of the solutions in the large-$|\vartheta|$ limit are
completely determined by the corresponding characteristic Floquet exponents,
which makes it straightforward to select the proper boundary conditions and,
thus, defines a well posed eigenvalue problem in the ballooning space. To
solve this problem efficiently, a workflow has been developed on the basis of
FALCONFalessi  _et al._ (2019, 2020), which can be used to calculate the
frequencies and parallel mode structures of AEs in realistic tokamak
geometries. As an application of the code, we calculate the TAE in the
reference equilibria of DTT and demonstrate it has a small but finite damping
rate due to the coupling with the acoustic continuum. The existence of the
damping rate is also confirmed by the derived variational principle, which
highlights that the acoustic polarized component of the fluctuation carries a
finite amount of energy leaking out of the potential well in ballooning space;
i.e., toward radial singular structures, resulting in the dissipation of the
mode. Although the obtained TAE damping rate for the parameters used here is
small, due to the different polarization of the eigenmode and the ISW
continuum, one should note that, in fusion plasma with higher $\beta$, SAW and
ISW may be more strongly coupled even in the TAE frequency range, and the
mechanism that we have proposed could potentially become more relevantFalessi
_et al._ (2020). By comparing the numerical results in different DTT
equilibria with PT and NT, we find that simply changing the triangularity does
not have a significant influence on TAE frequency and mode structure within
the ideal MHD description.

Although we have limited ourselves to the local problem by neglecting the slow
variation of the radial envelope in Eqs. (6) for the sake of simplicity, the
self-consistent solution of the global problem is included in the general
theoretical framework Zonca and Chen (2014a, b). The general approach adopted
in the present work can be straightforwardly extended to kinetic description
due to the similar structure of the mode equationsZonca and Chen (2014b),
allowing to explore the effects of finite Larmor radius, wave-particle
interaction and diamagnetic frequency on the SAW-ISW fluctuations, as
anticipated in Ref. 31.

## Acknowledgements

The authors acknowledge Prof. Liu Chen (Zhejiang University, China) for
valuable discussions. This work was supported by the National Science
Foundation of China under Grant Nos. 12275236 and 12261131622, and Italian
Ministry for Foreign Affairs and International Cooperation Project under Grant
No. CN23GR02. This work was also supported by the EUROfusion Consortium,
funded by the European Union via the Euratom Research and Training Programme
(Grant Agreement No. 101052200 EUROfusion). The views and opinions expressed
are, however, those of the author(s) only and do not necessarily reflect those
of the European Union or the European Commission. Neither the European Union
nor the European Commission can be held responsible for them.

## References

  * Chen and Zonca (2016) L. Chen and F. Zonca, Review of Modern Physics 88, 015008 (2016).
  * Cheng, Chen, and Chance (1985) C. Cheng, L. Chen, and M. Chance, Ann. Phys. 161, 21 (1985).
  * Cheng and Chance (1986) C. Z. Cheng and M. S. Chance, The Physics of Fluids 29, 3695 (1986).
  * Heidbrink  _et al._ (1991) W. Heidbrink, E. Strait, E. Doyle, G. Sager, and R. Snider, Nuclear Fusion 31, 1635 (1991).
  * Heidbrink  _et al._ (1993) W. Heidbrink, E. Strait, M. Chu, and A. Turnbull, Phys. Rev. Lett. 71, 855 (1993).
  * Rosenbluth  _et al._ (1992) M. N. Rosenbluth, H. L. Berk, J. W. Van Dam, and D. M. Lindberg, Phys. Rev. Lett. 68, 596 (1992).
  * Zonca and Chen (1992) F. Zonca and L. Chen, Phys. Rev. Lett. 68, 592 (1992).
  * Zonca and Chen (1993) F. Zonca and L. Chen, Physics of Fluids B: Plasma Physics 5, 3668 (1993).
  * Chen (1994) L. Chen, Physics of Plasmas 1, 1519 (1994).
  * Chu  _et al._ (1992) M. S. Chu, J. M. Greene, L. L. Lao, A. D. Turnbull, and M. S. Chance, Physics of Fluids B: Plasma Physics 4, 3713 (1992).
  * Falessi  _et al._ (2019) M. V. Falessi, N. Carlevaro, V. Fusco, G. Vlad, and F. Zonca, Physics of Plasmas 26, 082502 (2019).
  * Falessi  _et al._ (2020) M. V. Falessi, N. Carlevaro, V. Fusco, E. Giovannozzi, P. Lauber, G. Vlad, and F. Zonca, Journal of Plasma Physics 86, 845860501 (2020).
  * Chen and Zonca (2017) L. Chen and F. Zonca, Physics of Plasmas 24, 072511 (2017).
  * Pucella  _et al._ (2022) G. Pucella, E. Alessi, F. Auriemma, P. Buratti, M. V. Falessi, E. Giovannozzi, F. Zonca, M. Baruzzo, C. D. Challis, R. Dumont, D. Frigione, L. Garzotti, J. Hobirk, A. Kappatou, D. L. Keeling, D. King, V. G. Kiptily, E. Lerche, P. J. Lomas, M. Maslov, I. Nunes, F. Rimini, P. Sirén, C. Sozzi, M. F. Stamp, Z. Stancar, H. Sun, D. V. Eester, M. Zerbini, and J. Contributors, Plasma Physics and Controlled Fusion 64, 045023 (2022).
  * Note (1) The access to the code repository can be obtained by writing to the code maintainer at<EMAIL_ADDRESS>
  * Piron  _et al._ (2023) L. Piron, P. Buratti, M. Falessi, M. Gambrioli, G. Graham, M. Lennhol, and D. Valcarcel,  (2023), 10.2139/ssrn.4511126.
  * Connor, Hastie, and Taylor (1978) J. Connor, R. Hastie, and J. Taylor, Phys. Rev. Lett. 40, 396 (1978).
  * Lu, Zonca, and Cardinali (2012) Z. Lu, F. Zonca, and A. Cardinali, Physics of Plasmas 19 (2012).
  * Albanese and Pizzuto (2017) R. Albanese and A. Pizzuto, Fusion Engineering and Design 122, 274 (2017).
  * Zeeland  _et al._ (2019) M. V. Zeeland, C. Collins, W. Heidbrink, M. Austin, X. Du, V. Duarte, A. Hyatt, G. Kramer, N. Gorelenkov, B. Grierson, D. Lin, A. Marinoni, G. McKee, C. Muscatello, C. Petty, C. Sung, K. Thome, M. Walker, and Y. Zhu, Nuclear Fusion 59, 086028 (2019).
  * Wang  _et al._ (2018) T. Wang, Z. Qiu, F. Zonca, S. Briguglio, G. Fogaccia, G. Vlad, and X. Wang, Physics of Plasmas 25, 062509 (2018).
  * Wang  _et al._ (2019) T. Wang, X. Wang, S. Briguglio, Z. Qiu, G. Vlad, and F. Zonca, Physics of Plasmas 26, 012504 (2019).
  * Zonca and Chen (2014a) F. Zonca and L. Chen, Physics of Plasmas 21, 072120 (2014a).
  * Zonca and Chen (2014b) F. Zonca and L. Chen, Physics of Plasmas 21, 072121 (2014b).
  * Floquet (1883) G. Floquet, Annales scientifiques de l’École Normale Supérieure 2e série, 12, 47 (1883).
  * Albanese  _et al._ (2003) R. Albanese, G. Calabrò, M. Mattei, and F. Villone, Fusion Engineering and Design 66-68, 715 (2003), 22nd Symposium on Fusion Technology.
  * Lütjens, Bondeson, and Sauter (1996) H. Lütjens, A. Bondeson, and O. Sauter, Computer Physics Communications 97, 219 (1996).
  * Qiu, Zonca, and Chen (2012) Z. Qiu, F. Zonca, and L. Chen, Physics of Plasmas 19, 082507 (2012).
  * Chen and Hasegawa (1974) L. Chen and A. Hasegawa, The Physics of Fluids 17, 1399 (1974).
  * Hasegawa and Chen (1974) A. Hasegawa and L. Chen, Phys. Rev. Lett. 32, 454 (1974).
  * Li  _et al._ (2023) Y. Li, M. V. Falessi, P. Lauber, Y. Li, Z. Qiu, G. Wei, and F. Zonca, Plasma Physics and Controlled Fusion 65, 084001 (2023).
  * Camenen  _et al._ (2007) Y. Camenen, A. Pochelon, R. Behn, A. Bottino, A. Bortolon, S. Coda, A. Karpushov, O. Sauter, G. Zhuang, and the TCV team, Nuclear Fusion 47, 510 (2007).
  * Austin  _et al._ (2019) M. E. Austin, A. Marinoni, M. L. Walker, M. W. Brookman, J. S. deGrassie, A. W. Hyatt, G. R. McKee, C. C. Petty, T. L. Rhodes, S. P. Smith, C. Sung, K. E. Thome, and A. D. Turnbull, Phys. Rev. Lett. 122, 115001 (2019).
  * Nelson  _et al._ (2023) A. O. Nelson, L. Schmitz, C. Paz-Soldan, K. E. Thome, T. B. Cote, N. Leuthold, F. Scotti, M. E. Austin, A. Hyatt, and T. Osborne, Phys. Rev. Lett. 131, 195101 (2023).
  * Ghai  _et al._ (2021) Y. Ghai, D. A. Spong, J. Varela, L. Garcia, M. A. V. Zeeland, and M. E. Austin, Nuclear Fusion 61, 126020 (2021).
  * Oyola et al. (2023) P. Oyola et al., "Mitigation of Alfvén Eigenmodes in Negative triangularity plasmas at the TCV tokamak", 29th IAEA Fusion Energy Conference (2023).
  * Zonca  _et al._ (2015) F. Zonca, L. Chen, S. Briguglio, G. Fogaccia, G. Vlad, and X. Wang, New Journal of Physics 17, 013052 (2015).