Shear Alfv\'en and acoustic continuum in general axisymmetric toroidal geometry

TL;DR

This paper develops a general framework for analyzing shear Alfvén and ion sound wave continua in axisymmetric toroidal plasmas, extending previous models to include complex geometries and equilibrium effects, with applications to tokamak stability.

Contribution

It introduces a novel formulation in ballooning space that extends continuum spectrum analysis to include gyrokinetic and 3D effects, building on and generalizing earlier work.

Findings

Excellent agreement with MARS code results

Continuum structures linked to plasma inertia

Method applicable to complex 3D geometries

Abstract

The equations describing the continuous spectrum of shear Alfv\'en and ion sound waves propagating along magnetic field lines are introduced and solved in the ballooning space for general geometry in the ideal MHD limit. This approach is equivalent to earlier analyses by Chu et al. 1992 [Phys. Fluids B 4, 3713 (1992)] but the present formulation in the ballooning space allows to readily extend it to include gyrokinetic and three-dimensional equilibrium effects. In particular, following Chen and Zonca 2017 [Phys. Plasmas 24, 072511 (2017)], the MHD limit is adopted to illustrate the general methodology in a simple case, and the equations are solved within the framework of Floquet and Hill's equation theory. The connection of shear Alfv\'en and ion sound wave continuum structures to the generalized plasma inertia in the general fishbone like dispersion relation is also illustrated and discussed. As an application, the continuous frequency spectrum is calculated for a reference equilibrium of the Divertor Tokamak Test facility. The results are compared with those obtained by the MARS code adopting the standard methodology, demonstrating excellent agreement.

Equations

$${\mathbf{B}}_0=F(\psi )\nabla \phi +\nabla \phi \times \nabla \psi ,$$

$$q=\frac{{\mathbf{B}}_0\cdot \nabla \zeta }{{\mathbf{B}}_0\cdot \nabla \theta }=q(r)$$

$$\delta {\mathbf{\xi }}_{\perp }=\frac{c}{B_0}\mathbf{b}\times \nabla {\Phi }_s,$$

$${\Phi }_s(r,\theta ,\zeta )=\sum _m\exp (in\zeta -im\theta ){\Phi }_{sm}(r)$$

$${\Phi }_s(r,\theta ,\zeta )$$

$$\left({\partial }_{\vartheta }^2-\frac{{\partial }_{\vartheta }^2|\nabla r|}{|\nabla r|}+\frac{{\omega }^2{J}^2{B}_0^2}{v_A^2}\right)y_1$$

$$\left(1+\frac{c_s^2}{{\omega }^2}\frac{1}{J^2{B}_0^2}{\partial }_{\vartheta }^2\right)y_2$$

$$y_1\equiv \frac{{\hat{\varphi }}_s}{{\left(\overline{\beta }{q}^2\right)}^{1/2}}\frac{c{k}_{\vartheta }}{{\bar{B}}_0{R}_0},y_2\equiv i\frac{\delta {\hat{P}}_{comp}}{(2\Gamma {)}^{1/2}{P}_0},$$

$${\kappa }_g=-\frac{1}{J{B}_0}\frac{F}{|\nabla \psi |}\frac{\partial }{\partial \vartheta }\ln B_0,$$

$$\eta \equiv 2\pi \frac{{\int }_0^{\vartheta }d{\vartheta }^{\prime }|\nabla r{|}^{-2}}{{\int }_0^{2\pi }d{\vartheta }^{\prime }|\nabla r{|}^{-2}}.$$

$$|\nabla r{|}^2\mathcal{C}{\frac{\partial }{\partial \vartheta }|}_{\psi }={\frac{\partial }{\partial \eta }|}_{\psi },$$

$$J_{\eta }=J\mathcal{C}|\nabla r{|}^2.$$

$$g_1=\frac{1}{|\nabla r|}{y}_1,g_2=\frac{c_s}{\omega q{R}_0|\nabla r|}\frac{|s\vartheta |}{s\vartheta }{y}_2,$$

$$\left({\partial }_{\eta }^2+\frac{{\omega }^2{J}_{\eta }^2{B}_0^2}{v_A^2}\right)g_1$$

$$\left({\partial }_{\eta }^2-|\nabla r|{\partial }_{\eta }^2|\nabla r{|}^{-1}+\frac{{\omega }^2{J}_{\eta }^2{B}_0^2}{c_s^2}\right)g_2$$

$${\kappa }_g=-\frac{1}{J_{\eta }{B}_0^2}\frac{F}{|\nabla \psi |}\frac{\partial }{\partial \eta }{B}_0.$$

$${\mathbf{x}}_i=e^{i{\nu }_i\eta }{\mathbf{P}}_i(\eta ),$$

$${\nu }_i={\nu }_i(\omega ,r).$$

$${\nu }_i^2(\omega ,r)=(nq(r)-m{)}^2.$$

$$g_1^{(i)}(\eta ;{\nu }_i,r)$$

$$g_2^{(i)}(\eta ;{\nu }_i,r)$$

$$\int d\eta |{\hat{g}}_1^{(i)}{|}^2=\int d\eta |{\hat{g}}_2^{(i)}{|}^2=1,$$

$$|e_1^{(i)}({\nu }_i,r){|}^2+|e_2^{(i)}({\nu }_i,r){|}^2=1.$$

$$i|s|\Lambda =\delta {\hat{W}}_f+\delta {\hat{W}}_k,$$

$$\mathbf{x}=w_1{\mathbf{x}}_1+w_2{\mathbf{x}}_2,$$

$$i\Lambda ={\frac{\left(i{\nu }_1{P}_1(\eta )+P_1^{\prime }(\eta )\right)w_1+\left(i{\nu }_2{P}_2(\eta )+P_2^{\prime }(\eta )\right)w_2}{P_1(\eta )w_1+P_2(\eta )w_2}|}_{\eta =2\ell \pi },$$

