Physics of drift Alfvén instabilities and energetic particles in fusion plasmas

Shear Alfvén wave (SAW)/drift Alfvén wave (DAW) fluctuations can be destabilized by energetic particles (EPs) as well as thermal plasma components, which play a key role in the EP energy and momentum transport processes in burning fusion plasmas. The drift Alfvén energetic particle stability (DAEPS) code, which is an eigenvalue code using the finite element method, was developed to analyze Alfvén instabilities excited by EPs. The model equations, consisting of the quasineutrality condition and the Schrödinger-like form of the vorticity equation, are derived within the general fishbone-like dispersion relation theoretical framework, which is widely used to analyze SAW/DAW physics. The mode structure decomposition approach and asymptotic matching between the inertial/singular layer and ideal regions are adopted. Therefore, the DAEPS code can provide not only frequency and growth/damping rate but also the parallel mode structure as well as the asymptotic behavior corresponding to the singular-layer contribution. Thus, it fully describes fluid and kinetic continuous spectra as well as unstable and damped modes. The model equations have been extended to include general axisymmetric geometry and to solve for the response of circulating and trapped particles by means of the action-angle approach. In this work, we discuss linear dispersion relation and parallel mode structure of drift Alfvén instabilities excited by EPs, computed with the DAEPS code with realistic experimental plasma profile and magnetic configuration. We compare DAEPS results with FALCON/LIGKA to provide a verification of the code. We then adopt the Dyson–Schrödinger model (DSM) to further analyze the EP energy and momentum flux. We will briefly discuss how the parallel mode structure of the drift Alfvén instabilities can be used in the DSM to calculate the nonlinear radial envelope evolution and the EP transport.


Introduction
In strongly magnetized reactor-relevant fusion plasmas, the input power density is dominated by fusion alpha particles, which are characterized by a nearly isotropic distribution function in velocity space and a decreasing in energy due to Coulomb collisions predominantly with thermal electrons. Thereby, collective instabilities in the range of the cyclotron frequency will not be excited, and the fluctuation spectrum, driven by configuration space gradients as free energy source, will be characterized by the dynamic frequencies of the particle guiding center motions [1]. In particular, fusion alpha particles and, more generally, suprathermal or energetic particles (EPs) can excite shear Alfvén waves (SAWs), whose group velocity is parallel to the ambient magnetic field and thus have a very efficient wave-particle power transfer. Due to the underlying universal instability drive [2] connected with spatial gradients, fluctuation growth rates are increasing with the mode number up to a maximum, when the perpendicular wavelength becomes comparable with the characteristic particle orbit width, and orbit averaging makes the waveparticle power transfer ineffective. Equilibrium magnetic field configurations and plasma nonuniformities play very important roles in determining the plasma stability properties via particle drifts and must be properly accounted for to make reliable predictions in realistic plasma conditions [1]. When thermal plasma particle drifts importantly affect plasma stability and fluctuation dynamics, fluctuation spectra are generally referred to as drift waves (DWs) or drift Alfvén waves (DAWs), depending on their predominant acoustic or Alfvénic polarization. In this work, we focus on finite-β (β being the ratio of kinetic to magnetic energy density) fusion plasmas and the interaction of EPs with SAW and (DAW) fluctuations.
The SAW/DAW fluctuation spectrum can be destabilized by EPs as well as thermal plasma components in fusion tokamak devices via wave-particle resonances. These fluctuations are characterized by a broad spectrum of wavelengths, frequencies, and growth rates, which can significantly influence the saturation level of the instabilities and thus play an important role in the EP transport processes in burning fusion plasmas [3][4][5]. In fact, SAW/DAW spatial scales cover the range from the thermal ion Larmor radius microscale to the system macroscales, including the various intermediate mesoscales, consisting of the EP Larmor radius 7 and the characteristic radial correlation length of structures that are formed within burning plasmas as complex self-organized system [6][7][8]. Meanwhile, the corresponding frequency spectra cover 7 The EP Larmor radius plays an interesting double role of microscale for EP-driven instabilities and mesoscale for DAW and electrostatic drift wave turbulence [5]. the whole range from low magnetohydrodynamic (MHD) frequencies up to the typical SAW frequency ∼ v A /L, with v A the Alfvén speed and L the characteristic system size that, in a tokamak, is L ∼ qR 0 , with q the safety factor value on the magnetic axis and R 0 the magnetic axis major radius. A particularly important role in the SAW/DAW physics is played by the continuous spectrum [9], which generally exists in magnetized nonuniform plasmas and, within the ideal MHD description, is connected with the spontaneous formation of radial singular structures due to phase mixing [10,11]. This property naturally links the macroscales that are involved in the SAW resonant excitation by EPs with the microscales typical of DAW. Thus, the necessity of a gyrokinetic analysis becomes evident. Further, for a proper assessment of the nonlinear behavior of SAW/DAW and the corresponding EP transport, it is important to accurately describe both the stable and unstable components of the fluctuation spectrum [10,12,13].
The Alfvén eigenmode (AE)-induced EP transport has been widely observed in tokamak experiments [14][15][16][17][18][19], which may lead to the degradation of plasma confinement and damage to the first wall in fusion devices. Therefore, it is crucial to fully understand the AE physics and, on this basis, identify the various modes with different characteristics and estimate their contributions to EP transport. The resonant character of EP transport induced by SAW/DAW fluctuations calls for the analysis of EP fluxes in the phase space [1,4,5]. These transport processes can be fully accounted for by the nonlinear evolution of phase-space zonal structures (PSZSs) [4,20], which are the component of the toroidally symmetric (n = 0) particle response that is constant along equilibrium orbits and thus is undamped by collisionless processes. These PSZSs are the phase-space counterpart of the zonal field structures (ZFSs) [21], viz., the well-known zonal flows and currents [22][23][24][25], which play crucial roles in the regulation of fluctuation level and fluctuation-induced transport because they scatter primary instabilities into a shorter-wavelength stable domain [26]. Together with the finite ambient fluctuation level, PSZS and ZFS constitute the zonal state, whose evolution can be understood as a succession of neighboring nonlinear equilibria [27], which is described by a general transport theory written in conservative form [20,21]. The drift Alfvén energetic particle stability (DAEPS) code, which is a non-perturbative eigenvalue code using the finite element method (FEM), was developed to analyze Alfvén instabilities excited by EPs [28]. The model equations were derived within the general fishbone-like dispersion relation (GFLDR) theoretical framework [29,30], which is widely used to analyze SAW/DAW physics. The mode structure decomposition (MSD) method [31] and asymptotic matching between the inertial/singular layer and ideal regions were adopted. Therefore, the DAEPS code can provide frequency and growth/damping rates as well as the asymptotic behavior of the parallel mode structure corresponding to the singular-layer contribution. Thus, it fully describes fluid and kinetic continuous spectra as well as unstable and damped modes. Recently, the DAEPS model equations have been extended to include general axisymmetric geometry and to solve for the response of circulating and trapped particles by means of the action-angle approach [32]. In this paper, recent updates to the DAEPS code will be reported, including theoretical model and application studies to the interaction between EPs and toroidal Alfvén eigenmode (TAE) instability [33] in realistic geometry. In particular, section 2 will introduce the DAEPS governing equations. Section 3 introduces the realistic plasma equilibria adopted as reference cases in DAEPS applications in the present work. These are the divertor tokamak test (DTT) facility reference scenario [34] and the 15 MA ITER scenario [35]. Meanwhile, section 4.1 presents an in-depth verification of the DAEPS gyrokinetic calculation of continuous spectra in DTT and ITER. Verifications against well-tested results of the Floquet Alfvén Continuum (FALCON) code [36] in the ideal MHD fluid limit are given first for both DTT and ITER. Then, a discussion is provided on the verification of the gyrokinetic results for the ITER continuous spectrum against the Linear Gyrokinetic Shear Alfvén Physics (LIGKA) code [37]. After a complete verification of the capability of DAEPS to address the relevant gyrokinetic description on microscales, section 4.2 is devoted to a linear stability analysis of the high toroidal mode number (n = 20) TAE instability in the DTT reference case. Based on these results, the TAE-induced EP phase-space fluxes are computed in section 5 after a brief summary of the foundations of the PSZS transport theory, including the three-level hierarchy of different approximations that allow constructing reduced EP transport models within the same theoretical framework, i.e. the Dyson-Schrödinger model (DSM) [4,20,21]. As illustrative examples, fluxes are explicitly computed in the quasilinear limit, showing the capability of DAEPS to evaluate the evolution of the equilibrium distribution function as well as the corresponding EP transport. Finally, a summary and discussion are given in section 6, with a description on how the DAEPS results on parallel mode structures and gyrokinetic particle responses will be adopted in the future to self-consistently calculate the nonlinear radial envelope evolution and the EP transport in the DSM. Various technical details on the definition of the continuous spectrum (appendix A), the numerical implementation of action-angle variables (appendix B), and phasespace particle fluxes and transport equations (appendix C) are given in three additional appendices.

