Impact of shaping on microstability in high-performance tokamak plasmas

We have used the local-δf gyrokinetic code GS2 to perform studies of the effect of flux-surface shaping on two highly-shaped, low- and high-β JT-60SA-relevant equilibria, including a successful benchmark with the GKV code. We find that for a high-performance plasma, i.e. one with high plasma beta and steep pressure gradients, the turbulent outwards radial fluxes may be reduced by minimizing the elongation. We explain the results as a competition between the local magnetic shear and finite-Larmor-radius (FLR) stabilization. Electromagnetic studies indicate that kinetic ballooning modes are stabilized by increased shaping due to an increased sensitivity to FLR effects, relative to the ion-temperature-gradient instability. Nevertheless, at high enough β, increased elongation degrades the local magnetic shear stabilization that enables access to the region of ballooning second-stability.


Introduction
One of the performance-limiting factors on tokamak experiments is the significant turbulence-dominated radial transport of heat, momentum and particles from the core plasma out to the edge. One method of affecting this turbulent transport is to change the shape of the axisymmetric fluxsurfaces traced out by the confining magnetic field; this has been validated both experimentally and in numerical simulations [1][2][3][4][5][6][7][8]. Flux-surface shape is often characterised by the elongation κ and triangularity δ, where κ(r) ≡ Z(r, θ max )/r * Author to whom any correspondence should be addressed. 5 Present address: Tokamak Energy Ltd, Abingdon, Oxfordshire, OX14 4SD, United Kingdom of Great Britain and Northern Ireland Original content from this work may be used under the terms of the Creative Commons Attribution 3.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. and δ(r) ≡ [R 0 (r) − R(r, θ max )]/r, where R(r, θ) is the plasma major radius, θ is the poloidal angle, Z(r, θ) is the vertical distance above the midplane, r ≡ [R(r, 0) − R(r, π)]/2 is the plasma half-diameter at the flux-surface midplane, R 0 (r) ≡ [R(r, 0) + R(r, π)]/2 is the major radius of the flux-surface center and θ max is the poloidal angle of the maximum Z. A labelled example flux-surface shape is shown in figure 1. The above expressions are valid for the up-down symmetric plasmas studied in this work, although similar expressions can be used to account for up-down asymmetry. Plasma shaping is known to affect the stability of magnetohydrodynamic (MHD) modes. Specifically, increased elongation and triangularity increase the threshold 'Troyon' plasma β above which dangerous (e.g. external kink-ballooning) MHD instabilities are triggered: β Troyon = β N I p /(aB ref ), where β is the ratio of thermal to magnetic pressure, β N 2% is an empirically calculated scaling factor, a is the minor radius of the last closed flux-surface (LCFS) and B ref is the reference magnetic field, defined as the on-axis toroidal field [9]. The Troyon beta limit increases linearly with plasma current I p which, for fixed safety factor q, is increased by elongation and triangularity. Shaping can however have other important consequences, particularly for vertical displacement events (VDE) which are associated with disruptions. Whilst increased elongation can ameliorate the aforementioned MHD instabilities, more severe controls for VDE are then required. This necessitates the careful investigation of shaping effects from both MHD/disruption and transport perspectives.
The impact of shaping on transport due to microinstabilities has been less-thoroughly studied, and those studies that have been performed have focussed only on a relatively small region of the vast multi-dimensional shaping parameter space. Despite this, increased elongation is generally thought to be stabilizing [1][2][3][4], whereas triangularity can have varying effects. Both the TCV and DIII-D experiments have reported that negative triangularity yields improved performance, allowing performance comparable to H-mode whilst maintaining L-mode edge profiles [5][6][7]. Meanwhile, [2] has reported a dependence on elongation of the effect of triangularity, being stabilizing at large κ and somewhat destabilizing at moderate κ. There are various explanations for these observed effects on the turbulence; some suggest that flux-surface shaping modifies locally the driving gradients [8], whilst others point to its effect on local magnetic shear [10].
In this work we extend an earlier modelling study by Nakata [11] studying two JT-60SA-relevant magnetic equilibria [12]. Following a benchmark between the local, δ f -gyrokinetic codes GKV [13,14] and GS2, we present a study on the effects of triangularity and elongation in these two equilibria, both with and without electromagnetic effects. In section 2, we present how shaping influences turbulence in the gyrokinetic framework, as a combination of the local-magnetic-shear stabilization of the toroidal drive, as well as finite-Larmor-radius (FLR) stabilization. In section 3, we present the results of the benchmark, followed by electrostatic linear and nonlinear shaping scans in {δ, κ}. We use our model to explain the qualitative trends observed in these scans, including the novel result that increased elongation can be destabilizing. Finally, we present linear electromagnetic shaping scans, showing how kinetic-ballooning-modes (KBMs) react differently to fluxsurface shape. The work is then summarized in section 4.