$$\mathbf{X}(\vartheta +2\pi )=\mathbf{X}(\vartheta )\mathbf{M}.$$

$$\mathbf{M}={\mathbf{X}}^{-1}(0)\mathbf{X}(2\pi ).$$

$$\mathbf{M}=\mathbf{X}(2\pi ).$$

$${\rho }_i=e^{i2\pi {\nu }_i}.$$

$$\frac{d^{2}x}{d{\eta }^2}+V(\eta ;\omega )x=0,$$

$${\mathbf{x}}_i(\eta +2\pi ;\omega )={\rho }_i{\mathbf{x}}_i(\eta ;\omega ),$$

Full text

I Introduction

Computing the shear Alfvén wave (SAW) continuous frequency spectrumBarston (1964); Sedlác̆ek (1971); Grad (1969); Uberoi (1972); Tataronis and Grossman (1973); Hasegawa and Chen (1974); Chen and Hasegawa (1974a, b); Dewar et al. (1974); Appert, Gruber, and Vaclavik (1974); Goedbloed (1975); Pao (1975); Chance et al. (1977) in general toroidal geometry is important because it determines mode structures and dispersive properties of Alfvénic fluctuations excited by energetic particles (EPs) in fusion devices, which, in turn, play crucial role in determining the plasma collective behaviors and eventually may affect fusion performanceChen and Zonca (2016). Due to equilibrium magnetic field curvaturePogutse and Yurchenko (1978); D’Ippolito and Goedbloed (1980); Kieras and Tataronis (1982); Cheng, Chen, and Chance (1985); Cheng and Chance (1986), SAW and ion sound wave (ISW) continuous frequency spectrum are coupled and, thus, their structures should be considered self-consistently, taking into account realistic equilibrium magnetic field geometry and plasma nonuniformity.

In numerical applications adopting realistic equilibrium reconstruction, the typical approach to the coupled SAW and ISW spectra is to Fourier decompose the fluctuation structure in general poloidal and toroidal angles and to compute the null space (kernel) of the matrix of the highest order radial derivativeCheng and Chance (1986), e.g., of the radial plasma displacement. This is completely equivalent to solve for the radial singular structures of the coupled SAW and ISW propagating along the field lines, described by two second order coupled differential equationsChu et al. (1992). In particular, the latter method has been adopted to compute the SAW and ISW spectra in DIII-DChu et al. (1992), focusing on both shaping and finite pressure effects on the continuous spectrum structures at low frequencyChu et al. (1992); Huysmans et al. (1995); Goedbloed (1998); Van der Holst, Beliën, and Goedbloed (2000). If the ideal MHD plasma description is adopted, again, the two approaches described above are mathematically equivalent, although it has been pointed out in Ref. 20 that the second method has advantages when computing the high toroidal mode number continuous spectra.

Kinetic description becomes necessary at short wavelengths and/or low frequencies to account for finite parallel electric field and wave damping. For the SAW-ISW coupling at low frequency and corresponding frequency gapsChu et al. (1992); Huysmans et al. (1995); Goedbloed (1998); Van der Holst, Beliën, and Goedbloed (2000) with possible discrete modes localized thereinGorelenkov et al. (2007a, b, 2009), the need of kinetic theory for the damping assessment was pointed out in Refs. 27; 28; 29 and motivated gyrokinetic stability studies of low-frequency Alfvénic fluctuations exctited by EPs with realistic geometry and plasma profilesZhang et al. (2016); Liu et al. (2017); Bierwage and Lauber (2017). The radial singular structures characterizing the SAW and ISW continuous spectra are most easily isolatedPegoraro and Schep (1986); Newcomb (1990) adopting the ballooning mode representationConnor, Hastie, and Taylor (1978), or the mode structure decomposition approachLu, Zonca, and Cardinali (2012a), which applies for arbitrary mode numbers and reduces to the ballooning representation in the high mode number limit. Within this framework, the general description of SAW and ISW using gyrokinetic theory has been given in Refs. 14; 37; 38, and allows to include thermal plasmaZonca, Chen, and Santoro (1996); Zonca et al. (1998, 1999); Lauber (2013); Bowden, Hole, and Könies (2015); Lauber and Lu (2018) as well as EPs inducedZonca and Chen (1996); Zonca et al. (1999); Lauber (2013); Zonca and Chen (2014b); Bowden, Hole, and Könies (2015) kinetic effects in the analysis of the continuous spectra.

This premise provides a motivation of the present work, which is to show that the coupled second order differential equations describing SAW and ISW propagation along equilibrium magnetic field lines for general geometry are most easily expressed and solved in the ballooning space. In the ideal MHD limit, as already stated, this is equivalent to the approach adopted in Ref. 20. The present methodology, however, is readily extended to include gyrokinetic descriptionChen and Zonca (2016); Zonca and Chen (2014a, b) and 3D geometryDewar and Glasser (1983). Furthermore, this formulation not only allows us to discuss the structures of SAW and ISW continuous spectra in general toroidal geometry, but it provides the simplest and most direct way of computing the generalized inertia as an ingredient of the general fishbone-like dispersion relation (GFLDR)Chen and Zonca (2016); Zonca and Chen (2014a, b); that is, the unified framework for describing Alfvénic fluctuations excited by EPs in tokamaks. In Sec. II, in order to establish a link to prior work more closely, we take the ideal MHD limitChen and Zonca (2017) of the general descriptionChen and Zonca (2016); Zonca and Chen (2014a, b) and we write the coupled SAW and ISW equations for radial singular structures in general tokamak geometry. Along with providing a synthetic summary of prior theoretical analyses, which are applied here to realistic plasma equilibria, we illustrate the link of the generalized inertia in the GFLDR to the structures of the continuous spectra and show how it is computedChen and Zonca (2016); Zonca and Chen (2014a, b); Chen and Zonca (2017). Furthermore, the important role of polarization of physical fluctuations is analyzed with emphasis on its relevance for assessing their effective absorption by resonantly excited radial singular structures of the continuous spectrum. We then analyze the general form of the obtained equations, collocating that within the framework of Floquet theoryFloquet (1883) and Hill’s equationHill (1886); Magnus and Winkler (2013); Denk (1995). As original application, in Sec. III we solve for the SAW and ISW continuous spectra for a reference equilibrium of the Divertor Tokamak Test (DTT) facilityAlbanese et al. (2017, 2019). We discuss both (artificially) decoupled and (realistic) coupled cases to illustrate features of the obtained SAW and ISW spectra and directly compare results of the present approach with those obtained by the MARSBondeson, Vlad, and Lütjens (1992) code, which is equivalent to computing the null space of the matrix of the highest order radial derivativeCheng and Chance (1986) of the plasma radial displacement. In particular, we are able to identify the role of various equilibrium geometry effects in SAW-SAW, ISW-ISW and SAW-ISW couplings and corresponding structures in the continuous spectra. All these detailed features are naturally accounted for in the present formulation and can be included as boundary conditions for calculating parallel mode structures and GFLDRZonca and Chen (2014a, b) dispersion relations in general Tokamak geometry. Doing so, however, is beyond the scope of the present work and will be done elsewhere. Further extensions to gyrokinetic analyses and 3D geometry are discussed in Sec IV, where concluding remarks are provided along with considerations about future perspectives.