Governing equations
The model equations for the DAEPS code [28] originate from the GFLDR framework [29,30], which has been intensively used to analyze SAW and DAW physics by adopting the MSD method [31] and asymptotic matching between the inertial region/singular layer and the ideal region. The GFLDR, for a single n toroidal mode, takes the form of i|ŝ|Λ n = δŴ nf + δŴ nk , where the generalized inertial term Λ n is the normalized singular-layer contribution, including kinetic response;ŝ is the magnetic shear; and δŴ nf and δŴ nk are the fluid and kinetic contributions of the potential energy, respectively. The global dispersion relation can be formulated in a quadratic form due to the variational nature of the GFLDR. A trial function is needed to calculate the frequency, growth rate, fluid and kinetic potential energies, and their asymptotic behaviors, which need to be accurate enough for a precise calculation of the linear eigenvalues. For high toroidal mode number n, when the radial envelope width is broader than the characteristic perpendicular wavelength, the GFLDR can be cast as a local dispersion relation iΛ n = δW nf + δW nk , where δW n = δŴ n /|ŝ| for localized modes with radial width smaller than the characteristic scale length of equilibrium profiles.
To deal with the general tokamak geometry for arbitrary toroidal mode number n, we make use of the MSD method in Boozer coordinates (ψ p , θ, ζ) [38], where ψ p is the poloidal magnetic flux function, θ is the poloidal angle, ζ is the generalized toroidal angle, and the covariant and contravariant forms of magnetic field are given by Here, q(ψ p ) is the safety factor; g(ψ p ), I(ψ p ), and δ(ψ p , θ) are the contravariant components of the magnetic field. For a given toroidal mode number n, any fluctuation f(r, θ, ζ) can be written as [31] f(r, θ, ζ) = where r(ψ p ) is a radial-like flux coordinate; ϑ is the extended poloidal angle in ballooning space, corresponding to the normalized parallel length along the magnetic field line; and any residual time dependence is left implicit. Furthermore, in the second line of equation (2), we made use of the Poisson summation formula, i.e.
∑ ℓ∈Z δ(ϑ − θ + 2πℓ), and δ(. . .) is the Kronecker delta function. For high-n mode with finite magnetic shear, the most unstable modes are localized around the rational surface, and the slow variation of the envelope can be separated from the fast radial variation of the mode structure related to the ∼(nq − m) dependence on the parallel wave vector. In this case, the MSD method can be reduced to the ballooning representation, Note that, here,f n (ϑ) generally contains residual radial dependence on the plasma equilibrium scale, while the remaining implicit time dependence can all be included in the radial envelope A n (r) → A n (r, t).
Using the DAEPS code, considering the low-β plasma of interest here, we describe the plasma oscillations in terms of the electrostatic potential δϕ and the magnetic scalar potential δψ, which is related to the parallel to B component of the vector potential δA ∥ by δA ∥ = −icω −1 b · ∇δψ. The parallel electric field can then be expressed by δE ∥ = −b · ∇(δϕ − δψ). Meanwhile, the compressional component of the magnetic field fluctuation, δB ∥ , is being solved explicitly by means of the perpendicular (to B) pressure balance BδB ∥ + 4π δP ⊥ = 0, with δP ⊥ the overall perpendicular pressure response. The model equations are derived within the GFLDR theoretical framework, consisting of the ballooning space vorticity equation in general geometry, and the quasineutrality condition where we have dropped the subscript n and the 'hat' on field variables in the ballooning space to simplify notations, Ψ = κ ⊥ δψ corresponds to the magnetic scalar potential with is related to the parallel electric field, and ∑ s denotes summation on particle species s. The gyrocenter distribution function δKs of species s, accounting for kinetic particle compression response, is solved through the linear gyrokinetic equation, In equation (4), the lhs accounts for the fluid contributions, consisting of field line bending, inertia, and ballooning interchange terms, and the rhs represents the kinetic particle compression response. Furthermore, is the Alfvén frequency, with V A0 = B 0 / √ 4πρm the Alfvén speed on the magnetic axis, B 0 is the on-axis magnetic field, ρm is the equilibrium plasma mass density, and r = B 0 /B(r, ϑ). Meanwhile, ω T * ps = (k θ cTs)/(esBLps) and ω T ds = (k θ cTs)/(esBR 0 ) correspond to the diamagnetic frequency and magnetic drift frequency, respectively, withB(r) = qψ ′ p /r, prime denoting derivation with respect to r, Lps = −∂r ln Ps, and R 0 is the magnetic axis major radius. We have also introduced the notation g = κn +ŝϑκg, with κn = −R 0 ∂r ln B and κg = gq gq+I R0 r ∂ ϑ ln B denoting the normal and geodesic curvature, respectively. Finally, βs = 8π Ps/B 2 0 is the ratio of plasma pressure to magnetic pressure, and x = v/vts is the particle velocity normalized by the thermal speed, vts = (2Ts/ms) 1/2 . In equation (6), J b = (∇ψp × ∇ϑ · ∇ζ) −1 = (gq + I)/B 2 is the Jacobian of Boozer coordinates, and ∂ ∥ = (J b B) −1 ∂ ϑ represents the parallel derivative along the magnetic field. Furthermore, ω ds = ω T ds msE Ts g is the magnetic drift frequency, with E = v 2 /2 = x 2 Ts/ms the particle kinetic energy per unit mass and λ = rv 2 ⊥ /v 2 the pitch angle. Finally, QFs = (ω∂ E +ω * s) Fs accounts for the free energy provided by the phase-space gradient of the equilibrium distribution function Fs, andω * sFs = Ω −1 cs k × b · ∇Fs with Ωcs the cyclotron frequency. We note that equation (4) is an extension to general geometry of the vorticity equation used in our previous model [28]. It can be reduced to equation (3) of [28] by adoptingŝ-α equilibrium model and the ballooning representation [39].