Model
To understand how plasma shaping affects the local, δ f gyrokinetic framework in which GS2 works, we write the Fourier-analysed (perpendicular to the equilibrium magnetic field) collisionless gyrokinetic equation and field equations as: where h s is the non-adiabatic part of the fluctuating distribution function for a species labelled s (note that we have omitted the wavenumber index on fluctuating quantities to simplify notation), v is the particle velocity, J n is the n th -order Bessel function of the first kind, Z s e is the charge, T s is the temperature, F 0,s is the assumed-Maxwellian equilibrium distribution function, is the gyro-averaged gyrokinetic potential, δφ is the fluctuating electrostatic potential, δA and δB are the parallel components of the magnetic vector potential and field, respectively, {. . .} indicates the Fourier-transformed Poisson bracket, contains the non-E × B magnetic drifts, indicates an rderivative, p is the pressure, Ω s is the Larmor frequency, V χ ≡ i(cχ/B 2 )B × k ⊥ is the drift velocity due to the interaction between the equilibrium magnetic field and the fluctuating electric and magnetic potentials, u s ≡ |k ⊥ ||v ⊥ |/Ω s and k ⊥ ≡ k ψ ∇ψ + kα∇α where ψ is the poloidal magnetic flux, used as a radial coordinate, andα is the field-line label. One way that the flux-surface shape enters these equations is via ∇ψ and ∇α in k ⊥ . We note that the magnetic field strength and parallel streaming terms are also affected by the flux surface shape, but this effect is often small in the inverse aspect ratio ≡ r/R 0 . The perpendicular wavevector k ⊥ can strengthen (weaken) the toroidal drive term if (mis)aligned The effect of magnetic shear on a ballooning-type perturbation. The instability is driven strongly at the outboard midplane where V D,s · k ⊥ is maximal. The perturbations are advected along the field lines to a different poloidal (and toroidal) position. With moderate positive magnetic shear, the outboard side lags behind the inboard side poloidally, so the perturbation retains some major-radial extent and remains strongly driven. For negative magnetic shear, the inboard side lags and the perturbation gains more vertical extent, thus coupling more weakly to the predominantly vertical magnetic drifts. A similar argument is used to explain why large positive values of magnetic shear have a stabilizing effect.
with V D,s ; this is closely linked to local-shear stabilization as shown in figure 2. It can also provide FLR stabilization via the gyro-average, manifested as the Bessel functions. To calculate ∇ψ and ∇α, we work in toroidal {r, θ, ζ} and cylindrical {R, ζ, Z} coordinates, where r, θ, Z and R retain their previous definitions and ζ is the toroidal angle. By using the Clebsch (B = ∇α × ∇ψ) and axisymmetric (B = I∇ζ + ∇ζ × ∇ψ) representations of the equilibrium magnetic field, one can write an expression forα: whereq(r, θ) ≡ B · ∇ζ/B · ∇θ is the local safety factor. It follows that ∇α contains the local magnetic shears ≡ r(logq) and can be written as As shown in figure 2, the local magnetic shear affects the toroidal drive by shearing turbulent structures towards or away from the major radial direction in which they are most strongly driven [15,16]. To obtain expressions for the local magnetic shear as well as the other geometric coefficients appearing in the gyrokinetic system of equations, one can specify R(r, θ) and Z(r, θ) for a given flux-surface. A typical parametrization, and indeed the one we use here, is the Miller parametrization [17] which uses nine parameters to specify the poloidal cross section: where δ is the triangularity and κ is the elongation. Only first derivatives of r-dependent quantities are specified, and the Grad-Shafranov equation [18,19] is used to ensure that the equilibrium locally satisfies MHD force balance. Consequently, β ≡ β(log p) must also be specified. More details on how to obtain expressions fors and |k ⊥ | 2 are provided in appendix A.

