Shaping effects on the geodesic acoustic mode in tokamaks

The geodesic acoustic mode (GAM) is investigated with gyro-kinetic equations in Miller local equilibrium model for shaped tokamak plasmas with an arbitrary elongation κ and a finite triangularity δ. In particular, the effects of triangularity on GAM frequency and damping rate are analyzed both analytically and numerically. The asymptotic analytical and exact numerical results both show that the frequency almost linearly increases with the triangularity but increases relatively more slowly for a negative δ, which agree well with the TCV observation on the trend. The analytical results clearly claim that the triangularity effect strength is dependent on the inverse aspect ratio ϵ and Shafranov shift gradient Δ′ , while the numerical results indicate that the safety factor q also has a significant impact on the triangularity effects. In addition, the damping rate increases rapidly with triangularity when q is not too large and then saturates when δ is above about 0.3.

Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
analyses [8,[17][18][19][20]. It is well known that a simple GAM frequency performed in the ideal magnetohydrodynamic (MHD) model was first obtained by Winsor et al [17] tracing back to 1968, shown as , (1) in which q is the safety factor, c s is the sound speed, r is the flux label and R 0 is the major radius of tokamaks. The infinite aspect ratio and circular cross-section were employed in the derivation of equation (1) and were also adopted widely as basic assumptions in the later analytical literature [1]. However, it was found the GAM frequency weakly decreased with increasing inverse aspect ratio ϵ = r/R 0 in TEMPEST simulation [14]. More importantly, the GAM is more often observed at the edge region where the crosssection is far from circular due to the non-ignorable shaping effects of elongation κ and triangularity δ [1]. Already back in 2005, AUG [2] reported a significantly decreasing frequency of GAM with increasing κ, which was validated a year later in DIII-D [4] and reobserved in AUG [21,22] together with TCV [3]. Almost at the same time, Angelino et al [12] investigated the role of κ utilizing numerical simulations. Two years later, Angelino et al [15] also obtained the analytical results as with the Braginskii equation [23] and the Culham equilibrium [24]. For analytical work, a pioneering one was performed by Shi et al [25] based on a Solov'ev type equilibrium with MHD equations, including effects of ϵ, κ and δ, in principle giving ω 2 G ∝ (1 + κ −2 )/2. In 2008, Gao et al [18] obtained the frequency and damping rate of GAM from gyro-kinetic equations based on the widely used Miller equilibrium model [26], as ω 2 G ∝ 2/(κ 2 + 1). A series of papers [27][28][29] were further done for more shaping parameters and the comprehensive result was shown as [29] where v ti = 2T i /m i is the ion thermal velocity, τ = T e /T i is the temperature ratio between ions and electrons, s κ = r∂ r κ/κ stands for the elongation gradient with respect to the flux surface and ∆ ′ = ∂ r R 0 is the gradient of the Shafranov shift ∆. In addition, starting from the MHD model and a global shaping model different from Miller's, Wahlberg et al [30] obtained the dependence of frequency on not too large values of κ as ω 2 G ∝ (2 − κ − κs κ /4), consistent with equation (2) when |κ − 1| ≪ 1.
In short, different equilibria and models are employed, and the effects of κ, by far the most dominant shape parameter, are extensively analytically studied [15,18,25,[27][28][29][30] and qualitatively similar to experiments [2-4, 21, 22] and simulations [12][13][14][15][16]. Nevertheless, another critical shaping parameter, triangularity δ, has not been widely investigated experimentally and analytically [1]. Although most theoretical papers set δ = 0 directly, it has been proven to be related to the plasma turbulence especially for a negative δ [31]. In 2017, Sorokina et al [32] derived the GAM frequency in the MHD model with Miller equilibrium [26] including δ, κ and ϵ, Recently, TCV reported a nearly linear increasing GAM frequency with δ from − 0.2 to 0.3 measured by Z. Huang and S. Coda, and the dependence tended to be weakened for negative δ, as the figure 5(b) shown in [33]. Compared to TCV results, equation (3) indicates a fairly weak dependence on δ. By utilizing MHD equations and a global shaping model, another result containing the Shafranov shift gradient ∆ ′ was obtained by Wahlberg and Graves [33]. However, discrepancies still exist between the existing analytical results and TCV observations. More importantly, the damping mechanism and kinetic corrections are not taken into account in the framework of MHD.
In this paper, we present an exact solution to the dispersion relation of GAM in the presence of triangularity by revisiting the gyro-kinetic equation and employing Miller local equilibrium model [26]. The remaining part of the paper is organized as follows. In section 2, the Miller local equilibrium model with a non-circular cross-section, gyro-kinetic equations, governing equation of GAM and analysis with weak deformation are introduced. In section 3, the dispersion relation is obtained with an arbitrary κ and a finite δ; in the sufficiently large safety factor, the asymptotic solutions of GAM frequency and damping rate are derived. The triangularity effects on GAM are briefly discussed by the exact numerical solution and the asymptotic analytical solution in section 4. Finally, the summary is given in section 5.

Basic equations
We consider the widely used Miller local equilibrium model [26] with a flux surface (R, Z) written as where θ is the generalized poloidal angle. In Miller model, it is of basic difficulty to set priori values for s κ and s δ ≜ r∂ r δ. One common practice is adopting s κ = 0, s δ = 0, or employing empirical equations as s κ = (κ − 1)/κ and s δ = δ. The tokamak magnetic field ⃗ B = I(ψ)∇ξ + ∇ξ × ∇ψ can be described as [18] in which ξ is the toroidal angle, ψ is the magnetic flux, B 0 is the field at the magnetic axis R 0 (r 0 ), J = (∇r × ∇θ · ∇ξ) −1 is the Jacobian and dl/dθ = √ (dR/dθ) 2 + (dZ/dθ) 2 is the differential of the poloidal arc length with respect to the poloidal angle.
For simplicity, the electron response is ignored, and consequently the electrostatic potential remains constant on the flux surface, as ϕ =φ exp [i S(r) − iωt]. S(r) is adopted as S(r) = k r (r − r 0 ) and k r is the wavenumber in the flux coordinate system. It is acceptable to ignore the electron response since τ is not coupled with the shaping parameters (except s κ ) in equation (2) (also can be seen in equation (20) of [28]). We look forward to further improving this point in future work. The perturbed distribution of ions is determined by the gyrokinetic equation and can be solved as in which F i 0 is the equilibrium distribution of ions chosen to be the Maxwellian distribution here, E is the ion particle energy, J 0 is a zero order Bessel function, where is the parallel (perpendicular) ion velocity with respect to the magnetic field, and We stress that ρ i is the realistic length of the Lamor radius, which should be rewritten as ρ i |∇r| when transforming to the flux coordinate system, as shown in equations (6) and (7). In addition, without the electron response, the quasi-neutrality condition gives the governing equation of GAM, aŝ

Analysis with weak deformation
Equations (4)-(8) are a closed set of equations to describe GAM for arbitrary values of shaping parameters in Miller local equilibrium model [26]. To obtain analytical results, we adopt large aspect ratio and high safety factor (typically, q ≳ 3), and assume δ, , which are valid for most realistic tokamak plasmas. By retaining terms of O(ϵ 2 ), equation (7) can be rewritten as with Since the shaping effects introduce the complex dependence of ω t , ω d and k r |∇r|ρ i on the poloidal angle θ, it is a lengthy and tedious process to utilize the general solution equation (10) directly. Here we consider the case of k r ρ i ≪ 1, that is, neglecting the effects of high order finite-Larmor-radius and finite-orbit-width [35,36]. Through inserting equation (10) into equation (6) and by keeping the second-order terms of k r ρ i , then, the perturbed ion distribution equation (6) splits into , It is noticed that F I 1 corresponds to the term of O(k r ρ i ) and has no contribution after the magnetic surface average integral, which is consistent with the case of a circular cross-section. In addition, the correction of F IV 1 to the dispersion relation is lower than F II 1 or F III 1 by the order of O(ϵ/q 2 ), and consequently it is also not required to be included. For the terms of exp(k r f(θ)) and exp(ϵf(θ)) in equation (11), calculations can be simplified by using Taylor expansion directly instead of the Bessel function.

Dispersion relation
After inserting equation (11) into equation (8) and tedious algebraic calculations, the dispersion relation is obtained as in which we denote Here, Z(ζ) is the plasma dispersion function, and G n , D are relatively complex resulting from shaping effects, provided as It should be emphasized that the trapped particle effect is not considered here, which requires q ≫ 1 or a sufficiently large aspect ratio [37]. The original dispersion relation equation (12) is so involved that no intuitive analytical results can be obtained directly, although exact numerical solutions can be calculated. Fortunately, an asymptotic solution can be established with ζ ≫ 1, which requires a sufficiently large safety factor. Actually, previous work [37] shows sufficient agreement between asymptotic solution and numerical solution with a not too large safety factor (typically, q ≳ 2) in the case of a circular cross-section. Regardless, the plasma dispersion function can always be asymptotically expanded as with large enough q. Neglecting terms of order higher than O(ζ −4 ), the asymptotic analytical solutions of GAM frequency and damping rate are obtained, in which we denote When circular shaping parameters are specified (κ = 1, s κ = δ = s δ = ϵ = ∆ ′ = 0) in equations (13) and (14), previous results of GAM frequency and damping rate are recovered [37][38][39][40]. Compared with GAM frequency of Gao's result [29], the main differences are the absence of any τ dependence in equation (13) and the absence of any δ dependence in equation (2). The main contribution of τ may be in the term of 7 4 + τ as shown in equation (2). Therefore, when the dispersion relation derived here is to be utilized for quantitative comparison with experiment, ω (δ) /ω (δ = 0) is a better choice rather than ω (δ) to eliminate the discrepancy brought by τ as much as possible. The slight discrepancy of terms O(s κ ) between equations (2) and (13) is due to diverse selections of magnetic surface average integral in the quasi-neutrality equation, and the distinction disappears in comparison with Gao's another frequency result, equation (20) of [28]. For the same reason, when τ = 0 and s κ = 0, equation (14) coincides with Gao's damping rate, i.e. equation (13) in [27] (mistakes in writing: q 5 is missing; ω GAM → qR 0 ω GAM /v ti in square brackets). In addition, we stress that the assumption of (κ 2 − 1)/(κ 2 + 1) ∼ O(ϵ) is not adopted in this paper, and consequently, equations (13) and (14) are suitable for an arbitrary value of κ. Except for the inherent discrepancies caused by gyro-kinetic model and MHD model, equation (13) coincides with equation (3) in the shaping effects, the result of Sorokina et al [32], and is also qualitatively similar to the result of Wahlberg et al, equation (1.4) of [33].

Discussion
It should be pointed out in advance that ∆ ′ is associated complexly with other shaping parameters, global magnetic shear as well as pressure gradient [41]. Moreover, κ, δ and their derivatives also influence each other [42]. Considering these, utilizing a set of self-consistent shaping parameters to investigate δ effects is beyond the scope of this paper. Consequently, we have to follow the previous analytical work [26,33,[43][44][45][46][47][48] and treat shaping parameters including ∆ ′ as independent. However, when equations (12)- (14) are used to be compared with experiments or simulations in the future, the self-consistent shaping parameters according to experimental measurements or simulation calculations are desired to be adopted in equations (12)- (14). In short, in this section, we are limited to the case of shaping parameters independent with each other to briefly discuss triangularity effects on GAM. Now, we first focus on the triangularity effects on the frequency of GAM by utilizing equations (12) and (13). The asymptotic analytical result equation (13) clearly indicates that the nearly linear increasing frequency with triangularity is due to its coupling with the finite aspect ratio ϵ and Shafranov shift gradient ∆ ′ , while the increasing growth slows down for negative δ due to the existence of the term O(δ 2 ). This dependence indicated by equation (13) is consistent with the observation in TCV on the trend [33]. If the triangularity gradient s δ is set to s δ = δ rather than zero, it also slightly influences the triangularity effects when δ is not near zero. The frequency dependence on δ can also be numerically calculated from the original dispersion relation equation (12), and a similar trend to the analytical one is shown as figure 1. However, figure 1 indicates that the frequency dependences on triangularity in the exact numerical solutions are much stronger than those in the asymptotic analytical solutions. This discrepancies between them mainly result from the insufficient safety factor (q = 3.57 in figure 1). Although there is no evidence in equation (13) that the safety factor can influence the effect of triangularity on GAM, the direct numerical calculation of equation (12) indicates that a not too large safety factor can significantly enhance the triangularity effect, as shown in figure 2. This enhancement effect from safety factor is not predicted in equation (13) since it is only valid for a sufficiently large safety factor (q ≫ 3), which is consistent with the fact that the enhancement in the exact numerical solution almost disappears for q ⩾ 5, as shown in figure 2. Consequently, there exists a special safety factor as a so-called 'sweet spot' for experimental observations of triangularity effects on GAM frequency. The special q value depends on other shaping parameters and is about 3.5 for parameters used in figure 2. In addition, it should be pointed out that the discrepancy between analytical and numerical results still exists in the large q limit, about 1.3% for ω(δ = 0.3)/ω(δ = 0) at q = 10, as shown in figure 2. This discrepancy actually results from terms of order O(ϵ 3 ) neglected in the transition process from equations (12) to (13). Moreover, equation (3), the MHD result of Sorokina et al, is plotted in figure 1, and it is mainly the Shafranov shift  (12). The analytical curves are according to the asymptotic analytical solution equation (13). The TCV experimental data (blue triangle point) are the same as figure 5(b) in [33] (courtesy of Z. Huang and S. Coda). In the numerical and analytical curves, on the basis of best guess measurements of TCV result, q = 3.57, κ = 1.3, ϵ = 0.2, sκ = 0, ∆ ′ = −0.35 are adopted according to [33] in the case where TCV parameters uncertainty exists.   For triangularity effect on the GAM damping rate, additional harmonic transit resonances are induced by nθ in ω t , ω d and k r |∇r|ρ i due to shaping effects and thus introduce terms exp(−ζ 2 /n 2 ) in the damping rate equation (14), which increase quickly with n. In addition, ω itself is a function of triangularity as equation (13) and damping rate γ ∝ exp(−ζ 2 /n 2 ) is sensitive to the value of GAM frequency. Therefore, a better choice is direct numerical calculation in the original dispersion relation equation (12) especially for a small q, as shown in figures 4 and 5. In figure 4, it is found that the increasing triangularity significantly increases the damping rate when q is not too large, and the numerical damping rate saturates as δ increases to a certain value (about 0.3 in figure 4). In other words, GAM is more easily damped for a  (12) and (14), respectively. positive triangularity and more resistant for a negative triangularity. Figure 5 not only reveals the discrepancy between the analytical and numerical results, but also indicates that the damping rate rapidly decreases to zero with the increase of q, which is consistent with that under a circular cross-section.

Summary
In the present work, the shaping effects on GAM, in particular the triangularity effects, are investigated by revisiting gyro-kinetic equations and employing Miller local equilibrium [26] for tokamaks with the non-circular cross-section. The dispersion relation of GAM is obtained as equation (12), and the asymptotic analytical results, equations (13) and (14), are derived for a sufficiently large safety factor (typically, q ⩾ 5). Both the asymptotic analytical and the exact numerical results show a nearly linear increasing frequency with triangularity δ and a relatively weakened dependence for negative δ, which are consistent with TCV observations on the trend [33] and the previous analytical result equation (3) from ideal MHD model, as shown in figure 1. Although the analytical solution indicates that the strength of triangularity effect on the GAM frequency is only dependent on the inverse aspect ratio ϵ and Shafranov shift gradient ∆ ′ , numerical calculations show that it is also significantly influenced by the safety factor, as shown in figure 2. For rough shaping parameters of TCV as used in figure 2, the increase in frequency with increasing triangularity reaches its maximum for safety factor values of about 3.5, which means it is better in a not too large q when observing triangularity effects experimentally. In figure 4, the damping rate of GAM is found to increase quickly with the triangularity for a not too large safety factor and then saturate above a certain value of δ (about 0.3 in figure 4). The GAM frequency and damping rate presented analytically and numerically in this paper can provide a convenient and instructive estimation to the investigation of shaping effects, especially the triangularity effects on GAM applied to relative tokamak experiments.