II Theoretical framework

The general theoretical framework is discussed in Ref. 14, based on the detailed derivations given by Zonca and Chen in Refs. 37; 38. There, the general equations governing SAW and ISW fluctuation in low-$\beta$ magnetized fusion plasmas are given by quasineutrality condition and vorticity equation, which are derived from the gyrokinetic equationFrieman and Chen (1982). Meanwhile, Ref. 14 also provides a thorough discussion of how gyrokinetic vorticity equation and quasineutrality condition yield the proper reduced MHD limit when appropriate Brizard (1994); Qin, Tang, and Rewoldt (1998, 1999). The reduced MHD limiting form is also derived by the standard approach in Ref. 47. In the present application to general tokamak equilibria, we merely assume the aforementioned theoretical framework without derivation and only provide the essential elements for introducing the notation and solving the relevant equations. The general features of the mode structure decomposition of fluctuations in toroidal plasmasLu, Zonca, and Cardinali (2012a) are also discussed in Refs. 37; 38 along with its usefulness to generally describe continuous spectra.

II.1 Notation and fundamental equations

In this work we assume an axisymmetric equilibrium magnetic field $\bm{B}_{0}$ expressed in flux coordinates $(r,\theta,\varphi)$:

[TABLE]

where $r(\psi)$ is a radial-like flux coordinate, $\varphi$ is the physical toroidal angle and the angular coordinate $\theta$ can be chosen such that the Jacobian $J=(\bm{\nabla}\psi\times\bm{\nabla}\theta\cdot\bm{\nabla}\varphi)^{-1}$ has a convenient expression. Following the common practice, we also introduce the straight field line toroidal angle $\zeta$, defined such that the safety factor

[TABLE]

is a flux function. We introduce the leading order perpendicular plasma displacement as usual,

[TABLE]

with $\Phi_{s}$ being the perturbed stream functionChen and Zonca (2017). The fluctuating field

[TABLE]

expressed as a Fourier series, where $m$ and $n$ represent the poloidal and toroidal mode numbers, respectively, can be decomposed asLu, Zonca, and Cardinali (2012a):

[TABLE]

Here, time dependences are left implicit for simplicity of notation, and $\vartheta$ represents the extended poloidal angleConnor, Hastie, and Taylor (1978) coordinate following equilibrium magnetic field lines. Periodic poloidal angle-$\theta$ dependences of equilibrium quantities are replaced by the same periodic dependences on $\vartheta$. However, while fluctuations must be periodic in $\theta$ space for physical reasons, they are generally not periodic in $\vartheta$, as shown in Eq. (4). The radial singular structures corresponding to the continuous spectra are obtained from the limiting forms of, respectively, vorticity and pressure equations for $|\vartheta|\rightarrow\infty$Zonca and Chen (2014a, b); Chen and Zonca (2017). Following Refs. 14; 37; 38, it is possible to derive the expressions describing the structure of the continuum modes in flux coordinates. In particular, assuming $JB_{0}^{2}$ constant on a flux surface, yielding the so called “Boozer coordinates”, we obtainChen and Zonca (2017):

[TABLE]

where

[TABLE]

$\hat{\phi}_{s}(r,\vartheta)\equiv|s\vartheta||\bm{\nabla}r|\hat{\Phi}_{s}(r,\vartheta)$, $s=rq^{\prime}/q$ is the magnetic shear, $k_{\vartheta}=-nq/r$, and $\delta\hat{P}_{comp}(r,\vartheta)$ is the representation of the compressional component of the pressure perturbation according to Eq. (4)Chen and Zonca (2017). Furthermore, $\kappa_{g}$ is the geodesic curvature in Boozer coordinates:

[TABLE]

$\bar{\beta}=8\pi\Gamma P_{0}/\bar{B}_{0}^{2}$, $\bar{B}_{0}$ is the magnetic field on axis, $r=a\rho_{tor}$, where $\rho_{tor}$ is defined, as usual, as a function of the normalized toroidal magnetic flux and $a$ is the minor radius of the torus. Note, again, that periodic equilibrium field dependences on $\theta$ are replaced by the same periodic dependences on $\vartheta$ in the ballooning space.