Results
Following the work of Nakata [11], we present the results of a benchmark between GS2 and the analogous code GKV, in figure 3. The benchmark was performed at three radial (ρ, the normalized toroidal magnetic flux) positions with adiabatic electrons and a single ion species, using Miller parametrizations of two numerical equilibria: one low-β(1.5%) and one high-β(3.7%) [20]. The flux surface shapes and equilibrium parameters are given in figure 4 and table 1. There is reasonable agreement between the linear growth rates γ and real frequencies ω calculated by the two codes. There is a small discrepancy that was also observed in a previous benchmark between the two codes [21]-this may be due to algorithmic differences that affect the numerical dissipation. Despite this, the agreement was sufficient that we extended the study beyond the benchmark for both equilibria by adding kinetic electrons and then electromagnetic perturbations; these results are shown in figure 5. The addition of kinetic electrons approximately doubles all linear growth rates and destabilizes trapped electron modes for k y ρ ref 1, where the normalized poloidal wavenumber k y is related to kα via k y = kαB ref /ψ . Adding full δB perturbations destabilizes KBMs for both equilibria, though they are more potent in the high-β case, being the dominant mode for k y ρ ref > 1 at all the simulated radii. These modes are identified as KBMs due a step change in the real frequency. Studies were also performed including either δB or δA . With only δB compressive magnetic fluctuations, the results were almost identical to the electrostatic simulations. With only field line bending (δA ), a KBM-like mode appeared at low k y , but with a growth rate approximately 50% smaller than the full δB fluctuations. In both equilibria, the electromagnetic effects stabilize ITG.