Reference scenarios
In this section, we briefly introduce the reference scenarios that will be used for DAEPS verification and for its application to the calculation of EP fluxes in phase space. Because DAEPS [28] was originally developed to deal with (ŝ, α) model equilibria [39], its governing equations and code structure were modified to account for a realistic magnetic equilibrium, as discussed in section 2. The DTT reference scenario is briefly introduced in the following, while we only refer to the ASTRA-simulated ITER 15 MA scenario (131018; run50 in the public ITER IMAS scenario database) [35]. We will use these two scenarios in the upcoming sections devoted to DAEPS applications.
The DTT facility is a D-shaped superconducting device, which is being built at the ENEA Research Center in Frascati, Italy, which mainly aims at investigating viable and/or advanced power exhaust solutions for the demonstration power plant (DEMO), thereby filling gaps between present-day devices, ITER [40], and DEMO [41]. Although the DTT is designed as a relatively compact machine focused on divertor and power exhaust issues, it is capable of reproducing edge as well as core plasma conditions that are relevant for burning plasma studies.
We consider a single null setting, whose basic profiles are depicted in figures 1-3 as a function of the normalized toroidal radius ρtor = r/a, where a is the tokamak minor radius, and ρtor is the square root of the normalized toroidal flux. We note that, for the analyzed case, the kinetic pressure on axis is 948 kPa, the flux function F on axis is 13.7 V s, while the density is 2.3 × 10 20 m −3 , the electron and ion temperatures are, respectively, 15.3 and 9.4 keV. All the equilibrium quantities have been calculated by means of self-consistent transport studies, i.e. see [34]. For simplicity, we consider a pure deuterium plasma, i.e. ne = n i = n. Consistent with the above, for stability studies, the EP population will be assumed as deuterium with a Maxwellian distribution function (for details, see section 4.2).