These equations describe SAWs and ISW continuous spectra coupled by equilibrium non-uniformities and geodesic curvature. Both type of waves propagate along equilibrium magnetic field lines and, therefore, radial equilibrium nonuniformities are accounted for by $r$, which enters in Eqs. (5) as a parameter. Note also that the continuous spectrum does not depend on magnetic shear although it formally appears in Eqs. (5) and (6). This can be understood since the continuous spectrum is a local feature of singular radial fluctuation structures. More in depth discussion about this point and what happens for vanishing magnetic shear can be found in Refs. 14; 37; 38. Equations (5), obtained in the MHD limitChen and Zonca (2017), are valid for general tokamak equilibrium geometry and can be used to calculate the coupled SAW and ISW continuous spectra for cases of practical interest, such as the DTT reference scenarioAlbanese et al. (2017, 2019) introduced in Sec. III.2. These equations are readily extended to 3D plasma equilibria, such as stellarators, thanks to the generality of ballooning representationDewar and Glasser (1983). The same formulation can be adopted within a linear gyrokinetic description of coupled SAW and ISW continua, in which case Eqs. (5) are replaced by the short radial scale gyrokinetic vorticity equation, i.e., Eq. (6) of Ref. 38, and the quasineutrality condition, i.e., Eq. (36) of the same work. This theoretical framework is similar to the description of electron wave-packets in a 1D periodic lattice, see e.g. Ref. 59. The interesting analogyChen and Zonca (2016) is due to the fact that the governing equations of low-frequency plasma waves in magnetized plasmas describe the propagation of wave-packets along magnetic field lines with periodicity given by the plasma connection length. Therefore, this correspondence is independent of the axisymmetry of the equilibrium magnetic field and/or the fluid/kinetic nature of the theory. In particular, Eqs. (6), are Schrödinger-like and produce structures analogue to electronic bands in solid-state physics. Further discussion of this interesting analogy is given in Refs. 14; 37; 38.In the following, we will focus on the solution of Eqs. (5). By direct inspection, it is readily noted that, adopting Boozer coordinates, ISW continuous spectrum is more simply represented than SAW continuous spectrum, due to the modulation term $\propto\partial_{\vartheta}^{2}|\bm{\nabla}r|/|\bm{\nabla}r|$. Since there is usually more interest in the SAW than in the ISW branch due to the stronger kinetic damping of the latter, we can rewrite Eqs. (5), introducing an appropriate angular coordinate to express the SAW frequency continuum more transparently; i.e.,

[TABLE]

By construction, the Jacobian $J_{\eta}$ of this new set of coordinates, dubbed continuum flux coordinates (CFC), is such that $J_{\eta}B_{0}^{2}/|\bm{\nabla}r|^{2}$ is a flux function. We obtain the following expression for the derivative along $\eta$ (at constant $\psi$):

[TABLE]

where $\mathcal{C}=(2\pi)^{-1}\int_{0}^{2\pi}d\vartheta^{\prime}\left|\bm{\nabla}r\right|^{-2}$. The Jacobians in the two set of coordinates are related by the following expression:

[TABLE]

Introducing the normalizations:

[TABLE]

we can re-write Eqs. (5) as

[TABLE]

respectively, where the geodesic curvature in this new set of coordinates reads:

[TABLE]

Equations (12) are completely equivalent to Eqs. (5) and will be solved in the following, illustrating the properties of coupled SAW and ISW continuous frequency spectra in general tokamak plasma equilibria. They are in the form of two coupled second order differential equations with periodic coefficients. Therefore, all known results of Floquet theoryFloquet (1883) can be applied. These are summarized in the remaining part of Sec. II for the reader’s convenience.

II.2 Linear differential equations with periodic coefficients

II.2.1 Floquet theory

Equations (12) are a linear system of second order differential equations with periodic coefficients. Therefore, after re-writing it as four first order differential equations, we know from Floquet theoryFloquet (1883) that it must have solutions in the form:

[TABLE]

where $\mathbf{P}_{i}$ is a $2\pi$-periodic function, $i=1,2,3,4$ and the $\nu_{i}$ are complex numbers called characteristic Floquet exponents of Eqs. (12). Finding the solutions of Eqs. (12) for a given value of $r$ and, therefore, calculating $\nu_{i}$ for each $\omega$ value is equivalent to obtaining the (local) dispersion curves of the system:

[TABLE]

We emphasize that this relation involves only local quantities and describes wave packets propagating along magnetic field on a given flux surface. Recalling the physical meaning of the derivative along the $\vartheta$ coordinateChen and Zonca (2016); Zonca and Chen (2014a, b), we can relate the characteristic Floquet exponents to $\nabla_{\parallel}\equiv(JB_{0})^{-1}\partial_{\vartheta}=(J_{\eta}B_{0})^{-1}\partial_{\eta}$, and, in particular, to the toroidal and poloidal mode numbers $n,m$:

[TABLE]

Equation (16) readily follows from the definition of $\nabla_{\parallel}$ and the mode structure decomposition of Eq. (4)Chen and Zonca (2016); Zonca and Chen (2014a, b); Chen and Zonca (2017). Therefore, integrating Eqs. (12) for different $\omega$ values and using Eqs. (15) and (16), it is possible to calculate the continuous frequency spectrum for every value of the toroidal mode number. In particular, Eq. (16) explicitly shows that continuous spectra can be computed for arbitrary $n$, once the dispersion curves $\nu_{i}=\nu_{i}(\omega,r)$ are given. Extensions to 3D/stellarator geometry and gyrokinetic theory would be obtained in the same way from the more general governing equations as briefly discussed in Sec. II.1.

For given $\nu_{i}(\omega,r)$; i.e., for given dispersion curves of Eqs. (12), it is possible to calculate the corresponding vector solutions; i.e., eigenvectors of the fundamental matrix $\mathbf{X}(\vartheta)$, that is a matrix-valued function whose columns are linearly independent solutions of Eqs. (12). Furthermore, since the considered problem is linear, we can generally write the normalized eigenvector solution corresponding to $\nu_{i}(\omega,r)$ as

[TABLE]

Here, without loss of generality, we can impose

[TABLE]

where the $\eta$-space integration domain may be chosen as any convenient multiple of $2\pi$ due to the periodicity of Eqs. (12). Meanwhile,

[TABLE]