Electrostatic, linear shaping scans
To guide us towards the flux-surface shape that minimizes turbulent heat flux, electrostatic scans were performed in the triangularity and elongation Miller parameters for the ρ = 0.5 equilibria. Simultaneously scaling δ and κ had no impact on the qualitative trends observed; the effect on the maximum linear growth rates was less than 10% for the majority of the {δ, κ} space, but reached ∼25% at maximum triangularity. It was therefore deemed unnecessary to simultaneously scale δ and κ in any of the simulations. These electrostatic simulations are equivalent to β 0 electromagnetic simulations. We choose to fix β ≡ β(log p) , which corresponds to modifying (log p) to compensate for the 'reduction' in β. As such, until electromagnetic effects are included, we refer to the equilibria as low-and high-(log p) . The results of the scan are shown in figure 6(a)). For the low-(log p) equilibrium, increasing the elongation has an almost universal stabilizing effect, whereas the triangularity is destabilizing at low elongation and stabilizing at high elongation. The result is that maximal shaping minimizes the linear ITG instability. The maximum growth  rates for the high-(log p) equilibrium are, for all {δ, κ}, less than for the low-(log p) equilibrium. Additionally, a moderate increase in elongation is destabilizing even up to δ ∼ 0.2, with the result that circular flux surfaces are as stable as those with maximal shaping. Using figures 6(b) and (c) we show that the trends in figure 6(a) can be explained via a combination of local magnetic shear and FLR stabilization. To quantify the effect of local magnetic shear, we use its value at the outboard midplane where the ballooning modes are driven most strongly. Meanwhile, we quantify the FLR stabilization via the integral π/2 −π/2 |k ⊥ | 2 dθ, which gives a measure of how much the mode is restricted over the bad-curvature region which extends approximately from −π/2 θ π/2. For κ 1, the local magnetic Table 1. Equilibrium parameters for the two equilibria. ≡ r/R 0 is the inverse-aspect ratio,ŝ is the magnetic shear, Δ ≡ R 0 is the Shafranov shift and α ≡ −β(log p) q 2 R 0 where q is the safety factor. T e = T i for both equilibria. shear is rapidly made less negative by increasing triangularity. In contrast, |k ⊥ | changes relatively slowly with triangularity, so the overall effect is a strong destabilization. At larger κ 2, the local shear is insensitive to changes in triangularity, so the increased |k ⊥ | 2 gives a net stabilization with δ. The effect of elongation ons changes with triangularity and α. At large δ, increased elongation makes the local shear more negative whilst also increasing |k ⊥ | 2 ; these two effects work in tandem to provide strong stabilization. In figure 7 we show that for δ 0.2, increased elongation reduces the local shear for small α 0.5, but increases it for larger α 0.5. Both equilibria have sufficient α for increased elongation to increasẽ s, although the high-(log p) equilibria has higher α and sõ s is sufficiently sensitive to changes in κ that it outcompetes the increased FLR stabilization afforded by higher κ. In contrast, the low-(log p) equilibrium, with lower α, is dominated by the stabilizing effect of k 2 ⊥ -stabilization. In appendix A, a simplified analytical model is developed; this gives qualitatively identical results to what is observed in figure 7. We note that quasilinear heat-flux estimates were also calculated. These predicted the low-(log p) equilibrium to have a lower heat flux than the high-(log p) case, but this was not borne out in nonlinear simulations.

Electrostatic, nonlinear shaping scans
We also present nonlinear electrostatic scans in elongation for both equilibria in figure 8. All simulations were carried out with 32 parallel grid points, 22 k y values, 256 k x values, with Δk x Δk y = 0.047. The polar velocity space grid contained 16 radial energy grid points, 20 passing pitch angles and up to 33 trapped pitch angles [22,23]. The trends observed in heat flux are similar to those of the linear growth rates: at δ = 0, increased elongation monotonically stabilizes the low-(log p) equilibrium whilst the high-(log p) case is destabilized up to κ 1.6. We note that the different κ-position of the maximum (cf κ 1.2 in the linear results) is the main discrepancy between the linear and nonlinear results. At δ = 0.5, increased elongation monotonically stabilizes the high-(log p) equilibrium. These results again suggest that maximal shaping minimizes turbulent transport for a low-β plasma with moderate Figure 6. Scans in triangularity and elongation for two equilibria, showing (a) linear growth rates, (b) local magnetic shear at the outboard midplane, (c) integral of |k ⊥ | 2 over −π/2 θ π/2. The black crosses indicate the nominal shapes. We note that the differences in local shear and |k ⊥ | 2 between the two equilibria are small. This is because whilst the low-(log p) equilibrium has a normalized pressure gradient almost three times smaller than the high-(log p) one, the increased q and R 0 result in a comparatively similar value of α, which governs the overall effect of pressure gradient on the geometrical coefficients. The colours have been chosen such that blue indicates increased stability.
pressure gradient. However, in a low-β plasma with steep pressure gradient, similar transport levels could be achieved with circular flux-surfaces as with maximal shaping. The fraction of energy in the zonal flow is also plotted in figure 8(b). To calculate this, we take the ratio of the following proxies for the zonal flow energy: , and the total energy contained in both turbulence and zonal flows: where Γ k x ,k y (θ) denotes the gamma function evaluated at k ⊥ (k x , k y , θ)ρ ref and . . . θ indicates an average over poloidal angle [24]. For the low-(log p) equilibrium at δ = 0 and the high-(log p) at δ = 0.5 the zonal energy fraction increases monotonically with elongation, reflecting the decreased transport. Similarly, for the high-(log p) , δ = 0 simulation the minimum zonal energy fraction coincides with the maximum heat flux. These results suggest that the effect of elongation on turbulent transport is closely linked to its effect on the relative strength of the zonal flow.

