Ion-temperature-gradient stability near the magnetic axis of quasisymmetric stellarators
IIon-temperature-gradient stability near the magneticaxis of quasisymmetric stellarators
R. Jorge, M. Landreman
Institute for Research in Electronics and Applied Physics, University of Maryland,College Park, MD 20742, USAE-mail: [email protected]
Abstract.
The stability of the ion-temperature gradient mode in quasisymmetricstellarators is assessed. This is performed using a set of analytical estimates togetherwith linear gyrokinetic simulations. The peak growth rates, their corresponding realfrequencies and wave-vectors are identified. A comparison is made between a first-order near-axis expansion model and eleven realistic designs obtained using numericaloptimization methods. It is found that while the near-axis expansion is able to replicatethe growth rates, real frequencies and perpendicular wave-vector at the inner core(both using simplified dispersion relations and first-principle gyrokinetic simulations),it leads to an overestimation of the growth rate at larger radii. An approximate analyticsolution of the ITG dispersion relation for the non-resonant limit suggests growth ratescould be systematically higher in quasi-axisymmetric (QA) configurations comparedto quasi-helically (QH) symmetric ones. However except for very close to the axis,linear gyrokinetic simulations do not show systematic differences between QA and QHconfigurations.
1. Introduction
The feasibility of magnetic confinement fusion relies on a delicate balance betweenupholding a hot enough plasma in the core of the reactor and the inevitable transport ofparticles and energy to the wall. Such transport is termed neoclassical when associatedwith geometrical effects in standard collisional transport theories, and termed turbulentwhen associated with perturbations in the plasma driven by underlying instabilities. Thereduction of neoclassical transport to acceptable levels has been, to an extent, attainedwith the tokamak and omnigeneous (such as quasisymmetric and quasi-isodynamic)stellarator designs. The reduction of turbulent transport, however, has been a focus oftheory and modelling of magnetic confinement plasma physics in the past few decades. Inthe core, two microinstabilities are thought to be the main drive of turbulent transport:the Ion Temperature Gradient (ITG) mode and the Trapped Electron Mode (TEM)[1, 2]. In this work, we focus on the question of linear ITG stability and its parameterdependence in the core of a particular class of omnigeneous stellarator designs, namelyquasisymmetric stellarators. a r X i v : . [ phy s i c s . p l a s m - ph ] F e b TG stability near the magnetic axis B = | B | varies on a surface of constant toroidal flux only through a fixed combination of thepoloidal θ and toroidal ϕ Boozer angles B = B ( ψ, M θ − N ϕ ) for some integers M and N [3, 4, 5]. There are currently two ways of developing quasisymmetric stellaratordesigns: analytically, using the near-axis expansion [6, 7, 8] or a near-axisymmetryformulation [9]; or numerically, using one of the stellarator optimization methodsavailable [10, 11, 12]. In this work, the microstability of quasisymmetric designs near themagnetic axis will be assessed using both methods, and the advantages and shortcomingsof using either method are highlighted.Microinstabilities are driven by the presence of free energy associated withtemperature and density gradients and are, therefore, ubiquitous in tokamaks andstellarators. These develop locally in a plasma at the gyroradius scale, as opposedto macroinstabilities that affect the plasma globally. Microinstabilities have regions offavorable development in the device, usually assessed by their relative "good" or "bad"curvature properties, evaluated using the sign and magnitude of the magnetic driftfrequency ω d (cid:39) T B × ∇ B · k ⊥ / ( eB ) where T is the background temperature, e theelementary charge and k ⊥ the wave-vector of the mode perpendicular to the magneticfield B . If this frequency is positive (or, more accurately, if it has the same sign asthe diamagnetic drift-frequency), trapped particle modes are likely to be unstable [13].While ω d is positive in the outboard side of tokamaks, rendering unstable the regionwhere trapped particles spend most of their time, in stellarators, the locations of badcurvature need to be evaluated numerically and are not necessarily in the outboardregion. As a result, many studies aimed at assessing the influence of magnetic fieldgeometry on plasma stability and turbulence have been performed, particularly usinggyrokinetic theory [2]. The freedom associated with the length and location of the badcurvature region has also allowed for stellarator optimization criteria based on the signof ω d and the non-overlap of trapped and "bad" curvature regions [14]. While the firstnon-axisymmetric gyrokinetic studies were performed using the linear eigenvalue FULLcode [15, 16], only recently have nonlinear simulations been performed [17] together withstellarator optimization efforts aimed at reducing turbulent transport [18].In this work, the ITG instability is analyzed first by incorporating the near-axis expansion into several known analytical estimates and then by performing lineargyrokinetic simulations for a set of eleven quasisymmetric designs using the GS2 code[19]. This instability can either be destabilized by the presence of an equilibriumtemperature gradient which alters the temperature of a perturbation travelling at thesound speed that is in phase with a perturbed E × B drift, usually called the slab branchof the mode, or by the coupling of the grad- B and curvature drifts with the equilibriumtemperature gradient, usually called the toroidal branch of the mode [20]. We mentionthat, although this work focuses on the properties of the linearized gyrokinetic equation, TG stability near the magnetic axis
2. Physical Model and Geometry
In this work, we employ the gyrokinetic approximation, i.e., we assume that the plasmais strongly magnetized, and that fluctuations have small amplitude and low frequency[30]. These assumptions are, in general, satisfied in magnetic confinement nuclear fusiondevices. In here, we find that the particle (electron and ion) Larmor radius, ρ = v th / Ω with v th = (cid:112) T /m the thermal velocity,
Ω = qB/m the cyclotron frequency, T thetemperature, m the mass, q the charge and B the magnetic field strength, is muchsmaller than typical length scales of the macroscopic equilibrium, L . Furthermore, thefrequency of fluctuations, ω , is much smaller than the cyclotron frequency, Ω . This allowsus to average the particle’s equations of motion, effectively turning their Larmor orbitinto a guiding-center one via the gyroaverage operator (cid:104) ... (cid:105) . Using the small fluctuationassumption and focusing on an electrostatic regime, we split the distribution function TG stability near the magnetic axis f = F + δf , and the electrostatic potential φ = δφ , into an equilibrium(denoted with a subscript ) and a fluctuating (denoted with δ ) part, where δf (cid:28) F and qφ/T (cid:28) . Furthermore, the fluctuating distribution function is split into anadiabatic and non-adiabatic part as δf = h − ( qφ/T ) F . For the equilibrium, we take F to be a Maxwell-Boltzmann distribution with thermal velocity v th and density n ,and both the background density n = n ( ψ ) and temperature T = T ( ψ ) are taken to befunctions of the magnetic toroidal flux ψ only, with characteristic lengths L n and L T ,respectively. Due to its importance in assessing regions of stability in the gyrokineticspace of parameters, we define η = L n /L T as the ratio between the two characteristiclengths. In the following, for simplicity, we denote (cid:104) φ (cid:105) and (cid:104) h (cid:105) as φ and h , respectively,we assume that electrons are adiabatic with density n and temperature T e and thata single ion species is present, with its gyroaveraged distribution function given by h ,density n and temperature T i = τ T T e .The nonlinear gyrokinetic equation solved by GS2, in the linear electrostatic limit,can be written as [31] ∂h∂t + ( v (cid:107) b + v d ) · ∇ h + b × ∇ φB · ∇ F = (cid:104) C (cid:105) + qF T ∂φ∂t . (1)The particle’s grad-B and curvature drifts are contained in v d = b Ω × (cid:18) v (cid:107) b · ∇ b + v ⊥ ∇ BB (cid:19) . (2)We note that, in Eq. (1), the spatial derivatives are taken at constant energy E = mv / and magnetic moment µ = mv ⊥ / B and collisions are modeled via the collision operator C . The system of equations is closed via the quasineutrality equation which, in the limitof adiabatic electrons, reads (cid:90) d v (cid:104) h (cid:105) r = n (1 + τ T ) qφT , (3)where (cid:104) h (cid:105) r denotes a gyroaverage performed at constant position vector r as opposedto (cid:104) h (cid:105) which is performed at constant guiding-center locaion R . A relation betweenthe two can be found by Fourier decomposing h into a parallel and perpendicularwave-vectors k = k (cid:107) b + k ⊥ with b = B /B , allowing us to write, in Fourier space, (cid:104) h (cid:105) r = J ( k ⊥ v ⊥ / Ω) (cid:104) h (cid:105) .In order to perform a stability analysis, we write Eq. (1) in the field-line followinglimit. Additionally, due to the high temperature of the plasma in the core, we look atthe collisionless limit ν e (cid:28) ω where ν e is the electron-ion collision frequency, and set C = 0 . In the field-line following limit, we write the perturbed quantities h and φ as h = ˆ h ( l ) e i ( S − ωt ) , (4)where l is the distance along a magnetic field line and b · ∇ S = 0 , leading us to define k ⊥ ≡ ∇ S . In this limit, and normalizing φ by q/T , Eq. (1) can be written as iv (cid:107) b · ∇ ˆ g + ( ω − ˜ ω d )ˆ g = ˆ φF ( ω − ˜ ω ∗ ) J (cid:18) v ⊥ k ⊥ Ω (cid:19) , (5) TG stability near the magnetic axis ˜ ω d = k ⊥ · v d the magnetic drift-frequency. When a low β limit is taken, ˜ ω d can bewritten as ˜ ω d = ω d ( v (cid:107) + v ⊥ / /v th with ω d its velocity-independent counterpart. Thevelocity-dependent diamagnetic frequency ˜ ω ∗ is given by ˜ ω ∗ = ω ∗ (cid:20) η (cid:18) v v th − (cid:19)(cid:21) , (6)where ω ∗ = b × k ⊥ · ∇ F / ( m Ω T F ) usually written as ω ∗ = ( T k α /q ) d ln n/dψ when aClebsch representation for the magnetic field is used.As a final note, we mention here the normalizations used to output GS2 quantities.Lengths are normalized to a fixed length a (a measure of the minor radius of the device),time is normalized to a/v th and perturbed quantities are scaled up by a/ρ th . In theaxisymmetric case, a is half of the diameter of the last closed flux surface measuredat the elevation of the magnetic axis. Here, however, we identify a with the effectiveminor radius computed by VMEC [32], named Aminor_p in its wout*.nc file, defined as π (Aminor_p ) = (1 / π ) (cid:82) π S ( φ ) dφ where φ is the toroidal angle associated with thecylindrical coordinate system ( R, Z, φ ) ( φ is not to be confused with the electrostaticpotential) and S ( φ ) is the area of the outer surface’s cross-section in the R − Z plane. We now rewrite the differential operators in the linearized gyrokinetic equation, Eq. (5),in terms of the straight-field-line magnetic coordinates commonly referred to as Boozercoordinates [33] and employ the Garren-Boozer near-axis expansion formalism [6, 34].We note that such formalism uses an inverse coordinate approach, writing r as a functionof Boozer coordinates, as opposed to formalisms using a direct coordinate approach,which write the toroidal flux as a function of the spatial coordinates r [35, 36, 37].First, we express the magnetic field B using a Clebsch representation of the form B = ∇ ψ × ∇ α, (7)where πψ is the toroidal flux and α is a field line label. Taking as a spatialcoordinate system the set ( ψ, α, z ) with z any quantity independent of ψ and α , theeight independent geometrical quantities Q needed to solve Eq. (5) are [26] Q = (cid:8) B, b · ∇ z, |∇ ψ | , |∇ α | , ∇ ψ · ∇ α, ( b × ∇ B ) · ∇ α, ( b × ∇ B ) · ∇ ψ, ( b × κ ) · ∇ α } . (8)From here onward, we employ Boozer coordinates ( ψ, θ, ϕ ) with θ and ϕ the poloidaland toroidal angles, respectively, and choose z = ϕ . For convenience, a helical angle ϑ = θ − N ϕ with N an integer is introduced, such that the magnetic field B can bewritten as B = ∇ ψ × ∇ ϑ + ι N ∇ ϕ × ∇ ψ, (9) = β ∇ ψ + I ∇ ϑ + ( G + N I ) ∇ ϕ, (10) TG stability near the magnetic axis ι N = ι with ι the rotational transform and I, G and ι are constants on ψ surfaces,i.e., I = I ( ψ ) , G = G ( ψ ) and ι = ι ( ψ ) . A relation between α and the angles ϑ and ϕ can be found by comparing Eq. (7) and Eq. (9), yielding α = ϑ − ι N ϕ = θ − ιϕ .The geometry coefficients Q in Eq. (8) are evaluated at lowest order in (cid:15) = r/ R ,where r is defined as π | ψ | = πr B and R a scale length representing the majorradius of the device (e.g., the inverse of the maximum axis curvature). Furthermore,we focus on quasisymmetric magnetic fields, i.e., we use magnetic geometries where B = B ( ψ, M ϑ − N ϕ ) with M and N integers. As shown in Ref. [26], this yields thefollowing set of geometric coefficients B = B (1 + rη cos ϑ ) , (11) b · ∇ z = s G πL (1 + rη cos ϑ ) , (12) |∇ ψ | = r B η κ (cid:2) η sin ϑ + κ (cos ϑ − σ sin ϑ ) (cid:3) , (13) |∇ α | = 1 r η κ (cid:2) η cos ϑ + κ ( σ cos ϑ + sin ϑ ) (cid:3) , (14) ∇ ψ · ∇ α = s ψ B η κ (cid:0)(cid:2) η + κ (cid:0) σ − (cid:1)(cid:3) sin 2 ϑ − κ σ cos 2 ϑ (cid:1) , (15) ( b × ∇ B ) · ∇ α = s ψ r B η cos ϑ, (16) ( b × ∇ B ) · ∇ ψ = rB η sin ϑ, (17) ( b × κ ) · ∇ α = s ψ r η cos ϑ. (18)In the set of Eqs. (11) to (18), B is the constant magnetic field on-axis, r = (cid:112) ψ/B isan effective minor radius, η is a constant that is related to the elongation of the ellipticalflux-surfaces, κ = κ ( ϕ ) is the axis curvature, s ψ = sgn ( ψ ) , s G = sgn ( G ) , and σ = σ ( ϕ ) isthe solution of dσdϕ + ( ι − N ) (cid:18) η κ + 1 + σ (cid:19) − s G Lη πκ (cid:18) I B − s ψ τ (cid:19) = 0 , (19)where τ = τ ( ϕ ) is the axis torsion, ι is the on-axis rotational transform, I is the lowestorder toroidal current I (cid:39) I r and L is the total length of the magnetic axis curve r ( φ ) L = (cid:90) π (cid:12)(cid:12)(cid:12)(cid:12) d r dφ (cid:12)(cid:12)(cid:12)(cid:12) dφ. (20)Alternatively, due to Eq. (11), η can be regarded as the constant parameter that reflectsthe magnitude of variation in B .We now introduce dimensionless x ( ψ ) and y ( α ) coordinates that correspond toscaled versions of ψ and α : x = dxdψ ( ψ − ψ ) , (21) TG stability near the magnetic axis y = dydα ( α − α ) , (22)with dx/dψ and dy/dα constant. Here, we take ψ = α = 0 and define dx/dψ and dy/dα such that x reduces to the usual minor radius in the cylindrical limit and that theproduct xdx/dψ equals the inverse of a reference magnetic field, B ref . More specifically,we take x = ar/r a and y = αar/r a with r a the value of r at the VMEC boundary and a = Aminor_p the effective minor radius. The linearized gyrokinetic equation, Eq. (5),is then Fourier decomposed along the x and y directions with k ⊥ = k x ∇ x + k y ∇ y theperpendicular wave-vector.An example first-order near-axis stellarator geometry is shown in Fig. 1, namelythe quasi-axisymmetric design of Section 5.1 of Ref. [7]. A field line with r = 0 . , α = 0 and − π < ϑ < π is shown, together with the location of ϑ = 0 . A poloidalprojection of a longer field-line with − π < ϑ < π and ten poloidal cross-sections ofthe flux surface are shown in the top-right of Fig. 1. The flux-tube approach employedin GS2 to simulate plasma conditions in the field-line following limit can then thoughtof as a series of field-lines similar to the one in Fig. 1 with a finite (but small) extensionalong the radial ( ψ ) and poloidal ( α ) directions.The stellarator geometries used in this study, apart from a near-axis benchmarkcase, consist of eleven quasisymmetric stellarator designs. These were developed byseveral independent research teams using different optimization codes, with varyingdegrees of quasisymmetry. By order of decreasing number of field periods, these arethe NZ1988 design of Ref. [38], the Drevlak design of Ref. [39] recently developedat the Max Planck Institute for Plasma Physics in Greifswald, Germany, HSX [40],the KuQHS48 design of Ref. [41], the WISTELL-A configuration of Ref. [42] recentlydeveloped at the University of Wisconsin-Madison, USA, NCSX (configuration LI383,[43]), ARIES-CS configuration N3ARE [44], the QAS2 configuration of Ref. [45] with avanishing on-axis current, ESTELL [46], CFQS [47] and the Henneberg design of Ref.[48]. The properties of each stellarator design are listed in Table 1 of Ref. [26], togetherwith the B and η constants used for each configuration and a comparison betweenthe original and the resulting near-axis geometric coefficients. The original coefficients Q are obtained from the VMEC file corresponding to each design using the geometrymodule of the stella code [49].Finally, we review the main input and output parameters of the model equationswhen combined with the near-axis geometry described above. For each configuration, weuse its axis shape to find the axis curvature κ and torsion τ . The parameters B and η arefound using the fitting procedure described in Ref. [26] allowing us to find the σ functionfor each configuration, hence to fully specify the geometry coefficients in Eqs. (11)to (18). The physical input parameters used by GS2 to solve the linear gyrokineticequation, Eq. (5), are the x and y components of the perpendicular wave-vector, k x and k y , respectively, and the physical parameters related to the particle species. We notethat, at α = 0 , choosing a flux tube with k x = 0 yields the largest growth rate [50]. TG stability near the magnetic axis - - R [ meters ] Z [ m e t e r s ] Figure 1.
Top left: top view of a first-order quasi-axisymmetric surface of constanttoroidal flux ψ (Section 5.1 of Ref. [7]). Colour indicates B on the outermost fluxsurface. An example field-line with − π < ϑ < π , α = 0 and r = 0 . is shown inblack together with the location of ϑ = 0 . Top right: ten poloidal planes and theirrespective axis location, together with a field line with − π < ϑ < π and α = 0 .Bottom: side view of the toroidal flux surface. Therefore, in the following, all simulations are performed at α = k x = 0 . As electronsare considered adiabatic, the input characteristic gradient lengths consist only of theion density, L n , and temperature, L T , lengths. The atomic number is set to Z = 1 . Thetemporal and spatial resolution parameters are set using convergence scans.
3. Models of Near-Axis ITG stability
In this section, we use the geometric coefficients Q in the set of Eqs. (11) to (18)in order to simplify the characteristic ITG frequencies ω ∗ and ω d , together with theperpendicular wave-vector k ⊥ . Then, we plug these simplified expressions into knownlimits of the gyrokinetic equation commonly employed in ITG stability studies. Theresulting models allow for a more intuitive understanding of the properties of the ITGmode near the magnetic axis and are used in the subsequent sections to interpret the TG stability near the magnetic axis k α = ark y /r a , thediamagnetic frequency ω ∗ can be written as ω ∗ Ω = aL n ρ th r a k y ρ th , (23)where ρ thi = v th / Ω and L n = d ln n/dr (and similarly for L T ), while the drift-frequencyreads ω d Ω = s ψ ar a ηρ th k y ρ th cos ϑ. (24)In the following, we further assume that the most unstable linear modes peak at k ψ = 0 and write the square of the normalized perpendicular wave-vector b as b = k ⊥ ρ th = a r a k y ρ th (cid:20) η κ cos ϑ + κ η ( σ cos ϑ + sin ϑ ) (cid:21) . (25)Taking z = ϑ , the parallel derivative operator can be written as b · ∇ φ = ∇ (cid:107) φ = s G πL ι N (1 + rη cos ϑ ) ∂φ∂ϑ . (26)Furthermore, by defining the normalized quantities ˜ k y = k y ρ th , ˜ L n = L n π/L , ˜ η = ηL/ π , ˜ κ = κ/η , ˜ ω = ω/ω t with the transit frequency ω t = 2 πv th /L (similarlyfor ω d and ω ∗ ), we rewrite the ITG parameters above as ˜ ω ∗ = ˜ k y ˜ L n , (27)and ˜ ω d = s ψ ˜ k y ˜ η cos ϑ, (28)and b = ˜ k y ˜ κ (cid:2) cos ϑ + ˜ κ ( σ cos ϑ + sin ϑ ) (cid:3) , (29)where the approximation a (cid:39) r a was used. This reduces the set of dimensionlessparameters entering the linearized gyrokinetic equation to τ T , ˜ k y , ˜ L n , η, ˜ η, ˜ κ and σ . Wenote that for stellarator-symmetric designs σ ( ϕ = 0) = 0 and that according to Eq. (19) σ (cid:48) (0) = 1( ι − N ) dσdϕ (cid:12)(cid:12)(cid:12)(cid:12) = 2˜ κ (0) [ ˜ I − s ψ ˜ τ (0)]( ι − N ) [1 + ˜ κ (0) ] , (30)where ˜ I = I L/ (2 πB ) is the normalized current and ˜ τ = τ L/ π is the normalizedtorsion.One limit in which further analytic simplification is possible is the case in whichthe eigenfunction is localized near ϑ = 0 . We can then expand the ITG parameters ω d and b around ϑ = 0 as ˜ ω d = ˜ ω d (cid:18) − ϑ (cid:19) , (31) TG stability near the magnetic axis b = b (1 + sϑ ) , (32)where we find ˜ ω d = s ψ ˜ k y ˜ η , b = ˜ k y / ˜ κ , the effective shear parameter s = − (cid:34) κ (cid:48)(cid:48) (0)( ι − N ) ˜ κ (0) − ˜ κ (0) (cid:18) σ (cid:48) (0) ι − N (cid:19) (cid:35) , (33) Following Ref. [51], we consider the limit where the transit frequency ω t is small whencompared with the mode frequency ω . Furthermore, we adopt the non-resonant limitwhere ω/ω T ∗ = δ (cid:28) and ω d /ω ∼ ω t /ω ∼ k ⊥ ρ ∼ δ , yielding the following secondorder differential equation for the electrostatic potential [52] ω T ∗ v th ω ∇ (cid:107) φ = (cid:20) τ T + (1 − bη ) ω ∗ ω + ω d ω T ∗ ω (cid:21) φ, (34)where ω T ∗ = ηω ∗ and ∇ (cid:107) φ (cid:39) ( ι − N ) (4 π /L ) ∂ φ/∂ϑ . In the highly localizedapproximation of Eqs. (31) and (32), Eq. (34) turns into a Schrödinger equation ofthe Weber type whose bound solutions have the form φ ∼ H n ( √ δϑ ) exp( δϑ / where H n is a Hermite polynomial of order n = 0 , , , ... and the width δ is given by δ = − ˜ ω ˜ k y (cid:32) s ψ ˜ η + 2˜ k y s ˜ κ ˜ ω (cid:33) , (35)where ω is the root of the dispersion relation τ T ˜ L n ˜ k y η ˜ ω + (cid:32) − η ˜ k y ˜ κ (cid:33) ˜ ω η + s ψ ˜ k y ˜ η ˜ ω + (cid:18) n + 12 (cid:19) δ + = 0 , (36)for which Re [ δ ] > (denoted as δ + ). We note that, in the dispersion relation above, theeffective shear parameter s only contributes to δ as an s factor yielding eigenfunctionsthat are symmetric with respect to s .An example numerical solution for the dispersion relation in Eq. (36) with s =0 . , ˜ η = 1 and ˜ κ = 1 is shown in Fig. 2. There, the peak growth rate γ and thecorresponding value of ˜ k y and δ + are shown for a range of L n and L T values satisfying η = L n /L T > . It is seen that, in general, the ITG growth rate γ and its correspondingmode width δ + increase with increasing a/L n and a/L T while the opposite is found for k y . Besides a scan over a/L n and a/L T , the dispersion relation in Eq. (36) can be usedto assess the impact of s and η b (which are quantities directly related to the near-axisexpansion parameters) in the stability of the ITG mode. In Fig. 3 a parameter scan over s and η b is shown for a fixed value of η = 3 and a/L n = 1 . It is seen that the results aresymmetric with respect to the sign of η and that all four parameters increase as either s or η are increased. TG stability near the magnetic axis γ k y ρ th ω δ + Figure 2.
Solution of the ITG dispersion relation Eq. (36) with η > for the case s = 0 . , ˜ η = 1 and ˜ κ = 1 , namely the peak growth rate γ (top left), the corresponding k y ρ th (top right), real frequency ω (bottom left) and the eigenfunction width δ + (bottom right). The frequencies γ and ω are normalized to the transit frequency ω t = 2 πv th /L . Finally, we apply the dispersion relation in Eq. (36) to the eleven quasisymmetricdesigns used in this study. For each design, the value of η is read from Table 1 in Ref.[26], the axis parameters are read from each VMEC file (used to compute the curvatureand torsion) and σ is calculated using Eq. (19), yielding the values of ι in Table 1 of Ref.[26]. The values of the effective shear parameter s for each design is shown in Fig. 4.It is found that quasi-helically symmetric stellarators have, in general, values of s thatare one to two orders of magnitude above the ones of quasi-axisymmetric stellarators.We also show in Fig. 4 the resulting growth rate estimate and the corresponding valuesof ηL/ π and s . The designs with ηL/ π < are quasi-axisymmetric, while the oneswith ηL/ π > are quasi-helically symmetric. The results in Fig. 4 suggest thatquasi-helically symmetric configurations are prone to having higher values of ηL/ π , TG stability near the magnetic axis - - η s γ - - η s k y ρ th - - η s ω - - η s δ + Figure 3.
Solution of the ITG dispersion relation Eq. (36) for the case η = 3 and a/L n = 1 , namely the peak growth rate γ (top left), the corresponding k y ρ th (topright), real frequency ω (bottom left) and the eigenfunction width δ + (bottom right).The frequencies γ and ω are normalized to the transit frequency ω t = 2 πv th /L while η is normalized to L/ π . lower values of s and lower growth rates (normalized to ω t ) when compared with quasi-axisymmetric designs. We now show how a bounce-averaged metric Q bounce , which takes into account the roleof good and bad curvature regions in the particle motion, can be simplified using thenear-axis expansion. This parameter, although more focused on TEM rather than ITG,allows us to show how metrics based on bounce-averaged quantities, in the near-axisexpansion formalism, can be reduced to one-dimensional elliptic integrals. Here, we TG stability near the magnetic axis .
00 0 .
25 0 .
50 0 .
75 1 .
00 1 .
25 1 .
50 1 .
75 2 . ηL/ π − l og | s | WISTELL-A NZ1988HSXKuQHS48 DrevlakNCSXARIES-CS QAS2ESTELLCFQSHenneberg Designs (QH) NZ1988 ˜ γ = 0.21(QH) HSX ˜ γ = 0.22(QH) KuQHS48 ˜ γ = 0.22(QH) Drevlak ˜ γ = 0.27(QH) WISTELL-A ˜ γ = 0.33(QA) ESTELL ˜ γ = 0.56(QA) ARIES-CS ˜ γ = 1.10(QA) CFQS ˜ γ = 1.16(QA) NCSX ˜ γ = 1.26(QA) Henneberg ˜ γ = 2.16(QA) QAS2 ˜ γ = 2.75 . . . . . ˜ γ Figure 4.
Effective shear parameter s and η for each of the eleven quasisymmetricdesigns considered in this study computed using Eq. (33). The designs with hollowcircles are quasi-axisymmetric, while the ones with filled circles are quasi-helicallysymmetric. The color represents the growth rate γ computed from the dispersionrelation in Eq. (36) with ˜ L n = 0 . and η = 2 . . employ the metric used in Ref. [53] Q bounce = − (cid:90) /B min /B max ω d ( λ ) dλ, (37)where ω d ( λ ) is the bounce-averaged drift-frequency ω d ( λ ) = (cid:82) l l ˜ ω d dl/v (cid:107) (cid:82) l l dl/v (cid:107) . (38)The minimization of this metric aims at distributing a large fraction of trapped particlesover a favorable bounce-averaged curvature. Although Eq. (37) has been used in Ref.[53] to stabilize TEMs, similar metrics are often employed in the analysis of tokamakand stellarator microstability [54].In the near-axis expansion formalism, at first order in r , the metric Q bounce inEq. (37) can be simplified by noting that v (cid:107) = ± v √ − λB with B = B (1 + rη cos ϑ ) .The resulting integral in l is rewritten by parametrizing the field-line coordinate l using dϑ/dl = b ·∇ ϑ = | ι N | B/ ( G + ιI ) (cid:39) (2 π/L ) | ι N | (1+ rη cos ϑ ) , yielding the bounce points l and l = − l at ϑ = − ϑ = − /a ) with a = − rηλB / ([1 − λB (1 + rη )] > .By performing a variable substitution x = a sin ϑ/ and assuming s ψ = 1 , the followingexpression for Q bounce is found Q bounce = 8 √ k y ρ th ηv v th (cid:90) [ B (1 − rη )] − [ B (1+ rη )] − dλ (cid:112) − λB (1 + rη ) Π c (1 /c, t, /a )Π( t, /a ) , (39)with Π c and Π the elliptic integrals Π c (1 /c, t, /a ) = (cid:90) (1 − x /a ) (cid:112) − x /c √ − x (cid:112) − x /a dx − tx , (40) TG stability near the magnetic axis r η q b o un ce Exact formulaQuadratic fit .
00 0 .
25 0 .
50 0 .
75 1 .
00 1 .
25 1 .
50 1 .
75 2 . ηL/ π . . . . . . . . . Q b o un ce WISTELL-ANZ1988HSXKuQHS48 DrevlakNCSXARIES-CSQAS2ESTELLCFQSHenneberg
Figure 5.
Left: Analytical form and the corresponding quadratic fit . x − . x forthe function q bounce entering the heat flux proxy Q bounce in Eq. (42). Right: Q bounce normalized to k y ρ th v /B v th as a function of ηL/ π for the eleven configurationsconsidered in this study and r = 0 . . The designs with hollow circles are quasi-axisymmetric, while the ones with filled circles are quasi-helically symmetric. and Π( t, /a ) = (cid:90) √ − x (cid:112) − x /a dx − tx , (41)respectively. In Eq. (39), t = 2 rη/ [(1+ rη ) a ] and c = [2 − λB (1+ rη )] / [1 − λB (1+ rη )] .Noting that the integrand of Eq. (39) is a function of rη and λB only, we rewrite Q bounce by turning the integral over λ to an integral over λB , yielding. Q bounce = 8 √ k y ρ th ηv v th B q bounce ( rη ) . (42)As shown in Fig. 5, the function q bounce can be reasonably approximated by a quadraticfunction of rη , leading to the following approximate model for Q bounce close to themagnetic axis Q bounce (cid:39) . k y ρ th ηv B v th rη [1 − . rη ) ] . (43)The positive sign of Q bounce in Eq. (43) suggests that, for quasisymmetric stellarators,the majority of trapped particles experience bad curvature. Furthermore, Eq. (43) showsthat there is little freedom to optimize quasi-symmetric stellarators for trapped electronmodes besides varying the value of η . Finally, an estimate for Q bounce for the elevenquasisymmetric designs considered in this study using Eq. (43) is shown in Fig. 5. Wenote that in Fig. 5 there are are no systematic differences between quasi-axisymmetricand quasi-helically symmetric stellarators. Using a quasilinear model, in Ref. [14, 18] the following proxy function Q prox for ITGtransport was derived Q prox = − c D B γk ψ |∇ ψ | ψ TL T , (44) TG stability near the magnetic axis .
00 0 .
25 0 .
50 0 .
75 1 .
00 1 .
25 1 .
50 1 .
75 2 . ηL/ π . . . . . . . A v e r ag e ( Q p r o x Q ) WISTELL-ANZ1988HSXKuQHS48 DrevlakNCSXARIES-CSQAS2ESTELLCFQS Henneberg
Figure 6.
Normalized proxy function Q prox /Q computed using Eq. (46) for each ofthe eleven quasisymmetric designs considered in this study averaged along a field linewith α = 0 . The designs with hollow circles are quasi-axisymmetric, while the oneswith filled circles are quasi-helically symmetric. with c D a constant, B a reference magnetic field and γ the growth rate resulting fromthe following simplified ITG dispersion relation γ (cid:39) ω ∗ L n | τ T κ ( κ p − κ cr ) | / H ( κ p − κ cr ) H ( − κ ) , (45)where κ is the radial component of the curvature vector b · ∇ b , κ p = (1 + η ) /L n isthe pressure gradient length and κ cr a critical pressure gradient length stemming froma contribution of the ITG slab branch. Minimization of this proxy was shown to be aneffective way to reduce turbulent transport levels, raising the prospect of a new class ofstellarators with improved overall confinement [18].Using the expression for |∇ ψ | in Eq. (13) and choosing B = B , we can write theproxy function in Eq. (44) as Q prox = Q √ η cos ϑB H (cos ϑ ) (cid:20) η κ sin ϑ + κ η (cos ϑ − σ sin ϑ ) (cid:21) (46)where Q = − c D k ψ TL T k y ρ th c s a (cid:114) τ T L n (cid:112) η + κ cr L n H (1 + η + κ cr L n ) . (47)The function Q prox in Eq. (46) can then be used as a proxy for the reduction of ITGturbulence near the axis of quasisymmetric stellarators.In Fig. 6, the normalized proxy function Q prox /Q is evaluated for the configurationsconsidered in this study, averaged along a field line with α = 0 . Due to the fact thatthe parameter Q in Eq. (47) is not directly dependent on the geometry used, it isconsidered constant in Fig. 6. The resulting averaged proxy function in Fig. 6 suggests TG stability near the magnetic axis - - R [ meters ] Z [ m e t e r s ] Figure 7.
Top left: top view of a surface of constant toroidal flux ψ for the quasi-helically symmetric NAQS stellarator (Section 5.4 of Ref. [55]) at first order in (cid:15) . Colour indicates B on the outermost flux surface. An example field-line with − π < ϑ < π , α = 0 and r = 0 . is shown in black together with the location of ϑ = 0 . Top right: ten poloidal planes and their respective axis location, together witha field line with − π < ϑ < π and α = 0 . Bottom: side view of the toroidal fluxsurface. that there is not a significant difference in Q between quasi-axisymmetric and quasi-helically symmetric configurations. Indeed, the range of values of Q prox in Fig. 6 canvary two orders of magnitude in each set of quasisymmetric configurations. We notethat this analysis assumes Q is constant across all stellarator designs, which may notbe true in practice due to possible variations of characteristic wave-vectors and criticalgradients between configurations.
4. Numerical Benchmark with a Near-Axis Geometry
Before assessing the linear stability of the eleven quasisymmetric designs and itscomparison with the near-axis expansion, we start with a comparison between a solutionof the near-axis equations and a VMEC configuration based on this solution. Thisallows us to examine the details of the input parameters and respective convergencetests, the comparison between eigenfunctions, the dependence of the growth rate onthe perpendicular wave-vector and the variation of the maximum growth rate with thedensity and temperature gradient scale lengths before presenting the main results in thenext section. For this analysis, we take the quasi-helically symmetric solution of Section
TG stability near the magnetic axis ψ = 0 . T m , B = 1 . T and η = 2 . ,with a resulting rotational transform on-axis of ι = 1 . when solving Eq. (19). Thecomparison between the geometrical quantities is shown in Fig. 8 using flux surfaceswith s = ψ/ψ a (with ψ a the toroidal flux at the plasma boundary) ranging from 0.01 to0.8 and a field line with α = 0 . Only six of the eight quantities are shown as Eqs. (11)and (12) and Eqs. (16) and (18) are independent in the limit of near-axis quasisymmetricfields used here. Although, overall, there is good agreement, some differences betweenthe two methods arise, even at s = 0 . . The fact that a magnetic field solution obtainedfrom setting the boundary value of B from a first-order near-axis solution at a particularradius r is not perfectly quasisymmetric is a phenomenon that has been recently beenobserved in Ref. [56] and addressed using second-order solutions in Ref. [55].We proceed by using the GS2 code to assess the properties of ITG modes in boththe VMEC and the near-axis solution. For each simulation, the temporal evolution of | ˆ φ | is fitted to an exponential of the form e γt with γ the growth rate. As a base casescenario, we use ngauss =3 (half of the number of untrapped pitch-angles moving inone direction along field line), deltat =0.4 (the temporal increment between time steps), nstep=150 (the total number of time steps), ngrid =10 (the total number of energypoints), nzgrid =100 (the total number of points in z along the field line), nlambda =24(number of trapped-pitch angles along the field line) and nfp =4 (number of field periods,a measure of the total length of the field line). In Fig. 9, a convergence test is performedwhere the base case growth rate for both the VMEC and the near-axis geometries iscomputed, along with the growth rate resulting from doubling the resolution along eachparameter. As physical input parameters, we use k y ρ = 1 . , a/L T = 3 and a/L n = 1 .A base case scenario is considered to be converged if the scatter plot cluster of allsimulations has a deviation smaller than 5%, a criterion that is satisfied for all radii inFig. 9.For the base case scenario at s = 0 . , we show in Fig. 10 the resulting real andimaginary electrostatic potential eigenfunctions ˆ φ at the end of the simulation periodfor both the VMEC (left) and near-axis (right) geometries. Overall, both eigenfunctionslook practically indistinguishable from each other for every simulation of the convergencetest. Finally, we note that these are smooth functions of the parallel coordinate z and reach insignificant values before the end of the simulation domain, which providesadditional evidence for the identification of a converged simulation scenario.Next, we perform a scan over k y in order to determine the peak growth rate γ foreach value of L n and L T . For the values of a/L T = 5 and a/L n = 2 considered above, weshow in Fig. 11 the dependence of γ = γ ( k y ) and ω = ω ( k y ) for ≤ k y ρ ≤ for bothconfigurations at s = 0 . . The peak growth rate is found at k y ρ = 2 , a value within theexpected range of the ITG instability. Also, the near-axis expansion appears to yield alarger value of the growth rate across the majority of the spectrum, with a particularincidence at the peak values of k y ∼ .Finally, we vary the density and temperature gradient scale lengths, a/L n and a/L T , TG stability near the magnetic axis NAQS - - ϕ - - ϕ - - ϕ - - ϕ - - - - - ϕ - - - - - ϕ - - Near - Axis s = = = = = = = = = = = = Figure 8.
Comparison between the fitted and the original near-axis geometrycoefficients for the NAQS configuration of Section 5.4 of Ref [56] along the cylindricaltoroidal angle φ at several radial locations s = ψ/ψ a with ψ a the toroidal flux at theplasma boundary. All quantities are expressed in SI units. TG stability near the magnetic axis . . . Near-Axis γ . . . V M E C γ NAQS s=0.01 base casedouble ngausshalf dtdouble nstepdouble negriddouble nzgriddouble nlambdadouble nfp . . . Near-Axis γ . . . V M E C γ NAQS s=0.05 base casedouble ngausshalf dtdouble nstepdouble negriddouble nzgriddouble nlambdadouble nfp . . . Near-Axis γ . . . V M E C γ NAQS s=0.1 base casedouble ngausshalf dtdouble nstepdouble negriddouble nzgriddouble nlambdadouble nfp . . . Near-Axis γ . . . V M E C γ NAQS s=0.3 base casedouble ngausshalf dtdouble nstepdouble negriddouble nzgriddouble nlambdadouble nfp . . . Near-Axis γ . . . V M E C γ NAQS s=0.5 base casedouble ngausshalf dtdouble nstepdouble negriddouble nzgriddouble nlambdadouble nfp . . . Near-Axis γ . . . V M E C γ NAQS s=0.8 base casedouble ngausshalf dtdouble nstepdouble negriddouble nzgriddouble nlambdadouble nfp
Figure 9.
Growth rate convergence test of the base case scenario for the VMEC andnear-axis NAQS geometries using the gyrokinetic GS2 code at different radii s = ψ/ψ a and k y ρ = 1 , a/L T = 3 and a/L n = 1 . − − z . . . . . ˆ φ / ˆ φ Re( ˆ φ/ ˆ φ )Im( ˆ φ/ ˆ φ ) − − z . . . . . ˆ φ / ˆ φ Re( ˆ φ/ ˆ φ )Im( ˆ φ/ ˆ φ ) Figure 10.
Real and imaginary eigenfunctions ˆ φ at the end of the simulation period forboth the VMEC (left) and near-axis (right) geometries at s = 0 . . Each eigenfunctionis normalized to its value at z = 0 . and pick the fastest growing mode by selecting the maximum value of γ ( k y ) for each L n and L T . The resulting peak growth rates, associated ω and k y within ≤ a/L n ≤ and ≤ a/L T ≤ are shown in Fig. 12. The locations of ω and k y where γ ≤ are shownin white. As is characteristic of the ITG mode, the peak growth rate is negative for arange of values where η = L n /L T > η crit with η crit ∼ [57], such as the condition forinstability η > / of the toroidal branch of the ITG mode found in [58]. A more detailedunderstanding of the mechanisms that set the ITG critical gradient using first-principleapproaches, such as the one in Ref. [59], are left for future work. TG stability near the magnetic axis − − k y . . . . . . . . γ Near-AxisVMEC − − k y . . . . . . ω Near-AxisVMEC
Figure 11.
Dependence of γ = γ ( k y ) and ω = ω ( k y ) for ≤ k y ρ ≤ for the VMECand near-axis NAQS configuration for a/L T = 5 and a/L n = 2 at s = 0 . .
5. Numerical results
In this section, we carry out the analysis performed in Section 4 for a single benchmarkcase, now for the eleven realistic quasisymmetric designs under study. In particular,we perform a convergence scan to determine a base case scenario that satisfies theconvergence criterion for all stellarator designs, perform a scan over k y , L n and L T andpick the fastest growth rate in order to assess the critical gradients η crit and peak growthrates γ max .We start with the convergence scan. A comparison between the growth rateobtained with the near-axis expansion and VMEC is shown in Fig. 13 for several radii s = ψ/ψ a = (0 . , . , . , . , . , . and k y ρ = 1 , a/L T = 3 and a/L n = 1 . Inorder to obtain converged base case scenarios, the spatial resolution was increasedto nzgrid =150, nlambda increased to nlambda =28 and, due to the presence of lowervalues of γ with respect to the NAQS benchmark, the simulation time was increased to nstep=250 . Concerning the top row in Fig. 13 (radii with s < . ), the growth rate γ obtained with a near-axis model is found to be in very good agreement with the onesusing the corresponding VMEC design. However, it is found that for larger radii, thenear-axis expansion provides an overestimation of the growth rate. This behaviour isobserved for all cases studied here.We distinguish between the quasi-axisymmetric and quasi-helically symmetricconfigurations in Fig. 13 by using both hollow and filled circles, respectively. Itis seen that, in general, there is no systematic difference between the two for thephysical parameters simulated here, except for a small tendency of quasi-axisymmetricstellarators close to the magnetic axis to exhibit higher growth rates. However, thistendency seems to reverse further away from the axis where, at s = 0 . , the majority ofthe quasi-axisymmetric configurations have considerably lower growth rates than quasi- TG stability near the magnetic axis a/L n a / L T VMEC . . . . . γ a/L n a / L T Near-Axis . . . . . γ a/L n a / L T VMEC . . . . . . . ω a/L n a / L T Near-Axis . . . . . . . ω a/L n a / L T VMEC k y a/L n a / L T Near-Axis k y Figure 12.
Comparison of the peak growth rates γ and associated real frequencies ω and wave-numbers k y between the VMEC and near-axis geometries for the NAQSstellarator at s = 0 . . The gradient lengths are scanned over ≤ a/L n ≤ and ≤ a/L T ≤ . TG stability near the magnetic axis . . . Near-Axis γ . . . V M E C γ WISTELL-ANZ1988 HSXKuQHS48DrevlakNCSXARIES-CS QAS2ESTELLCFQS Henneberg s=0.01 base casedouble ngausshalf dtdouble nstepdouble negriddouble nzgriddouble nlambdadouble nfp . . . Near-Axis γ . . . V M E C γ WISTELL-ANZ1988 HSXKuQHS48DrevlakNCSXARIES-CSQAS2ESTELLCFQSHenneberg s=0.05 base casedouble ngausshalf dtdouble nstepdouble negriddouble nzgriddouble nlambdadouble nfp . . . Near-Axis γ . . . V M E C γ WISTELL-ANZ1988 HSXKuQHS48DrevlakNCSXARIES-CSQAS2ESTELLCFQS Henneberg s=0.1 base casedouble ngausshalf dtdouble nstepdouble negriddouble nzgriddouble nlambdadouble nfp . . . Near-Axis γ . . . V M E C γ WISTELL-ANZ1988 HSXKuQHS48DrevlakNCSXARIES-CSQAS2ESTELL CFQSHenneberg s=0.3 base casedouble ngausshalf dtdouble nstepdouble negriddouble nzgriddouble nlambdadouble nfp . . . Near-Axis γ . . . V M E C γ WISTELL-ANZ1988HSXKuQHS48DrevlakNCSXARIES-CSQAS2 ESTELLCFQSHenneberg s=0.5 base casedouble ngausshalf dtdouble nstepdouble negriddouble nzgriddouble nlambdadouble nfp . . . Near-Axis γ . . . V M E C γ WISTELL-ANZ1988HSXKuQHS48DrevlakNCSXARIES-CSQAS2 ESTELLCFQSHenneberg s=0.8 base casedouble ngausshalf dtdouble nstepdouble negriddouble nzgriddouble nlambdadouble nfp
Figure 13.
Growth rate convergence test of the base case scenario for the VMEC andnear-axis geometries using the gyrokinetic GS2 code at different radii s = ψ/ψ a and k y ρ = 1 , a/L T = 3 and a/L n = 1 . Quasi-axisymmetric stellarators are representedwith hollow circles while quasi-helically symmetric ones are represented with filledcircles. helically symmetric ones, widening the gap between the maximum and minimum valuesof the computed γ .We now comment on the accuracy of the estimate using a simplified dispersionrelation in Eq. (36) when compared with the growth rates using first-principlegyrokinetic simulations in Fig. 13. This comparison is shown in Fig. 14, where we showthe growth rate estimate γ normalized to ω t computed using the dispersion relationin Eq. (36) (denoted as γ DR ) on the horizontal axis and the growth rate estimate γ normalized to ω t computed using the gyrokinetic GS2 code and the correspondingVMEC equilibrium files (denoted as γ GS , corresponding to the values in Fig. 13multiplied by L/ πa ). It is seen that the estimate γ DR resulting from the simplifieddispersion relation in Eq. (36) consistently overestimates the growth rate which, insome cases, can lead to a one order of magnitude difference between the two growthrates for quasi-axisymmetric stellarators with higher values of γ . However, the generaltrend and the relative growth rates between different configurations is well predicted bythe simplified dispersion relation.Finally, we assess the accuracy of the near-axis expansion at predicting the peakgrowth rates when performing a scan over a/L n and a/L T . Two examples are shown inFig. 15 and Fig. 16 for the NCSX and HSX stellarators, respectively. These confirm the TG stability near the magnetic axis . . . . . . ˜ γ DR . . . . . . ˜ γ G S WISTELL-ANZ1988 HSXKuQHS48Drevlak NCSXARIES-CS QAS2ESTELLCFQS Henneberg
Figure 14.
Horizontal axis: growth rate estimate γ normalized to ω t computedusing the dispersion relation in Eq. (36) (denoted as γ DR ). Vertical axis: growthrate estimate γ normalized to ω t computed using the gyrokinetic GS2 code and thecorresponding VMEC equilibrium files (denoted as γ GS , corresponding to the valuesin Fig. 13 multiplied by L/ πa ). Quasi-axisymmetric stellarators are represented withhollow circles while quasi-helically symmetric ones are represented with filled circlesand r = 0 . . general trend observed in Fig. 12 and the accuracy of the near-axis expansion for valuesof s < . . In general, the linear critical gradient for the ITG instability near the axis ofquasisymmetric stellarators lies at values of η ∼ . We also find that the peak growthrate increases with increasing a/L T and is maximum for values of < a/L n < . Thevalues of the real frequency ω at the peak growth rate increase with increasing a/L T while peaking at either a/L n = 0 or a/L n (cid:39) . The wave-vector k y at the peak growthrate is observed to increase mainly with increasing a/L n .
6. Conclusions
In this work, the ion-temperature-gradient stability of quasisymmetric stellarators isassessed near the magnetic axis by employing the gyrokinetic approximation togetherwith a near-axis expansion framework. It is found that the near-axis expansion is notonly able to simplify analytical estimates of ITG growth rates and turbulent fluxes, butalso able to reproduce the properties of eleven quasisymmetric configurations designedusing numerical optimization methods. In general, as the radius is increased beyondthe inner core, it is found that the near-axis model overestimates the growth rateswhen compared with the corresponding optimized VMEC design. Such overestimationcould arise due to the absence of global magnetic shear ˆ s = − ( r/ι ) ι (cid:48) ( r ) in the near-axis expansion formalism at the order considered here which could contribute as astabilizing factor at larger radii where higher order geometry effects are more relevant.Global magnetic shear may also provide additional finite-Larmor radius stabilization TG stability near the magnetic axis a/L n a / L T VMEC . . . . . . γ a/L n a / L T Near-Axis . . . . . . γ a/L n a / L T VMEC . . . . . . . ω a/L n a / L T Near-Axis . . . . . . . ω a/L n a / L T VMEC k y a/L n a / L T Near-Axis k y Figure 15.
Comparison of the peak growth rates γ and associated real frequencies ω and wave-numbers k y between the VMEC and near-axis geometries for the NCSXstellarator at s = 0 . . The gradient lengths are scanned over ≤ a/L n ≤ and ≤ a/L T ≤ . TG stability near the magnetic axis a/L n a / L T VMEC . . . . . γ a/L n a / L T Near-Axis . . . . . γ a/L n a / L T VMEC . . . . . . . . ω a/L n a / L T Near-Axis . . . . . . . . ω a/L n a / L T VMEC k y a/L n a / L T Near-Axis k y Figure 16.
Comparison of the peak growth rates γ and associated real frequencies ω and wave-numbers k y between the VMEC and near-axis geometries for the HSXstellarator at s = 0 . . The gradient lengths are scanned over ≤ a/L n ≤ and ≤ a/L T ≤ . TG stability near the magnetic axis k ⊥ . Furthermore, the presence of inhomogeneities in the magnetic field of the VMECconfigurations, which could lead to departures from quasisymmetry, may allow thepresence of a larger fraction of trapped particles at the good curvature side and,therefore, decrease γ .While the analysis performed here pertained the ITG instability only, a similaranalysis can be carried out to other instabilities, such as the trapped-electron mode,the electron-temperature-gradient mode and the micro-tearing mode. Furthermore,as future work, we intend to generalize the numerical methods employed herefor quasisymmetric designs to maximum-J devices [60] and have a more completeunderstanding of the supression of electrostatic microinstabilities recently found in suchconfigurations [61].
7. Acknowledgments
We wish to thank I. G. Abel, W. Dorland, G. T. Roberg-Clark and G. G. Plunk for manyfruitful discussions regarding the content of the manuscript. This work was supportedby a grant from the Simons Foundation (560651, ML) and from the U.S. Departmentof Energy, Office of Science, Office of Fusion Energy Science, under award number DE-FG02-93ER54197.
References [1] P. Xanthopoulos and F. Jenko. Gyrokinetic analysis of linear microinstabilities for the stellaratorWendelstein 7-X.
Physics of Plasmas , 14(4):042501, 2007.[2] P. Helander, T. Bird, F. Jenko, R. Kleiber, G. G. Plunk, J. H.E. Proll, J. Riemann, andP. Xanthopoulos. Advances in stellarator gyrokinetics.
Nuclear Fusion , 55(5):053030, 2015.[3] J. Nuhrenberg and R. Zille. Quasi-helically symmetric toroidal stellarators.
Physics Letters A ,129(2):113–117, 1988.[4] A. H. Boozer. Quasi-helical symmetry in stellarators.
Plasma Physics and Controlled Fusion ,37(11A):A103, 1995.[5] P. R. Garabedian. Stellarators with the magnetic symmetry of a tokamak.
Physics of Plasmas ,3(7):2483, 1996.[6] D. A. Garren and A. H. Boozer. Existence of quasihelically symmetric stellarators.
Physics ofFluids B , 3(10):2822, 1991.[7] M. Landreman and W. Sengupta. Direct construction of optimized stellarator shapes. Part 1.Theory in cylindrical coordinates.
Journal of Plasma Physics , 84(6):905840616, 2018.[8] R. Jorge, W. Sengupta, and M. Landreman. Construction of quasisymmetric stellarators using adirect coordinate approach.
Nuclear Fusion , 60(7):076021, 2020.[9] G. G. Plunk. Perturbing an axisymmetric magnetic equilibrium to obtain a quasi-axisymmetricstellarator.
Journal of Plasma Physics , 86(4):905860409, 2020.[10] B. E. Nelson, L. A. Berry, A. B. Brooks, M. J. Cole, J. C. Chrzanowski, H. M. Fan, P. J. Fogarty,P. L. Goranson, P. J. Heitzenroeder, S. P. Hirshman, G. H. Jones, J. F. Lyon, G. H. Neilson,W. T. Reiersen, D. J. Strickler, and D. E. Williamson. Design of the national compact stellaratorexperiment (NCSX). In
Fusion Engineering and Design , volume 66-68, page 169, 2003.[11] J. M. Canik, D. T. Anderson, F. S.B. Anderson, K. M. Likin, J. N. Talmadge, and K. Zhai.
TG stability near the magnetic axis Experimental demonstration of improved neoclassical transport with quasihelical symmetry.
Physical Review Letters , 98(8):085002, 2007.[12] M. Drevlak, C. D. Beidler, J. Geiger, P. Helander, and Y. Turkin. Optimisation of stellaratorequilibria with ROSE.
Nuclear Fusion , 59(1):016010, 2019.[13] P. Helander, J. H.E. Proll, and G. G. Plunk. Collisionless microinstabilities in stellarators. I.Analytical theory of trapped-particle modes.
Physics of Plasmas , 20(12):122505, 2013.[14] H. E. Mynick, N. Pomphrey, and P. Xanthopoulos. Reducing turbulent transport in toroidalconfigurations via shaping.
Physics of Plasmas , 18(5):056101, 2011.[15] G. Rewoldt, W. M. Tang, and M. S. Chance. Electromagnetic kinetic toroidal eigenmodes forgeneral magnetohydrodynamic equilibria.
Physics of Fluids , 25(3):480, 1982.[16] G. Rewoldt, W. M. Tang, and R. J. Hastie. Collisional effects on kinetic electromagnetic modesand associated quasilinear transport.
Physics of Fluids , 30(3):807, 1987.[17] P. Xanthopoulos, F. Merz, T. Görler, and F. Jenko. Nonlinear gyrokinetic simulations of ion-temperature-gradient turbulence for the optimized wendelstein 7-X stellarator.
Physical ReviewLetters , 99(3):035002, 2007.[18] H. E. Mynick, N. Pomphrey, and P. Xanthopoulos. Optimizing stellarators for turbulent transport.
Physical Review Letters , 105(9):095004, 2010.[19] W. Dorland, F. Jenko, M. Kotschenreuther, and B. N. Rogers. Electron Temperature GradientTurbulence.
Physical Review Letters , 85(26):5579, 2000.[20] A. Zocco, P. Xanthopoulos, H. Doerk, J. W. Connor, and P. Helander. Threshold for thedestabilisation of the ion-temperature-gradient mode in magnetically confined toroidal plasmas.
Journal of Plasma Physics , 84(1):715840101, 2018.[21] P. Xanthopoulos, G. G. Plunk, A. Zocco, and P. Helander. Intrinsic turbulence stabilization in astellarator.
Physical Review X , 6(2):021033, 2016.[22] C. C. Hegna, P. W. Terry, and B. J. Faber. Theory of ITG turbulent saturation in stellarators:Identifying mechanisms to reduce turbulent transport.
Physics of Plasmas , 25(2):022511, 2018.[23] I. J. McKinney, M. J. Pueschel, B. J. Faber, C. C. Hegna, J. N. Talmadge, D. T. Anderson, H. E.Mynick, and P. Xanthopoulos. A comparison of turbulent transport in a quasi-helical and aquasi-axisymmetric stellarator.
Journal of Plasma Physics , 85(5):905850503, 2019.[24] H. Y. Wang, I. Holod, Z. Lin, J. Bao, J. Y. Fu, P. F. Liu, J. H. Nicolau, D. Spong, and Y. Xiao.Global gyrokinetic particle simulations of microturbulence in W7-X and LHD stellarators.
Physics of Plasmas , 27(8):082305, 2020.[25] M. Landreman. Optimized quasisymmetric stellarators are consistent with the Garren–Boozerconstruction.
Plasma Physics and Controlled Fusion , 61(7):075001, 2019.[26] R. Jorge and M. Landreman. The Use of Near-Axis Magnetic Fields for Stellarator TurbulenceSimulations.
Plasma Physics and Controlled Fusion , 63(1):014001, 2020.[27] B. J. Faber, M. J. Pueschel, P. W. Terry, C. C. Hegna, and J. E. Roman. Stellaratormicroinstabilities and turbulence at low magnetic shear.
Journal of Plasma Physics ,84(5):905840503, 2018.[28] M. F. Martin, M. Landreman, P. Xanthopoulos, N. R. Mandell, and W. Dorland. The parallelboundary condition for turbulence simulations in low magnetic shear devices.
Plasma Physicsand Controlled Fusion , 60(9):095008, 2018.[29] J. A. Baumgaertel, E. A. Belli, W. Dorland, W. Guttenfelder, G. W. Hammett, D. R. Mikkelsen,G. Rewoldt, W. M. Tang, and P. Xanthopoulos. Simulating gyrokinetic microinstabilities instellarator geometry with GS2.
Physics of Plasmas , 18(12):122301, 2011.[30] G. G. Howes, S. C. Cowley, W. Dorland, G. W. Hammett, E. Quataert, and A. A. Schekochihin.Astrophysical Gyrokinetics: Basic Equations and Linear Theory.
The Astrophysical Journal ,651(1):590, 2006.[31] E. A. Frieman and L. Chen. Nonlinear gyrokinetic equations for low-frequency electromagneticwaves in general plasma equilibria.
Physics of Fluids , 25(3):502, 1982.[32] S. P. Hirshman and J. C. Whitson. Steepest-descent moment method for three-dimensional
TG stability near the magnetic axis magnetohydrodynamic equilibria. Physics of Fluids , 26(12):3553, 1983.[33] A. H. Boozer. Plasma equilibrium with rational magnetic surfaces.
Physics of Fluids , 24(11):1999,1981.[34] D. A. Garren and A. H. Boozer. Magnetic field strength of toroidal plasma equilibria.
Physics ofFluids B , 3(10):2805, 1991.[35] C. Mercier. Equilibrium and stability of a toroidal magnetohydrodynamic system in theneighbourhood of a magnetic axis.
Nuclear Fusion , 4(3):213, 1964.[36] L. S. Solov’ev and V. D. Shafranov.
Reviews of Plasma Physics 5 . Consultants Bureau, New York- London, 1970.[37] R. Jorge, W. Sengupta, and M. Landreman. Near-axis expansion of stellarator equilibrium atarbitrary order in the distance to the axis.
Journal of Plasma Physics , 86(1):905860106, 2020.[38] J. Nuhrenberg and R. Zille. Quasi-helically symmetric toroidal stellarators.
Physics Letters A ,129(2):113, 1988.[39] M. Drevlak. Stellarator configuration design using ROSE. In , page 2, 2017.[40] F. S. B. Anderson, A. F. Almagri, D. T. Anderson, P. G. Matthews, J. N. Talmadge, andJ. L. Shohet. The Helically Symmetric Experiment, (HSX) Goals, Design and Status.
FusionTechnology , 27(3T):273, 1995.[41] L. P. Ku and A. H. Boozer. New classes of quasi-helically symmetric stellarators.
Nuclear Fusion ,51(1):013004, 2011.[42] Aaron Bader. Dataset for Wistell-A stellarator: zenodo.org/record/3952862
Plasma Physics and Controlled Fusion , 43(12A):A237, 2001.[44] F. Najmabadi, A. R. Raffray, S. I. Abdel-Khalik, L. Bromberg, L. Crosatti, L. El-Guebaly, P. R.Garabedian, A. A. Grossman, D. Henderson, A. Ibrahim, T. Ihli, T. B. Kaiser, B. Kiedrowski,L. P. Ku, J. F. Lyon, R. Maingi, S. Malang, C. Martin, T. K. Mau, B. Merrill, R. L. Moore, R. J.Peipert, D. A. Petti, D. L. Sadowski, M. Sawan, J. H. Schultz, R. Slaybaugh, K. T. Slattery,G. Sviatoslavsky, A. Turnbull, L. M. Waganer, X. R. Wang, J. B. Weathers, P. Wilson, J. C.Waldrop, M. Yoda, and M. Zarnstorff. The ARIES-CS compact stellarator fusion power plant.
Fusion Science and Technology , 54(3):655, 2008.[45] P. R. Garabedian. Three-dimensional analysis of tokamaks and stellarators.
Proceedings of theNational Academy of Sciences of the United States of America , 105(37):13716, 2008.[46] M. Drevlak, F. Brochard, P. Helander, J. Kisslinger, M. Mikhailov, C. Nührenberg, J. Nührenberg,and Y. Turkin. ESTELL: A Quasi-Toroidally Symmetric Stellarator.
Contributions to PlasmaPhysics , 53(6):459, 2013.[47] A. Shimizu, H. Liu, M. Isobe, S. Okamura, S. Nishimura, C. Suzuki, Y. Xu, X. Zhang, B. Liu,J. Huang, X. Wang, H. Liu, and C. Tang. Configuration property of the Chinese first quasi-axisymmetric stellarator.
Plasma and Fusion Research , 13:3403123, 2018.[48] S. A. Henneberg, M. Drevlak, C. Nührenberg, C. D. Beidler, Y. Turkin, J. Loizu, and P. Helander.Properties of a new quasi-axisymmetric configuration.
Nuclear Fusion , 59(2):026014, 2019.[49] M. Barnes, F. I. Parra, and M. Landreman. stella: An operator-split, implicit–explicit δ f-gyrokinetic code for general magnetic field configurations. Journal of Computational Physics ,391:365, 2019.[50] J. H.E. Proll, P. Xanthopoulos, and P. Helander. Collisionless microinstabilities in stellarators. II.Numerical simulations.
Physics of Plasmas , 20(12):122506, 2013.[51] G. G. Plunk, P. Helander, P. Xanthopoulos, and J. W. Connor. Collisionless microinstabilities in
TG stability near the magnetic axis stellarators. III. the ion-temperature-gradient mode. Physics of Plasmas , 21(3):032112, 2014.[52] A. Zocco, G. G. Plunk, and P. Xanthopoulos. Geometric stabilization of the electrostatic ion-temperature-gradient driven instability. II. Non-axisymmetric systems.
Physics of Plasmas ,27(2):022507, 2020.[53] J. H.E. Proll, H. E. Mynick, P. Xanthopoulos, S. A. Lazerson, and B. J. Faber. TEM turbulenceoptimisation in stellarators.
Plasma Physics and Controlled Fusion , 58(1):014006, 2015.[54] O. Beeke, M. Barnes, M. Romanelli, M. Nakata, and M. Yoshida. Impact of plasma shaping ontokamak microstability. arXiv:2012.02669 , 2020.[55] M. Landreman and W. Sengupta. Constructing stellarators with quasisymmetry to high order.
Journal of Plasma Physics , 85(6):815850601, 2019.[56] M. Landreman, W. Sengupta, and G. G. Plunk. Direct construction of optimizedstellarator shapes. Part 2. Numerical quasisymmetric solutions.
Journal of Plasma Physics ,85(1):905850103, 2019.[57] F. Romanelli. Ion temperature-gradient-driven modes and anomalous ion transport in tokamaks.
Physics of Fluids B , 1(5):1018–1025, 1989.[58] H. Biglari, P. H. Diamond, and M. N. Rosenbluth. Toroidal ion-pressure-gradient-driven driftinstabilities and transport revisited.
Physics of Fluids B , 1(1):109, 1989.[59] G. T. Roberg-Clark, G. G. Plunk, and P. Xanthopoulos. Calculating the linear critical gradient forthe ion-temperature-gradient mode in magnetically confined plasmas. arXiv:2010.13441 , 2020.[60] M. N. Rosenbluth. Low-frequency limit of interchange instability.
Physics of Fluids , 11(4):869,1968.[61] J. A. Alcuson, P. Xanthopoulos, G. G. Plunk, P. Helander, F. Wilms, Y. Turkin, A. Von Stechow,and O. Grulke. Suppression of electrostatic micro-instabilities in maximum-J stellarators.