Thus, the $(e_{1}^{(i)},e_{2}^{(i)})$ vector can be considered as definition of fluctuation polarization in the linear (Hilbert) space of Floquet solutions $[g_{1}^{(i)}(\eta;\nu_{i},r),g_{2}^{(i)}(\eta;\nu_{i},r)]$, corresponding to continuous spectrum fluctuations for a given $\nu_{i}$. Polarization is a necessary element that uniquely identifies the continuous spectrum, described by solutions of Eqs. (12), together with the characteristic Floquet exponent, $\nu_{i}(\omega,r)$, and the corresponding parallel fluctuation structures, $\hat{g}_{1}^{(i)}(\eta;\nu_{i},r)$ and $\hat{g}_{2}^{(i)}(\eta;\nu_{i},r)$. A pictorial representation of parallel propagation of radial singular structures of the continuous spectrum and of their polarization in the orthogonal Hilbert space is given in FIG.1.

Here, it is also worthwhile pointing out the connection of the Floquet characteristic exponent to the generalized plasma inertia, $\Lambda$, which is an essential element of the general fishbone-like dispersion relation (GFLDR)Chen and Zonca (2016); Zonca and Chen (2014a, b):

[TABLE]

where $\delta\hat{W}_{f}$ and $\delta\hat{W}_{k}$ account for, respectively, the generalized potential energy due to fluid-like plasma response and the kinetic plasma behavior due to, e.g., EPs. At a given radial location, the short scale (large-$\eta$) radial structure of any considered physical fluctuation (antenna driven and/or eigenmode) satisfies the same Eqs. (12) describing the continuous spectrum. Thus, the large-$|\eta|$ local physical solution must be a linear superposition of the solutions of Eq. (14). Among them, generally, only two satisfy boundary conditions (outgoing/decaying wave in $\eta$-spaceChen and Zonca (2016); Zonca and Chen (2014a, b)). Renumbering them as $i=1,2$, without loss of generality, we have

[TABLE]

where $w_{1}$ and $w_{2}$ are weights of the linear combination, which we can generally assume satisfying the normalization condition $w_{1}^{2}+w_{2}^{2}=1$. The actual value of $(w_{1},w_{2})$ is determined consistently with the mode structure and fluctuation dispersion relation. Adopting the GFLDR theoretical frameworkChen and Zonca (2016); Zonca and Chen (2014a, b), it is possible to showZonca and Chen (2014a, b)

[TABLE]

where $P_{i}(\eta)|_{\eta=2\ell\pi}(i=1,2)$ is referred to the $g_{1}$ solution of Eqs. (12) represented as in Eq. (14) for $\nu_{i}(i=1,2)$, respectively, and, for simplicity, we have assumed that the equilibrium is general but still up-down symmetric. Since the continuum fluctuation structures, Eq. (17), bear the information of mode polarization, this relation shows the importance of mode polarization to assess the actual coupling of any physical fluctuation to the continuous spectrum; and further illuminates the profound connection between the structures of SAW and ISW continuous spectra and the GFLDR, and the usefulness of the present approach when extended to 3D/stellarator geometry and/or gyrokinetic descriptions. More on this important point is discussed in Section III below.

In order to calculate $\nu_{i}$ as a function of $\omega$ for each $r$, we note that the fundamental matrix $\mathbf{X}(\vartheta)$ of the first order system associated to Eqs. (12) must satisfy the following relation:

[TABLE]

It can be shownFloquet (1883) that $\mathbf{M}$ does not depend on $\vartheta$ and, therefore, is conveniently calculated for $\vartheta=0$:

[TABLE]

Choosing the initial conditions such that $\mathbf{X}(0)=\mathbf{I}$, the 4x4 identity matrix in the present case, we obtain:

[TABLE]

The eigenvalues of $\mathbf{M}$ are the characteristic multipliers of the system and can be expressed in terms of the characteristic Floquet exponents by the following relation:

[TABLE]

Corresponding eigenvectors can be used to calculate parallel mode structures and polarization vectors, Eq. (17); e.g., to calculate $P_{i}^{\prime}(2\pi)/P_{i}(2\pi)$ and the value of the generalized inertia from Eq. (22). The same information can be adopted as boundary conditions for the calculation of fluid and kinetic potential energies in the GFLDRChen and Zonca (2016); Zonca and Chen (2014a, b). However, doing so is outside the scope of the present work.

II.2.2 Hill’s equation

In the decoupled case ($\kappa_{g}=0$), Eqs. (12) become a system of two independent differential equations in the form:

[TABLE]

where $x=g_{1},g_{2}$; which are a particular cases of Hill’s equationHill (1886); Magnus and Winkler (2013); Denk (1995). Note that, consistent with previous notations, parametric dependences on $r$ are left implicit for notation simplicity, while the dependence on $\omega$ of the potential function $V(\eta;\omega)$ is indicated explicitly. Reducing each one of these equations to a first order system and using the results summarized in the previous subsection, we know that they will have Floquet solutions satisfying the following property:

[TABLE]

where $i=1,2$. We can show that the following relation holdsAbramowitz and Stegun (1965):

[TABLE]

with:

[TABLE]

where $\bm{x}_{10}$ and $\bm{x}_{01}$ are the particular solutions of the first order system associated to each Hill’s equation such that $\bm{x}_{10}(0)=(1,0)$ and $\bm{x}_{01}(0)=(0,1)$. In the particular case of even $V(\eta;\omega)$, it can be shown that $A(\omega)=x_{10}(2\pi;\omega)$. Solving Eq. (7) gives the following expression for the multipliers:

[TABLE]

If $|A^{2}|<2$, roots are complex conjugates with unity absolute value while, if $|A^{2}=2|$, they coincide and are equal to $1$. Finally, if $|A^{2}>2|$, roots are real and reciprocal to each other. In particular, the first case describes the frequency continuum while the third the frequency gaps. From this expression is possible to calculate the characteristic exponent from the value of $A(\omega)$. In the first case we obtain:

[TABLE]