Electromagnetic shaping scans
In figure 5 electromagnetic effects are observed to have a significant effect on the linear growth-rate spectra. We therefore include electromagnetic effects in the shaping studies. Figure 9 shows the dependence of the maximum linear growth rates on triangularity and elongation for both equilibria. In the high-β equilibrium, a KBM is excited for all {δ, κ}, and the growth rates are everywhere higher than the electrostatic case. Furthermore, the effect of plasma shaping becomes uniformly stabilizing. The nominal low-β equilibrium is on the threshold of KBM stability, as indicated by a step change in the wavenumber and real frequency of the fastest growing mode either side of the green line in figure 9(a). In contrast to the electrostatic results, the low-β equilibrium has everywhere lower linear growth rates than the high-β case, and triangularity is no longer destabilizing at low elongation.
These results can be explained by an increased sensitivity of the KBM to FLR stabilization, compared to the electrostatic ITG. The local magnetic shear is most sensitive to triangularity , showing (a) linear growth rates maximised over k y , (b) local magnetic shear at the outboard midplane, (c) integral of |k ⊥ | 2 over −π/2 θ π/2. The green and black crosses indicate the nominal parameters for the low-and high-β equilibria, respectively. The scan in α was performed at fixed q and R 0 , so it corresponds to a scan in (log p) . The nominal {α, κ} value for the low-(log p) equilibrium is also shown by a green cross; a similar scan with the low-(log p) equilibrium parameters gives qualitatively similar results. at κ = 1, and only here is it able to compete with the |k ⊥ | 2 stabilization, giving the almost constant linear growth rates at κ = 1. For all other elongations, the KBM is not sufficiently sensitive tos and is thus stabilized by |k ⊥ | 2 which increases with both δ and κ. As mentioned, the KBM is the dominant mode for all {δ, κ} at high β, but for the low-β case it becomes sub-dominant for strong shaping.
In figure 10 we show scans in β at fixed (log p) . As β increases,s(θ = 0) is reduced via the corresponding increase in α ∝ β. Two such scans are shown, at nominal and 1.5 times nominal (log p) . At low β, we observe the β-stabilization of ITG modes, and note that the KBM threshold β increases with elongation due to FLR stabilization, consistent with figure 9. We also observe the persistence of KBMs (albeit with significantly reduced linear growth rates) into the region of second MHD stability, as previously reported in [25], for sufficiently high β approaching 10%. The reduction in growth rates in this region can be attributed to the increasing localshear-stabilization due to larger α ∝ β. As κ is increased in the region of second stability, the local-shear stabilization is reduced rapidly enough that it outcompetes FLR stabilization. This has the effect of removing the second stability afforded by increased β.
To summarize, the differences in the results obtained from the nominal-and increased-(log p) scans are as follows: the growth rates are everywhere lower in the increased-(log p) case, and the KBM reacts more strongly to increased (log p) compared to the ITG. This is because the absolute change in  α scales linearly with β, so the local magnetic shear is less reduced at low-β. In the increased-(log p) case, the second stability region is accessed at lower β. This in turn lowers the threshold β at which κ becomes destabilizing. This suggests that to minimize turbulent transport in a high-performance plasma with steep pressure gradients and high β, the flux surface elongation should be minimized.