DAEPS applications
Here, we illustrate applications of the DAEPS code to the analysis of DAW in realistic plasmas of fusion interest. In particular, we consider the reference scenario [34] for the DTT facility [42,43], which has been described in detail in section 3. This DTT reference scenario is adopted in section 4.1 to describe the properties of the SAW and ion sound wave (ISW) continuous spectra when fully gyrokinetic particle compression responses are accounted for [3,30,[44][45][46][47], as  well as in the fluid limit [48][49][50][51]. Section 4.1 is based on comparisons of the DAEPS results with those by FALCON [36] and LIGKA [37] codes. Having demonstrated the ability of DAEPS to properly treat various DAW dynamics at small radial scales, where they are typically most challenging to describe for a numerical approach, section 4.2 is devoted to the highn stability analysis of TAEs in the DTT reference scenario presented in section 3. The results confirm the earlier findings by Wang et al [52], showing that plasma core region can be thought as being divided into two subregions: the inner core, where DAW excitation by EPs is strongly influenced by the weak magnetic shear and is dominated by frequencies lower than the TAE gap, and the outer core, with typically finite magnetic shear and characterized by fluctuations in the TAE frequency range. The results of section 4.2 also confirm that the EP-driven DAW spectrum in DTT is expected to have toroidal mode numbers n ≳ O(10) [52], similar to ITER and typical of burning plasmas [1].

Kinetic continuous spectra
In the following section, we will compare the continuous spectra computed with DAEPS with those computed with the FALCON code [51], which are obtained by solving equations (A.4) and (A.5). Details about the FALCON code can  be found in [51,53] and will be omitted here. Similarly, we will compare the kinetic continuous spectra computed with DAEPS with those obtained by the LIGKA code [37]. In this way, we will provide a verification of the DAEPS model and the governing equations concerning the capability of DAEPS to handle the various physics processes in the kinetic layer region.
The comparison of the MHD low-frequency electromagnetic continuous spectrum computed with DAEPS and FALCON for n = 5 using the DTT reference equilibrium discussed in section 3 is indistinguishable. For this reason, figure 4 illustrates the comparison of the fluid (FALCON) and kinetic (DAEPS) continuous spectra for the same DTT equilibrium. Because the diamagnetic frequency and β are small, the overall effect of ion sound branches is correspondingly small in the kinetic calculation. In particular, figure 4 shows the kinetic continua, where TAEs can be located, or beta-induced Alfvén eigenmodes (BAEs), as well as beta-induced Alfvén acoustic eigenmodes (BAAEs), kinetic ballooning modes (KBMs), and ISWs. The comparison of the fluid and kinetic continuous spectra suggests that the kinetic structures of low-frequency spectra are consistent with the fluid results, especially for the strong Alfvénic polarization branches, which correspond to the slow sound approximation due to decoupling of SAWs and ISWs [53]. This is also consistent with the electron temperature being larger than the ion temperature across the whole plasma in DTT, and with the fact that the kinetic continuous spectrum in the fluid limit is expected to recover the MHD result for Γ = (7/4 + Te/T i )/(1 + Te/T i ) [54]. The frequency and radial location of the n = 20 TAE in the DTT reference equilibrium are also shown, which are the modes that will be taken as reference in the discussions of sections 4.2 and 5.
To further investigate these issues connected with the fluid and kinetic structures in the continuous spectra, we have compared the DAEPS and FALCON results with those of the LIGKA code [37] using the 15 MA ITER reference case as equilibrium [32]. Similar to the DTT case discussed earlier, the MHD continua computed with DAEPS and FALCON are indistinguishable. Thus, the comparison of the fluid low-frequency electromagnetic continuous spectrum with n = 5 of DAEPS and FALCON is shown in figure 5 (left panel), where FALCON uses the original equilibrium metric tensor, while DAEPS uses the analytiĉ s − α model, reconstructed by the EQUIPE post-processing code [36], i.e. a circular geometry approximation of the actual ITER equilibrium, which suggests that the low-frequency continuous spectrum is sensitive to realistic magnetic geometry. The continuous spectrum using the actual numerical equilibrium for the ITER reference case is shown in figure 5 (right panel) and suggests that the frequency of the kinetic structure of the low-frequency spectra is qualitatively consistent with the fluid result due to the small diamagnetic frequency in the n = 5 case. Figure 6 (left panel) shows the comparison of the frequency of the kinetic continuous spectrum calculated by DAEPS (colored lines) and LIGKA (black markers) for the ITER 15 MA equilibrium. It is suggested that the frequencies of the TAE, BAE, and KBM branches are highly consistent, while the frequency of the BAAE branch is qualitatively consistent near the accumulation point. The growth/damping rates, shown in figure 6 (right panel), are qualitatively consistent and illuminate the existence of multiple strongly damped modes identified by LIGKA. In fact, it is known that the kinetic expression of the low-frequency continuous spectrum is represented by a transcendental function with infinitely many and heavily damped roots [55]. These branches have practically no relevance in applications due to their strongly damped nature. Nonetheless, their features can be a useful benchmark test in numerical calculations, and thus, they are being further analyzed in the benchmarking efforts of the DAEPS and LIGKA codes. Of more important implications can be the marginal unstable features of the low-frequency (yellow) modes identified by DAEPS near mode rational surfaces. These are clearly connected with the low-frequency modes discussed recently concerning experimental observations in DIII-D [47,[56][57][58] and will be analyzed in a separate work.

TAE excitation by EPs
After the analysis of the SAW/ISW continuous spectra in DTT, presented in section 4.1, we consider the excitation of high-n AEs in DTT by a finite population of EPs. DTT stability analyses using the DAEPS code confirm earlier findings [52] that properties of the Alfvénic fluctuation spectrum in DTT are characterized by an inner core region, where magnetic shear is small and dominant modes are typically located in the lowfrequency region below the TAE gap, and an outer core region with finite magnetic shear and typical AEs in the TAE frequency range.
For the sake of conciseness, we focus on high-n TAE modes in the outer core region, and as a paradigmatic example, we analyze an n = 20 TAE excited by transit resonance with EPs. Because the details of additional heating in DTT are still being discussed [59], the focus here is not on the actual use of a realistic EP distribution function, but rather on the demonstration that n ≳ O (10) are effectively excited in DTT [52]. As a reasonable guess, we assume that EPs are Maxwellian, with a local density of 1% w.r.t. that of the thermal plasma and the same temperature profile as the thermal ions, resulting in a constant EP-to-thermalion-energy-density ratio β EP /β i = 0.435 (cf section 3 for details).
We first analyze the TAE mode structure and stability in the absence of EP drive in the DTT reference scenario, where the plasma is assumed to be of pure deuterium. The n = 20 TAE mode localized near the q = 1.1 surface has a two-dimensional mode structure that is shown in figure 7 (left panel) and mode frequency ω/ω Ti = 6.5714 − 0.0080 668i. Thus, as expected, it is marginally stable due to weak Landau damping. The corresponding parallel mode structure in the ballooning space is shown in figure 8, suggesting that the parallel electric field is negligibly small compared to the magnetic scalar potential, which indicates the predominant Alfvénic polarization of the instability. This case corresponds to the rhs of equation (10) being subdominant. If needed, finite compression effects on the continuous spectrum can be included by solving equation (11) algebraically (by dropping the ∼ ∂ 2 ϑ term), that is, adopting the slow sound approximation [51]. The radial envelope in figure 7 , is obtained for the ground state of the linear TAE spectrum, where r 0 is the radius of the q = 1.1 surface, and [60,61] D(ω, r, θ k ) denotes the local TAE dispersion function, and θ k ≡ kr/(nq ′ ), with kr the radial envelope wave number. Note that ∆r ∼ r 1/2 0 |nq ′ | −1/2 , and thus, the frequency shift of the global mode with respect to the solution of the local dispersion relation at r = r 0 is ∆ω/ω ∼ |nr 0 q ′ | −1 or, more precisely [60,61] This well-known result [39] shows that the local solution of the dispersion relation well reflects the global mode frequency for this kind of radially localized modes, although this simple result does not hold in the most general case [60,61]. Next, we repeat the stability analysis for the same n = 20 TAE mode localized near the q = 1.1 surface but in the presence of 1% of deuterium particles that represent the Maxwellian EP population with β EP /β i = 0.435. The twodimensional mode structure is qualitatively similar to the case with no EPs, shown in figure 7 (left panel), with a mode frequency of ω/ω Ti = 6.86 + 0.0779i, confirming the result in [52] that n ≳ O (10) are effectively excited in DTT. The parallel mode structure in the ballooning space is shown in figure 9,  showing the predominant Alfvénic polarization of the instability also in this case. The frequency shift with respect to the previous case with no excitation by EPs is of the order of the mode frequency separation from the upper accumulation point, as shown in figure 4, and the comparison of parallel mode structures in figures 8 and 9 suggests that n = 20 TAE excitation by EP is in the transition region from perturbative to non-perturbative EP behavior. We refer interested readers to the discussions in [52,62] for more details on this interesting point, which is one of the characteristic features of DTT that makes it relevant for burning plasma studies.

Phase-space fluxes
In the previous sections, we have introduced the DAEPS code, and we have shown some applications to the continuous spectrum and AE stability studies. In this section, the numerical results computed with DAEPS, consisting of the parallel linear mode structure, frequency, and growth rate, are exploited to calculate the EP phase-space fluxes induced by the presence of electromagnetic fluctuations. A detailed derivation of the phase-space flux expressions is given in appendix C for the reader's convenience.    The PSZS dynamics (see [20,21]) is described by a continuity equation in the phase space where e iQz (. . .) describes the average of a scalar function in the phase space over a particle orbit in the equilibrium magnetic field, that is, a bounce/transit average with period τ b,t combined with a shift operator, i.e. e iQz , accounting for particle drifts with Qz = F(ψp) [4,20]. Using the MSD approach for the fluctuating fields [31], i.e. equation (3) in section 2, and the convolution expression, as discussed in [1,30], we obtain, after integration along the equilibrium orbits, the expressions for the fluxes appearing in the PSZS equation and, in particular, the radial particle flux appearing in the second term of equation (9) that can be immediately rewritten changing integration variable from η to ϑ (11) Note that, for the sake of completeness and clarity, we have provided expressions in terms of quadratures over the canonical angle η as well as in ballooning angle ϑ. Corresponding expressions can be written for the phase-space energy flux. By direct inspection, we note that the wave intensity dependence is isolated, and it bears the mesoscale information of the phase-space fluxes. We also note that the scaling of the fluxes with |A(r, t)| 2 is universal and has been noted in the numerical simulation of DW-induced turbulent fluxes [63] as well as in the corresponding theoretical studies [26,64,65]. Meanwhile, the parallel mode structure is left implicit and is integrated over, and it reflects the equilibrium macroscale structures of phase-space fluxes. Furthermore, the fine-scale corrugations produced by fluctuation-induced fluxes are accounted for by the convolution sum and by the ∼ exp (2π inq(r)ℓ) phase factors. Therefore, the present expression for the phase-space fluxes illuminates the characteristic spatiotemporal scales of turbulent fluxes produced by the fluctuating fields. We can obtain an even more detailed analysis and insight into the underlying physics processes that are responsible for the phasespace transport by separating the resonant from the nonresonant particle responses, as shown in appendix B. The flux expression appearing in equation (11) can be immediately rewritten using the own notation of DAEPS. After normalizing to the bounce time and restoring the subscript n and the 'hat' on field variables in the ballooning space, which were dropped after equation (5) to simplify notations, we obtain where, for simplicity, we have assumed a single particle species s, we have specialized the phase-space flux expression to a single toroidal mode number n and noted δLn = (esδLn/Ts) = cĀ ∥ , δr = δr/R 0 , δĒ = 2δE/v 2 ts . As a particular application, we now evaluate the radial flux of fast particles induced by the resonantly excited TAE in the DTT reference case and described in the previous section. For the sake of simplicity, we focus our analysis on the numerical calculation of the quasilinear limit of the phase-space particle flux [66]. In fact, a more detailed analysis of this case requires separate work and will be presented elsewhere. We show that the kinetic compression term dominates the resonant particle contribution, as shown in figure 10. The nonresonant particle contribution, meanwhile, appears to be almost negligible. These two terms are plotted for the radial position highlighted in blue in figure 7, while we refer interested readers to appendix B for the separation and classification of the different contributions to the general expression of the phase-space particle flux (equation (11)). The dominant contribution from the resonant particle response in the kinetic compression term is qualitatively consistent with theoretical predictions [1,4] and numerical simulation results [62]. Meanwhile, the peak of the particle flux at E ≃ (5/2)T/m can be understood as a maximization of the wave particle power transfer for a Maxwellian distribution and is consistent with the discussions of [67,68] on the construction of an equivalent EP Maxwellian mimicking a slowing-down EP distribution function. These results demonstrate the robustness of the DAEPS analysis of parallel mode structures and give confidence in its applicability to the truly nonlinear calculation of phase-space fluxes within the DSM approach [4,21].

Summary and discussion
In this work, the local linear gyrokinetic description of finite-β nonuniform magnetized plasmas [69,70] has been implemented in the DAEPS code [28] in general geometry using straight magnetic field line toroidal flux coordinates [71,72]. More precisely, we have implemented the low-β limit of these equations, where the dynamics can be described in terms of two independent scalar fields, viz. the scalar and parallel vector potentials, because the parallel magnetic field perturbation can be explicitly solved for by means of the perpendicular pressure balance. We have also adopted a formulation based on the quasineutrality condition and gyrokinetic vorticity equation [22,70,73], which are equivalent to Poisson's and parallel Ampère's laws but are more convenient and transparent in their physical content, as they readily recover the MHD and extended kinetic MHD descriptions in the proper limit [1].
For our intended applications, it is crucial that DAEPS is an eigenvalue code, which allows us to investigate unstable as well as stable parts of the low-frequency electromagnetic fluctuation spectra that are both essential players in the nonlinear plasma evolution. Meanwhile, solving the aforementioned local gyrokinetic equations in ballooning space and general geometry gives an accurate description of the plasma kinetic response at micro-and mesoscales, which is mandatory. In fact, short-wavelength DAWs are naturally excited in nonuniform magnetized fusion plasmas due to the ubiquitous nature of continuous spectra. At the same time, EP excitations of n ≳ O (10) Alfvénic oscillations generally involve linear as well as nonlinear structure formations on mesoscales [1,4,5].
Here, we have verified plasma kinetic descriptions on the microscales in DAEPS by detailed comparisons against kinetic continuous spectra computed with the LIGKA code [37] for the 15 MA ITER scenario. Good agreement is found for the least damped continuous spectra, while some remaining discrepancy for strongly damped modes may be due to the transcendental nature of the local radially singular plasma kinetic response. These issues will be investigated further in the future, although their impact on practical applications is negligible, as fluctuations of this kind are suppressed. Further verifications of DAEPS against continuous spectra computed with the FALCON code [51] for both 15 MA ITER and DTT reference scenarios are also discussed with emphasis on realistic geometry and plasma compressional responses.
Alfvénic mode stability is investigated with DAEPS for the DTT reference scenario, showing that n ≳ O (10) Alfvénic oscillations are effectively excited by EPs, consistent with earlier analyses [52,62]. We focus on the n = 20 TAE mode excited by circulating EPs in the DTT outer core region, characterized by finite magnetic shear. A significant finding for this instability is that, for the present parameters, the mode is in the transition region between perturbative and non-perturbative EP effects, as described in [52,62] and the general theoretical framework [1], where non-perturbaviness is referred to the fact that the EP effects, measured as complex frequency shift, are comparable to the distance of the mode frequency from the continuous spectrum accumulation point. This is consistent with the core of DTT plasmas being in a significant regime for burning plasma application experiments, where core-edge coupling allows addressing particle and power exhaust issues in reactor-relevant conditions [43,74].
Finally, the n = 20 TAE parallel mode structures have been used to compute the EP fluxes in the phase space in the quasilinear limit [66]. This is the second-level approximation in the DSM theoretical framework for computing phase-space transport. Consistent with the general theoretical framework [4,21], at the zeroth level where mode structures and particle responses are obtained from any nonlinear gyrokinetic code, the EP fluxes are fully consistent with the corresponding turbulent transport [75,76]; at the first level, only parallel mode structures are approximated by their linear limit, but radial envelope and evolving reference equilibrium are fully nonlinear [4], while the second level of approximation recovers the quasilinear description [77,78]. More in-depth applications to illustrate in detail this multilevel approximation will be reported in future works. The advantage of the present approach consists in making explicit the microscale structures of EP phase-space fluxes, as well as their mesoscale behaviors due to radial envelope/intensity dependence. Thus, an accurate calculation of phase-space transport requires a coarser radial grid than if these explicit dependences were left implicit, allowing for a more accurate resolution of wave-particle resonances in action space, which is of crucial importance for the reliable description of nonlinear behaviors [1]. The high-n spectrum of unstable AEs and DAWs in DTT recalls the suggestive description of the Alfvénic fluctuation spectrum in reactor-relevant plasmas as a 'dense population of eigenmodes (lighthouses) with unique (equilibrium-dependent) frequencies and locations', causing 'significant multiple-TAE nonlinear interactions, yielding a diffusive redistribution' of EPs [3]. This process is expected to be well described by a quasilinear diffusive relaxation [79,80], and this is the main reason why we have chosen to focus our initial investigation of TAE-induced EP phase-space fluxes in DTT by means of the DSM in the quasilinear limit.
Nonetheless, it is also known that EP transport may occur in avalanches [67,81]. For self-consistent description of EP avalanches, the PSZS evolution occurs on the same timescale as the nonlinear radial envelope evolution, which, thus, must be computed self-consistently. However, the parallel mode structure remains essentially unaltered in this process [4,5], and thus, the first level of approximation in the DSM for the description of phase-space transport can be used to fully account for this dynamics. A detailed analysis of this important point requires a dedicated study that will be reported in a separate work.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.

Acknowledgments
The authors are indebted to Prof. Liu Chen for enlightening discussions. This work was supported by the EUROfusion Consortium, funded by the European Union via the Euratom

Appendix A. Continuous spectrum
The asymptotic behavior of magnetic scalar potential as ϑ → ±∞, which is critical for calculating the eigenfrequency and growth rate and identifying whether a mode is physical or spurious, is determined by the outgoing wave boundary condition with causality constraints [29,30,50]. By taking the limit ϑ → ∞, one can show that equation (4) reduces to a second-order differential equation with periodic coefficients. Thus, consistent with the Floquet theory, the asymptotic behavior of Ψ can be formally expressed as where P(ϑ) is a 2π periodic function corresponding to the variation of equilibrium in the poloidal direction, and ν is the Floquet 'characteristic exponent', which is connected with asymptotic behavior of the potential as well as with the inertial layer contribution in the GFLDR, i.e. Λ. In particular, ν = Λ for |ω| ≪ |ω A0 |, corresponding to the kinetic low-frequency electromagnetic continuous spectrum. Considering only the circulating particles for simplicity, our model equations for the calculation of kinetic continuous spectrum consist of gyrokinetic equation whereσ is the sign of v ∥ , and ω Ts = vtsB 0 /(J b B 2 ); the vorticity equation in the inertial layer ( where ϕs =κ ⊥ Φs, δP = −i(J b B 2 )δPcomp/(qk θ cP) denotes the normalized compressional perturbed pressure response, with δPcomp = −ΓP∇ · ξ, Γ being the ratio of specific heats, and ω 2 corresponds to the frequency of the sound wave 8 . Thus, the clear correspondence of equations (A.3) and (A.4) suggests the one on one connection of kinetic compression and MHD compression terms in the gyrokinetic and MHD descriptions, respectively. We will further comment on this important point in the following, when comparing MHD and kinetic continuous spectra.
As mentioned earlier, the DAEPS code uses the FEM to solve model equations in the ballooning space. For eigenmode analysis, the cubic B-spline FEM is used to discretize parallel electric potential and magnetic scalar potential, while for continuous spectrum studies, the spectral FEM is used to express the different harmonic components of the wave. The asymptotic behavior of the magnetic scalar potential at ϑ± can be properly handled via Neumann boundary condition By discretizing the mode structure by means of a linear combination of the finite elements, the vorticity equation can then be cast into a matrix form, which is solved by an iterative algorithm in the DAEPS code. The frequency of the mode is calculated from the eigenvalue of the matrix, while the value of Λ is fitted from the asymptotic behavior of the eigenvector of the matrix. By means of this approach, the DAEPS code can solve the fully kinetic continuous spectrum, by solving equations (8) and (9) closed by the quasineutrality condition, as well as the fluid continuous spectrum, by solving the corresponding ideal MHD equations (10) and (11). The same holds for the numerical solution of the boundary-value problem yielding eigenmode structure and frequency, as shown in section 4.2.

Appendix B. Action-angle approach
In the DAEPS code, action-angle variables are often used to simplify computations of integrals in velocity space. For magnetically trapped and passing particles, which are identified by whether the parallel velocity v ∥ =σ √ 2E(1 − λ/r) vanishes or not along the magnetic field line, we can define the bounce/transit frequency, ω b,t = 2π as the invariant of motion (action) J =¸v ∥ dl and the corresponding angle η = ω b,t´θ 0 dθ ′ /θ ′ [4]. Here,¸denotes integration along the closed particle orbit. The DAEPS code is interfaced with the EQUIPE routine of the FALCON code [51,53] that allows constructing η(ϑ; P ϕ , λ, E) and J(P ϕ , λ, E) for given constants of motion, with P ϕ the toroidal canonical angular momentum. In the small orbit limit, the action-anglevariable expressions simplify, and, for example, the action angle and transit frequency of circulating particles can be written as To deal with the singularity in the bounce/transit time arising at the critical boundary between trapped and passing particles, a new pitch-angle variable is introduced as follows, where r min denotes the minimal value of r(θ), where the magnetic field is maximum and the parallel velocity is minimum. Note that ∆ min = log (1 − r min ) − log (r min ) corresponds to well-circulating particles, whereas ∆ → ∞ to the barely passing particles. In the numerical model, the functions N (∆), η (ϑ, ∆) are calculated with cubic spline interpolation. Figure B1 shows flux surfaces with triangularity and rectangularity, respectively. This figure suggests that the function N, reflecting the particle bounce/transit time, behaves linearly with respect to the pitch-angle variable introduced earlier in the range corresponding to circulating particles. Using the new pitch-angle variable, the pitch-angle integration is transformed aŝ which converges as ∆ → ∞ and can be truncated at a proper ∆.

Appendix C. Phase-space transport equation
In [4,20], it has been shown that plasma transport can be described in terms of the collisionless undamped component of the distribution function, i.e. the PSZSs. The PSZSs are linearly stable and depend only on the particle (slowly evolving, nonlinear) equilibrium invariants of motion. Thus, they naturally extend the concept of zonal structures to the phase space and, together with the ZFSs [20,21], they fully characterize the nonlinear plasma equilibrium in the presence of a finite level of fluctuations, that is, the zonal state [20,21]. The PSZS theory allows to describe gyrokinetic transport over long timescales and EP dynamics in burning plasmas. The theoretical framework of the PSZS transport has been verified by hybrid kinetic MHD codes, such as HMGC [82] and HYMAGYC [83]. More recently, different gyrokinetic codes, e.g. [75,76], have also calculated the PSZS dynamics. Being able to compare these results with phase-space fluxes computed by means of the DAEPS code is a fundamental step to verify the theoretical framework of the PSZS transport [1,4,20] and to construct a hierarchy of reduced transport models, connecting various levels of approximations within the same unified theoretical framework [20,21]. Briefly recalling what these different approximations are, all consistently described by the PSZS transport theory, we have the following three-level hierarchy. The highest-fidelity approach, consisting in the assumptions involved in the nonlinear gyrokinetic description [84][85][86], is based on the construction of PSZS fluxes from the results provided by any global electromagnetic gyrokinetic code. A first level of further approximation, meanwhile, consists in assuming that nonlinear fluctuations are still well represented by linear parallel mode structures, while radial envelope structures are nonlinear and determined self-consistently with particle responses and the evolution of the corresponding PSZS [1,5]. It should be noted that the parallel mode structure, denoted byfn in equation (3), is formed on a timescale of approximately 1/ω, making it nearly insensitive to the nonlinear interaction that occurs on a much longer timescale of τ NL ∼ 1/γ L ≫ 1/ω. This allows the parallel mode structure to be computed from linear theory, while the nonlinear radial envelope must be treated fully nonlinearly. This reduced description, called the DSM [4,21], accounts for renormalized particle response, ZFSs, resonance broadening, and avalanches in magnetized plasmas, as it was shown recently [77,78]. The final level of approximation consists in assuming a linearized particle response, yielding expressions for phasespace fluxes that are fully consistent with the well-known quasilinear description.
In the following, we will derive the phase-space fluxes introduced with equation (9). It is worth noting that the derivation can be divided into two parts. In the first one, the gyrokinetic fluxes are written without introducing any simplifying assumptions, corresponding to the highest-fidelity description of fluctuation-induced transport introduced in section 5 as well as to the first level of approximation. In the second part of the derivation, we substitute the solution of the linear gyrokinetic equation in the PSZS flux expressions, thus obtaining the usual quasilinear limiting expressions of phase-space fluxes; see e.g. [66]. As we will see, the semi-analytic nature of this approach allows to discriminate between different contributions of the fluxes based on the underlying dynamics, e.g. related to resonant particles, kinetic compressibility, advective terms, etc.
The governing equation for the equilibrium orbit averaged distribution function is equation (13) given in the main text. The equilibrium orbit averaging can be written in the following way: where dl represents the arc length element along the particle orbit, and v ∥ is the parallel component of the particle velocity. The circle symbol in the first integration represents the periodic nature of particle orbits in the equilibrium magnetic field, with the canonical poloidal angle η ranging from 0 to 2π. Using the MSD approach [31], reducing to the well-known ballooning representation in the high toroidal mode number limit, the perturbed radial velocity and particle energy variations, after some algebra, can be written as Substituting these expressions into equations (10) and (11), we immediately obtain where, δr = δr/R 0 and δĒ = 2δE/v 2 ts . Furthermore, ignoring tearing parity and adopting the magnetic scalar potential δψ instead of the parallel vector potential, Assuming an equilibrium Maxwellian distribution, Fs, for simplicity but without loss of generality, we can introduce the following decomposition for the distribution function Fsψ + δKsn, (C.5) where the first term represents the adiabatic particle response, the second one is the advective part, while δKsn is the kinetic compression response that can be obtained by solving the following linear gyrokinetic equation ) (C.9) ×ˆdϑ r 3/2 √ r − λ e inq2π ℓ J 0s (ϑ)φ * (ϑ)ψ (ϑ − 2π ℓ) (C.10) (C.14) As stated previously, this approach allows for discussion of the physics underlying each contribution as well as their relative importance once the parallel mode structure and frequency of the modes are calculated. The decomposition introduced with equation (C.5) discriminates the contributions due to the kinetic compressibility δKs and the advective response of the plasma, while resonant and nonresonant responses are further identified by Fourier analysis, as we will see in the next paragraphs. As already stated in this section, these expressions for the fluxes are general and admit, as a particular case, the quasilinear limit if the linear particle response for δKs is adopted, i.e. the solution of equation (C.6). However, if the solution of the nonlinear gyrokinetic equation is used, the corresponding gyrokinetic fluxes are fully consistent with the results of a global gyrokinetic/hybrid code [76] when the self-consistent radial envelope A(r, t) is assumed. Recently, an intermediate approach has been introduced in [4,21], the socalled DSM, where the most important nonlinear effects are retained in the description of radial envelope, while the parallel mode structure is given by the solution of linear gyrokinetic equations. As anticipated, these three levels of simplification of the gyrokinetic fluxes within a unified theoretical framework constitute a hierarchy of reduced descriptions for plasma transport that can be used to develop advanced reduced transport models and for validation and verification purposes [87].
In this final part of appendix B, without losing generality, we focus on the calculation of the quasilinear fluxes that are evaluated numerically in section 5. As already stated, we proceed by substituting the solution of (C.6) inside equation (C.14), thus obtaining the complete expression for the fluctuation-induced fluxes. As an example, we show here one of the contributions due to the kinetic compression responsê dϑ r 3/2 √ r − λφ ⋆ (ϑ)δKsn(ϑ − 2π ℓ) whereΦ ∥ =φ −ψ. For the sake of notation simplicity, we can rewrite the preceding equation aŝ dϑ Here, x l and xr represent the lower and upper boundaries in the extended η angle domain that is used to cover the ballooning space in the actual calculation. By analyzing the flux in the Fourier space, we can immediately show, as expected, that the kinetic compression contribution is deeply related with the resonant particles. In fact, focusing on the ∝Φ ′ ∥ integrand, as the other contributions can be handled in a similar fashion dϑ