where $\nu$ is real. From this analysis we expect up to two (opposite) values of $\nu$ for each $\omega$ in the decoupled case, corresponding to waves propagating in opposite directions along the field lines. This result can be generalized to higher dimensionality, e.g. see Ref.51, to describe the coupled system, where we expect up to two (pairs of opposite) $\nu$ values for each $\omega$ corresponding, again, to waves propagating in opposite directions along the field lines.

III Numerical results

In this Section, in order to show the generality of the present approach, we calculate the frequency continuum of SAW and ISW waves in a Divertor Tokamak Test Facility (DTT) reference scenario. DTT is an Italian project Albanese et al. (2017, 2019) sponsored by the EUROfusion consortium with the goal of designing a new machine capable of eventually integrating all relevant physics and technological issues concerning power exhaust solutions for DEMO. The main DTT parameter are reported in FIG.2. In order to isolate the role of various equilibrium geometry effects in SAW-SAW, ISW-ISW and SAW-ISW couplings and discuss the corresponding structures in the continuous spectra, we illustrate first the decoupled SAW and ISW continuous spectra by artificially letting $\kappa_{g}=0$ (Section III.3). The realistic situation with coupled SAW and ISW spectra is then addressed in Section III.4. Numerical results obtained with the present approach are compared with those by the MARS code adopting the conventional method, showing excellent agreement.

III.1 Normalized equations

The set of equations describing SAWs coupled with ISW can be re-written in terms of the following dimensionless quantities:

[TABLE]

where $\overline{v}_{A0}$ is the Alfvén velocity on the magnetic axis. The non-dimensional form of Eqs. (12) then reads now

[TABLE]

Eqs. (34) are studied numerically in the remaining part of this work for the DTT reference scenario illustrated in the following subsection.

III.2 DTT reference scenario

The magnetic equilibrium has been originally calculated by means of the free bounddary equilibrium evolution code CREATE-NL Albanese et al. (2003) and further refined and mapped on flux coordinates by using the high-resolution equilibrium solver CHEASE Lütjens, Bondeson, and Sauter (1996). We consider a double null setting whose basic profiles are depicted in FIG.3 as a function of the normalized toroidal radius $r/a$. We note that, for the analyzed case, the kinetic pressure on axis is $1.0768\times 10^{6}$Pa, the flux function $F$, see Eq.(1), on axis is $12.45$Vs, while the density is $10^{19}$m*-3*.

Moreover, the normalized density profile has been obtained by studies of plasma scenario formation using the fast transport simulation code METIS Artaud et al. (2018). For the addressed scenario, the $(\psi,\vartheta)$ and $(\psi,\eta)$ isosurfaces are depicted in FIG.4.

III.3 SAW continuous spectrum

The equation describing SAW continuous spectrum can be obtained from Eqs. (34) by neglecting the coupling terms with the ISW branch. Thus, we have

[TABLE]

In order to illustrate with a practical application the procedure introduced in Sec. II.2.1, we first show the local dispersion curves $\nu(\Omega)$ obtained integrating Eq.(35) at fixed $r/a$, i.e., on a given $\psi$ isosurface, for different $\Omega$ values. These curves are plotted in FIG.5 for the reference value $r/a=0.58$, where, for convenience, $\nu$ is chosen as abscissa showing different branches and gaps.

Noting Eq. (4) and assuming a Fourier decomposition of the fluctuations mapped back from ballooning to real spaceLu, Zonca, and Cardinali (2012a), we can write

[TABLE]

Thus, two counter-propagating SAWs can form a standing wave and cause the formation of a frequency gap in the continuous spectrumChance et al. (1977); Cheng, Chen, and Chance (1985); Zonca and Chen (1992, 1993), when the Bragg reflection condition is satisfied; that is:

[TABLE]

Due to the periodicity of Floquet solutions, Eq. (14), Eq. (15) can be effectively reduced to the first Brillouin zone, $0\leq\nu\leq 1/2$, and frequency gaps are expected to appear atZonca and Chen (2014a, b) (cf. Appendix A)

[TABLE]

with

[TABLE]

consistent with Eqs. (36) and (37). As anticipated above, dispersion curves $\nu=\nu(\Omega,r/a)$, i.e., the dimensionless form of Eq. (16), include all the relevant information for the construction of continuous spectra. Using Eq.(14) and combining the results of different flux surfaces, we obtain the SAW continuous spectrum as a function of $r/a$ as depicted in FIG.6, where continua for various toroidal mode numbers are shown in different colors and illustrate the formation of toroidicity- as well as ellipticity-induced frequency gaps at, respectively, $\Omega\simeq 1/2$ (TAECheng and Chance (1986)) and $\Omega\simeq 1$ (EAEBetti and Freidberg (1991, 1992)). Higher frequency gaps could also be readily plotted by means of the dispersion curves reported in FIG.5 but are omitted here to illustrate results more clearly.

We now compare the results with the spectrum calculated by MARS codeBondeson, Vlad, and Lütjens (1992). The continuous spectrum is obtained assigning a guess real frequency $\omega_{c}=\omega_{\mathrm{guess}}$ and looking for converged solutions at $\omega=\omega_{\mathrm{c}}$, which correspond to a perturbed velocity whose toroidal contravariant $(m,n)$ Fourier components $v^{\phi}_{m,n}(r)$ exhibit a divergent radial dependence of the type $v^{\phi}_{m,n}(r)\propto 1/(r-r_{\mathrm{c}})$. When this condition is recognized, the points $(r_{\mathrm{c}}/a,\omega_{\mathrm{c}})$ in the plane $(r/a,\omega)$ are considered as the location of the continuous spectrum. The comparison for the toroidal mode number $n=2$ is depicted in FIG.7, showing good agreement between the two approaches.

