Critical gradient turbulence optimization toward a compact stellarator reactor concept

Integrating turbulence into stellarator optimization is shown by targeting the onset for the ion-temperature-gradient mode, highlighting effects of parallel connection length, local magnetic shear, and flux surface expansion. The result is a compact quasihelically symmetric stellarator configuration, admitting a set of uncomplicated coils, with significantly reduced turbulent heat fluxes compared to a known stellarator. The new configuration combines low values of neoclassical transport, good alpha particle confinement, and Mercier stability at a plasma beta of almost 2$\%$.

Introduction.-Aprimary obstacle for the success of magnetic confinement fusion is the transport caused by instabilities such as the ion temperature gradient (ITG) mode, which is thought to significantly reduce plasma confinement in experiments such as the Wendelstein 7-X stellarator [1][2][3].To overcome the losses from such turbulence, a given configuration can be scaled up in size and heating power.A less costly alternative, currently explored, is to shape the magnetic field to alleviate the turbulence.This option could be particularly appealing for reactor scenarios, in which it will likely be difficult to achieve density gradient stabilization of turbulence via pellet injections [4,5], since the penetration distance of pellets may be limited in comparison to the minor radius of a reactor.
To achieve turbulence optimization via shaping, several strategies have been developed to reduce the rate that turbulent transport increases ("stiffness") as a function of the ion temperature gradient [6][7][8][9][10].Another approach is to target the onset ("critical") gradient of significant turbulent transport [11][12][13], which relies primarily on linear physics of ITG modes themselves [14][15][16][17] and avoids the hard problem of solving turbulence in the full range of toroidal geometries.Here we demonstrate optimization using a critical gradient (CG) approach, which even leads to reduced stiffness of ITG turbulence in the nonlinear regime [12], albeit with some implied trade-offs for integrated stellarator optimization.
In this Letter, we first show CG optimization targeting the absolute threshold for ITG modes, producing the largest critical gradient of all stellarators known to us, while sacrificing magneto-hydrodynamic (MHD) stability.We then show that without compromising MHD stability, or other key properties needed by a stellarator design, one may target the CG of only the toroidal branch of the ITG mode, based on the assumption that turbulence intensity will be small below this threshold.This model highlights the familiar stabilizing effects of local shear, but gives greater emphasis to short connection lengths between regions of "good" and "bad" magnetic curvature.The resulting objective function is used * gar@ipp.mpg.devia optimization to produce a quasi-helically-symmetric configuration with strongly reduced ITG turbulence compared to a known stellarator experiment (HSX) [18], in addition to acceptable levels of neoclassical losses, alpha particle confinement, coil complexity, and MHD stability, thus completing the picture for an initial stellarator concept with improved ion confinement. Definitions.
-Following [19], we use the standard gyrokinetic system of equations [20] to describe electrostatic fluctuations destabilized along a thin flux tube tracing a magnetic field line.The ballooning transform [21,22] is used to separate out the fast perpendicular (to the magnetic field) scale from the slow parallel scale.The magnetic field representation in field following (Clebsch) representation reads, B = ∇ψ × ∇α, where ψ is a flux surface label and α labels the magnetic field line on the surface.The perpendicular wave vector is then expressed as k ⊥ = k α ∇α + k ψ ∇ψ, where k α and k ψ are constants, so the variation of k ⊥ (ℓ) stems from that of the geometric quantities ∇α and ∇ψ, with l the field-line-following (arc length) coordinate.
We assume Boltzmann-distributed (adiabatic) electrons, thus solving for the perturbed ion distribution g i (v ∥ , v ⊥ , l, t), defined to be the non-adiabatic part of δf i (δf i = f i − f i0 ) with f i the ion distribution function and f i0 a Maxwellian.The electrostatic potential is ϕ(l), and v ∥ and v ⊥ are the particle velocities parallel and perpendicular to the magnetic field, respectively.The gyrokinetic equation reads where ω is the mode frequency, is the diamagnetic frequency, and J 0 = J 0 (k ⊥ (l)v ⊥ /Ω(l)) is the Bessel function of zeroth order.The thermal velocity is v T = 2T /m, the thermal ion Larmor radius is , n and T are the background ion density and temperature, q is the ion charge, φ = qϕ/T is the normalized electrostatic potential, and Ω = qB/m is the cyclotron frequency, with B = |B|.The magnetic drift frequency in the low β approximation is For simplicity in this analysis, we set k ψ = 0. We rewrite the drift frequency as ω d (ℓ) ∝ K d (ℓ) ≡ a 2 ∇α • b × κ, referring to K d (ℓ) as the "drift curvature" and to individual regions of bad curvature along the field line (where K d > 0) as "drift wells".We define a radial coordinate r = a ψ/ψ edge , with a the minor radius corresponding to the flux surface at the edge, and ψ edge the toroidal flux at that location.The temperature gradient scale length is measured relative to the minor radius, a/L T = −(a/T )dT /dr.To study the most unstable ITG mode conditions, we have neglected certain stabilizing factors such as the density gradient [23,24] and plasma beta (electromagnetic effects) [25].
Finally, the gyrokinetic system is completed by quasineutrality, with τ = |q e |T /(qT e ), T e the electron temperature, and q e the electron charge.
Thresholds for ITG modes.As argued in [12], the CG can be estimated by using the model where R eff is the local effective radius of curvature determined by the profile of K d (ℓ), and the effective parallel connection length for modes near the absolute threshold, L ∥ Floquet [see discussion above eqn.( 4) in [12]], is determined by the relative size of good curvature outside the drift well, which may stabilize extended Floquet-like modes.The effect we seek to enhance is contained in the first term proportional to a/R eff and thus to the size of "bad" curvature on the outboard midplane.
It is expected, however, [13,26,27] that the onset of toroidal ITG modes, as can be inferred from linear spectra in gyrokinetic simulations [17], should lead to noticeable increases in nonlinear heat fluxes at a second, larger CG.The turbulence found below this onset (in the Floquet-like or slab-like regime) is thought to be more benign.We therefore also focus on the toroidal ITG mode [14,17,19,28] with strongly peaked eigenmode structure that decays within a single drift well.In the local (in l) theory of toroidal ITG modes [14,15,19], the CG is set by the drive parameter κ d = R eff /L T .This threshold can be computed for general parameter b = k ⊥ 2 ρ 2 by solving the local dispersion relation In realistic geometry, the threshold is controlled by the extent of drift wells, i.e. the parallel connection length L ∥ , but this can be related to finite Larmor radius (FLR) stabilization as follows: Note that a toroidal mode must have a drift frequency ω d ∝ k α that exceeds the parallel transit rate k ∥ v T ∼ πv T /L ∥ .Although this can always be satisfied by choice of k α , the increase of k α comes at the cost of increasing b as L ∥ is reduced.Thus, to estimate the critical gradient, we simply determine the minimum value of b for which the resonance condition is satisfied, namely that for which ω d ∼ πv T /L ∥ , yielding b min = (πa|∇α|R eff /L ∥ ) 2 , and with F (b) defined as above.R eff is determined by the peak of a quadratic fit to K d and L ∥ by the distance between points where the sign of K d reverses [12] within a drift well of "bad" curvature, while a|∇α| is evaluated at the center of the fitted drift well, effectively approximating it as a constant.In the small-b limit, F (b) is dominated by the constant term 2.84, close to the value of 2.66 in Eqn. 3, and as found other works [12,14,29] for the case τ = 1.In the large b limit, we find, ignoring the small constant ≃ 0.04, that the formula effectively predicts a/L T,crit ∼ a 2 |∇α|/L ∥ .Perhaps unsurprisingly, then, the threshold for toroidal modes in this regime is dominated by both the gradient of the bi-normal coordinate (linked to expansion of surfaces as well as local magnetic shear [30]) and the parallel connection length, which can be reduced by simply increasing the number of field periods in a configuration.Indeed, an experimental realization of this strategy can already be seen in the 10 field-period LHD heliotron, whose favorable ITG turbulence properties relative to W7-X have been demonstrated [31].More generally, equations ( 3) and ( 6) can now be used to rapidly estimate the absolute and toroidal ITG thresholds on a given magnetic field line.
Optimization results.-Weuse the SIMSOPT software framework [32] to generate two QHS vacuum stellarator configurations.Each stellarator magnetic field is described by a boundary surface given in the Fourier representation R(ϑ, ϕ) = m,n R m,n cos(mϑ − n f p nϕ), Z(ϑ, ϕ) = m,n Z m,n sin(mϑ − n f p nϕ). Optimization proceeds by treating the above-mentioned Fourier coefficients as parameters and varying them in a series of steps in order to find a least-squares minimization of the specified objective function f , increasing the number of boundary surface modes with each step.Both optimizations used the "warm start" configuration from SIMSOPT with approximate QHS and n f p = 4. Global zero-β equilibria are constructed at each iteration by running the VMEC [33] code, which solves the MHD equations using an energy-minimizing principle, setting the angular resolution to be M pol = N tor = 7.
For the first result, which we call "HSK", f = f QS + (A−4.10) 2 +f abs , where f QS is the quasisymmetry residual defined in [34] for QHS with n f p = 4 at the surfaces (r/a) 2 = [0.1,0.2, 0.3, 0.4, 0.5], f abs = (a/L T,crit,abs − 2.00) 2 is the critical gradient evaluated at the flux tube [(r/a) 2 = 0.5, α = 0], and f A = (A − 4.10) 2 is the aspect ratio target with A = R/a the aspect ratio output by VMEC.The boundary modes varied for the three optimization steps went up to m pol = n tor = [3,5,6].Linear flux tube gyrokinetic simulations with GENE [35] reveal that HSK has the largest critical gradient of any stellarator that we know of, a/L T,crit,abs = 1.75 [fig.(1)], as well as relatively low nonlinear ion heat fluxes above that threshold.Further details of HSK and the nonlinear simulations are presented in [30].The caveat is that the large "bad" curvature of destabilizing sign for HSK (a small R eff linked to enhancement of |∇α|) produces a vacuum magnetic hill, rendering it Mercier unstable [36,37] at all values of β tested.
In the second optimization we choose n f p = 6 and target the toroidal ITG threshold in the hopes of reducing turbulent transport while preserving MHD stability.The objective function is where f A = [Θ(A − 7.50)] 2 is the aspect ratio penalty, Θ(x) is defined to be xH(x) with H(x) the Heaviside step function, f QS is again the quasisymmetry residual but with n f p = 6, and [eq. ( 6)] is taken at the surface r/a = 0.5 and summed over the field lines α = [0, π/8, π/4], with each field line extending for 8 poloidal turns, in order to sample the surface.The vacuum magnetic well penalty is 2 with the surfaces (r/a) 2 = [0, 0.1, ..., 0.9] targeted and V ′′ (r) the second derivative of the flux surface volume with radius.We calculate the residual f QS on the surfaces (r/a) 2 = [0.1,0.2, ..., 0.9], and target the axis and boundary iota via , with ι the rotational transform.The "warm start" file was first optimized for increased |∇α| on the outboard side (see e.g.[30,38]), and increasing n fp to 6.The final optimization proceeded in two steps, with the boundary Fourier coefficients varied up to m pol = n tor = [3,4].
ITG turbulence.-Toevaluate the performance of QSTK in vacuum with regard to ITG turbulence we run full-surface nonlinear electrostatic gyrokinetic simulations using the GENE code [35,44] in comparison with the HSX stellarator.We assume adiabatic electrons, zero density gradient, T e = T , and temperature gradients a/L T = [1.5, 2.0, 2.5, 3.0] at half radius.In Figs. 4 (a)-(b) we plot the ion heat fluxes in gyro-Bohm units times A for each configuration, to adjust for the dependence of energy confinement time on aspect ratio implied by the   [49], which relies on the isomorphism between quasisymmetry and axisymmetry, was found to produce a similar bootstrap current profile for QSTK, although the results were not in as good agreement as in cases with "precise" quasisymmetry [34].Alpha losses are slightly reduced to roughly 4% at half radius for the case with bootstrap current [Fig.3(c)].Evaluations [50] of the Mercier criterion indicate that, for the configuration with pressure and bootstrap current included, the Mercier criterion is satisfied [fig.5(b)] and the bootstrap current produces a stabilizing up-shift in D M [51].Certain radii become Mercier unstable for larger values of β.We also apply the coil optimization features of SIMSOPT to produce the magnetic field of QSTK, finding that the relative maximum field error can be reduced to 5.1% and the relative mean error to 1.1% with four unique coils (48 coils in total) while penalizing coil length.The initial coil and MHD studies show that QSTK has potential for finite-β (reactor-relevant) operation scenarios.
Discussion and conclusions.-Othermicroturbulence, such as trapped electron modes and electron tempera- ture gradient (ETG) turbulence, remain to be studied in QSTK.ETG turbulence is likely to benefit from the ITG optimization for QSTK both linearly (from the isomorphism with ITG modes [14]) and nonlinearly (from the short connection length [52]) with regard to ETG losses.At high β values, turbulence is expected to transition from ITG to kinetic ballooning mode turbulence [25,53].This physics is delegated to a separate publication, in the framework of a QSTK reactor study.Despite these open challenges, the present work demonstrates the possibility of modifying the current design of magnetic flux surfaces, toward a drastic suppression of turbulence in the parameter range usually encountered in modern stellarator experiments.Our method integrates salient physics properties, such as good MHD stability and particle confinement, low neoclassical transport, and bootstrap current, together with the feasibility of modular coils.The QSTK configuration introduced here is thus a contender for a future compact fusion reactor based on the stellarator concept.
FIG. 1. Above: the boundary surface of HSK showing the contours of B, with one half of a field period removed.Below: figure (1) of Roberg-Clark et al. (2022) plotting critical gradients found with GENE versus the model eqn.(3), with the point for HSK added.

FIG. 2 .
FIG. 2. Boundary surface of QSTK (one field period removed), showing contours of B in color.Inset: cuts at constant toroidal angle of the boundary surface in the cylindrical (R,Z) plane.

FIG. 4 .
FIG. 4. Full-surface nonlinear gyrokinetic simulations of ITG turbulence comparing HSX to QSTK.(a) Log plot of ion heat flux multiplied by the respective aspect ratio of each configuration.(b) ITG density fluctuations ñ/n0 in one field period of HSX with a/LT = 2, plotted in Boozer toroidal (ϕ) versus poloidal (θ) angles.(c) Same as in (b) but for QSTK.

FIG. 5 .FIG. 6 .
FIG. 5. Bootstrap current and MHD stability of QSTK with an applied pressure profile (see text) at volume-averaged β = 1.65%.The dashed curves correspond to the case without boostrap current and the solid curves to the case with current.(a) Enclosed toroidal bootstrap current.(b) Rotational transform profile.(c) Mercier stability criterion with a dotted line at DM = 0.