Summary
In this work we have performed a successful benchmark between the local δ f -gyrokinetic codes GS2 and GKV for two equilibria with significant shaping effects. We have also performed, with the inclusion of electrons as a kinetic species, linear scans in the flux-surface shaping parameters δ and κ to study the effect on turbulent transport. We find that for low-β equilibria, triangularity is stabilizing at high elongation, and vice versa. We also find the novel result that elongations of κ ∼ 1.5 can be destabilizing in regions of steep normalized pressure gradient. This trend in linear growth rate is also observed in the nonlinear heat-flux. We explain these results as a competition between the effects of shaping on the local magnetic shear and on FLR stabilization. Electromagnetic scans were also performed, and KBMs were destabilized for both equilibria. Our results show that for the nominal equilibrium parameters, the KBM is stabilized monotonically by increased shaping; we infer from this that for the KBM, FLR stabilization is stronger relative to the local magnetic shear. Scans in {β, κ} show that at sufficiently high β, the local magnetic shear becomes sensitive enough to κ that it can compete with FLR stabilization and compromise second-stability. The threshold β at which this occurs decreases with increased normalized pressure gradient. Therefore, we suggest that for a high-performance plasma, i.e. one with high-β and steep pressure gradients, the turbulent outwards radial fluxes may be minimized by reducing elongation as much as possible. Conversely, for a moderate β plasma, we would suggest maximal elongation to FLR-stabilize the KBMs.

Appendix A. Calculation of the local magnetic shear and perpendicular wavenumber
In this section we flesh out the details of how to determine the effect of shaping parameters on stability. As an example, we show an analytic calculation of the effect of elongation to leading order in 1. As discussed in section 2, this involves determining the local magnetic shear and the perpendicular wavenumber.
The local safety factorq is the ratio of the toroidal to poloidal component of the magnetic field: where we used the axisymmetric form of the magnetic field: and defined the Jacobian J X for the transformation between {R, ζ, Z} → {X, θ, ζ} where X is any flux-surface label. This also allows us to define ψ in terms of the safety factor q: Therefore, the local magnetic shear is where J ≡ J r − J r ψ /ψ . To evaluate this, we use the Grad-Shafranov equation: By expanding the divergence and substituting explicit forms for ∇r and ∇θ using a prescribed R(r, θ) and Z(r, θ), we arrive at the following expression for J/J r : where superscript θ indicates a θ-derivative. To get an expression for I /I, we take the r-derivative of equation (13) and move all terms under the integral: This can be used in conjunction with equation (16) integrated over θ to give an expression for I : With these expressions, one can generate expressions fors and thus |k ⊥ | 2 .

A.1. Analytical expressions for concentric elliptical flux-surfaces
We next use a simplified Miller parametrization that includes only elongation to determine the effect of elongating a circular plasma: We choose to order κ ∼ κ/R 0 -this means that κ does not enter to leading order in small inverse aspect ratio (  1): We proceed in this limit to determines and |k ⊥ | 2 , beginning with the following quantities: |∇r| 2 ≡ R 2 J 2 r R θ 2 + Z θ 2 = 1 + (κ 2 − 1)cos 2 θ κ 2 (23) To determine I /I, we calculate each of the three integrals that appear in equation (18) Figure 11. Profiles calculated from our simple analytical model that includes just elongation and pressure-gradient effects in the small-inverse-aspect-ratio limit. We show (a) the local magnetic shear at the outboard midplane and (b) the integral of |k ⊥ | 2 over −π/2 θ π/2, as a function of α and κ.
Combining these, we find that where we will find that the leading order piece cancels with a term in J/J r , so we retain the O(1) terms. Using this, we find: κŝ − α cos θ 1 + (κ 2 − 1)cos 2 θ .
Comparing the size of the terms that comprises, the largest is J/J r . Then, to leading order, we find s = κŝ − α cos θ 1 + (κ 2 − 1)cos 2 θ .
Now we can uses to find ∇α and thus |k ⊥ | 2 : where and ϑ(θ, κ) ≡ arctan tan θ κ , which gives: where θ 0 ≡ k x /(ŝk y ). These expressions reduce to the typicalŝ − α results when κ = 1, since Λ(θ, 1) = sin θ and ϑ(θ, 1) = θ. In figure 11 we plot for this model the local shear at the outboard midplane, and the integral of |k ⊥ | 2 over the outboard side, as a function of α and κ. Comparing figures 11(a) and (b) to figures 7(b) and (c), we see good agreement despite ignoring the effects of triangularity, Shafranov shift and finite-aspect-ratio.