As anticipated above, one of the peculiar features of the approach proposed by the present work is that the dispersion curves are independent of $n$ and, therefore, the evaluation of the continuous spectrum for different toroidal mode numbers merely requires the numerical solution of an algebraic equation, Eq. (16). Another interesting property of the present methodology is the freedom in the choice of the radial mesh due to the (radial) local nature of the problem that we are solving. This feature can be exploited to more accurately describe regions where the continuous spectrum is more convoluted. Finally, the present approach is readily extended to 3D geometries and kinetic analyses, as noted in the Introduction and Section II.

The same analysis can be repeated for the ISW continuum. In fact, Eq. (34b) can be cast in the same form as Eq. 35 using Boozer coordinates and normalizing the frequency to the local sound speed rather than the Alfvén speed on axis. Results, thus, would be exactly the same as in FIG.6 upon a suitable rescaling of ISW continuum frequencies. In particular, a toroidicity induced gap would appear for $p=1$ in Eq. (38); an ellipticity induced gap would occur for $p=2$ and so on. These frequency gaps have not received significant attention in the literature since, as already pointed out, the ISW branch is typically affected by strong Landau dampingChavdarovski and Zonca (2009); Zonca et al. (2010); Chavdarovski and Zonca (2014); Liu et al. (2017); Zhang et al. (2016); Bierwage and Lauber (2017); Lauber (2013) and, therefore, discrete modes that are possibly located in these gaps are characterized by high excitation thresholdsChen and Zonca (2016); Zonca and Chen (2014a, b); Chen and Zonca (2017). Meanwhile, the features of ISW-ISW and SAW-SAW couplings due to equilibrium geometry are identical, as are the corresponding structures of continuous spectra. Therefore, a detailed calculation of the uncoupled ISW continuum is omitted here.

III.4 Coupled system

In the general case, we cannot neglect the SAW-ISW coupling due to equilibrium non-uniformities and we need to study the coupled system of Eqs. (34a) and (34b). In order to calculate the structures of SAW and ISW continuous spectra, we apply the procedure described in Section II. For the reference value $r/a=0.58$, the result is depicted in FIG.8, reproducing the expected behavior.

A change of the topology emerges with respect to the uncoupled case in the vicinity of the intersection points between SAW and ISW dispersion curves. In particular, it is worthwhile noting that standing waves and corresponding frequency gap formation can occur not only at the Bragg reflection condition, Eq. (37), due to the interaction of two counter-propagating SAWs or ISWs (cf. green and red solid circles in FIG.5), but also when

[TABLE]

where $\hat{\beta}\simeq\bar{\beta}$ is rigorously defined in Eq. (45). In FIG.5, regions where Eq. (40) is satisfied are denoted by the full orange circles, corresponding to a forward propagating ISW, modulated by the equilibrium non-uniformity within the magnetic flux surface, forming a standing wave by interaction with a backward scattered SAW. Modulated forward propagating ISW can also interact with a forward propagating SAW, when

[TABLE]

In this case, however, a standing wave cannot be formed, as it is clearly shown by the dashed orange circles in FIG.5. The gap formation at low frequency for $0.1<\nu<0.2$ is due to the interaction denoted by the full orange circle and satisfying Eq. (40) for $p=1$.

Combining the results of different flux surfaces, we finally get the behavior of the continuous spectra as shown in FIG.9 for $n=2$ and in FIG.10 for different toroidal mode numbers.

The richness of the continuous spectrum at low frequencyChu et al. (1992); Huysmans et al. (1995); Goedbloed (1998); Van der Holst, Beliën, and Goedbloed (2000) is readily noted, with the frequency gaps due to SAW-ISW and/or ISW-ISW couplings, dubbed Beta induced Alfvén Acoustic Eigenmode (BAAE) gapGorelenkov et al. (2009, 2007a, 2007b), as well as SAW-ISW and/or SAW-SAW couplings, dubbed Beta induced Alfvén Eigenmode (BAE) gapChu et al. (1992); Heidbrink et al. (1993); Turnbull et al. (1993), shown in FIG.10. The complex structures of the continuous spectrum are further illustrated by adding $n=10$ and $n=20$ toroidal mode numbers in FIG.10. For example, the weak dependence of upper branch of the continuous spectrum near the BAAE gap in the vicinity of mode rational surface is a reflection of the weaker plasma response to excitations by either fluid or EP effectsChen and Zonca (2017). By direct comparison between plots of the continuous spectrum and local dispersion curves, it is possible to identify the most important physical processes. As an example, comparing FIG.8 and FIG.9, we could deduce that the first three gaps observed in FIG.9 with $r/a\sim 0.58$ are produced by the interaction of two ISWs. In fact, for $n=2$ and the given safety factor profile, we obtain $\nu\sim 0.5$ from Eq. (16). Furthermore, the wave poloidal mode numbers producing the gaps in FIG.8 could be obtained by means of Eqs. (40) and (41). All this information could be routinely included in continuum spectrum plots. Nonetheless, for the sake of clarity, in this work we prefer to show only the frequency plot. Note, however, that the presence of continuous spectrum structures in a certain region of the $(r/a,\Omega)$ plane does not automatically imply continuum damping by resonant excitationChen and Hasegawa (1974b, a); Hasegawa and Chen (1974) of those same structures. As anticipated in Section II, a crucial role is played by fluctuation polarization, Eq. (17), which must be computed self-consistently accounting for fluctuation structure and dispersion relation, often affected non-perturbatively by EPsChen and Zonca (2016); Zhang et al. (2016); Liu et al. (2017); Bierwage and Lauber (2017); Zonca and Chen (2014a, b); Lauber (2013); Chen and Zonca (2017). Without solving for self-consistent mode structure and dispersion relation, a qualitative information about the effective role of resonant excitation of continuum structures can still be obtained from $|e_{1}|^{2}$ adopted as Alfvénicity parameter

[TABLE]

With this definition, consistent with that adopted in Ref. 47 for quantifying mode polarization in the considered MHD limit, ${\cal A}\simeq 1$ for SAW continuum, since $g_{2}\sim{\cal O}(\beta^{1/2})g_{1}$ at low frequency. Meanwhile, noting $g_{1}\sim{\cal O}(\beta^{1/2})g_{2}$ for ISW, ${\cal A}\sim{\cal O}(\beta)$ for the acoustic continuum. In general, $\mathcal{A}$ can be calculated for each point of the local dispersion curve and, therefore, for each point of the continuous spectrum plot. In fact, this corresponds to evaluate not only the Floquet characteristic exponent at each radial location, but the corresponding parallel fluctuation structure, $g_{1}^{(i)}(\eta;\nu_{i},r)$ (playing the role of “eigenvector”), and the corresponding polarization

[TABLE]

consistent with Eqs. (17)-(19). This allows, in principle, to represent the Alfvénicity by means of a color-map in these graphs. A detailed analysis along these lines is beyond the scope of this work, which mainly aims at the illustration of the methodology and the generality of its applications. As a particular case, we have computed the polarization and Alfvénicity at $r/a=0.58$ for the local dispersion curves represented in FIG.8. For reference, we have assumed $n=2$ such that $\nu=0.477$ corresponds to the value on $nq(0.58)$ reduced to the first Brillouin zone (cf. Sec. III.3 ). FIG.11 illustrates the relevant continuum frequencies, for which we have computed polarization and the values of Alfvénicity, ${\cal A}$, according to Eqs. (42)-(43). These are, ${\cal A}=0.02,0.06,0.59,0.28,0.61,0.34$ and $0.98$, from low to high frequencies, respectively. The first two values are very small, in consistency with the acoustic nature of the fluctuation that is predominant as expected from the green marking of the circle denoting SSW-SSW interaction in FIG.8. The mixed polarization nature (orange circle, denoting SSW-SAW coupling) is also clearly indicated by the Alfvénicity, which grows for the third to the sixth branch. We could actually state that ${\cal A}$ is a proper quantitative definition of the mixed polarization in these cases. Finally, the highest frequency branch is clearly Alfvénic, as suggested by both ${\cal A}=0.98$ and the red circle in FIG.8 denoting SAW-SAW coupling. In general, a SAW polarized fluctuation, characterized by $(e_{1},e_{2})=(1,0)$ in the uncoupled case, has weak interaction with the acoustic polarized continuum, which has $(e_{1},e_{2})=(0,1)$ in the uncoupled limit.

Meanwhile, an ISW polarization weakly interacts with the SAW continuous spectrum. This point is of particular importance when inspecting the structures of continuous spectra approaching the plasma edge, where, due to the decreasing ratio of sound to Alfvén speed, high order ISW sidebands may be generated by strong shaping equilibrium modulation effects and interact with SAW. This interaction generates a thick web of continuous spectrum structures, shown in FIG.12 again for $n=2$, where they are compared with numerical results obtained by MARS, showing excellent agreement also in this case. Clearly, not all structures of this web are physically relevant, because of the important role of wave polarization. More detailed illustration of this point requires a separate analysis computing the actual Alfvén wave fluctuation spectrum in DTT for the reference scenario (see, e.g., Refs. 70; 71).

IV Conclusions and future perspectives

This work, as an application of the general theoretical framework discussed in Refs. 14; 37; 38, illustrates the calculation of continuous spectrum structures in realistic tokamak geometry for the MHD limiting case of governing equationsChen and Zonca (2017). The use of ballooning formalism allows the calculation of local dispersion curves, isolating the physics information concerning the plasma response to radial local singular fluctuations, constituting the continuous spectrum. Meanwhile, toroidal and poloidal mode numbers are treated explicitly for arbitrary values, simplifying the numerical calculation of the complicated continuous spectrum structures at high mode number and/or near the plasma boundary. The present approach clearly illustrates the importance of fluctuation polarization for the assessment of damping by resonant absorption of radial singular continuous spectrum structures. As a qualitative estimate of the coupling strength of Alfvénic fluctuations to the acoustic continuum, without solving for the actual fluctuation structures, an Alfvénicity parameter is introduced, which is related to the fluctuation polarization and can be usefully adopted to visualize the importance of acoustic couplings.

Besides the simplicity of the present approach, based on the solution of two coupled second order ordinary differential equations for the local (singular) plasma response, the interesting aspect is that the present methodology could be readily extended to 3D plasma equilibria, such as stellarators, thanks to the generality of ballooning representationDewar and Glasser (1983). Following the the general theoretical framework of Refs. 14; 37; 38, the present formulation could be employed within a linear gyrokinetic description of coupled SAW and ISW continua. Relevant equations, in this case, would be the short radial scale gyrokinetic vorticity equation (Eq. (6) of Ref. 38) and the quasineutrality condition (Eq. (36) of the same work). Although the uderlying equations are different, the fundamental structure would be that of two second order periodic ordinary differential equations for the local (singular) plasma response, to be analyzed according to the well-known Floquet theory.

Within the present framework, the calculation of continuous spectrum, corresponding parallel fluctuation structures and polarization vectors provides the boundary conditions for computing any physical mode structures and corresponding dispersion relations in arbitrary geometry.

Acknowledgment

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

We thank DTT Team for the reference scenario.

Appendix A Derivation of useful integral quadratic forms

When we take the decoupled ($\kappa_{g}=0$) limit of the system of differential equations, Eqs. (34a) and (34b), it is possible to derive the following integral quadratic forms

[TABLE]

In the neighborhood of a given point of intersecting SAW and ISW continua, identified by a given $\nu(\Omega)$ for fixed $r/a$, we can write

[TABLE]

and, similarly,

[TABLE]

with $p\in\mathbb{Z}$. Here, we have noted that, near the given point of intersecting SAW and ISW continua, the $g_{1}$ and $g_{2}$ solutions are modulated by the equilibrium non-uniformity within the flux surface. The first of Eqs. (44) is what allows deriving Eq. (38), while the combination of the two is what readily yields Eqs. (40) and (41), with the following definition of $\hat{\beta}$:

[TABLE]