Optimization of quasi-symmetric stellarators with self-consistent bootstrap current and energetic particle confinement

Stellarators have many favorable attributes as a fusion reactor concept, including the absence of disruptions, no need for current drive and the associated recirculating power, intrinsic steady-state capability, and the ability to reliably extrapolate to future devices using computational design.¹ Stellarators with quasi-symmetry, in which the magnetic field strength has an invariant direction in Boozer coordinates, are of particular interest, as quasi-symmetry enables the confinement of fusion-produced alpha particles and reduction of neoclassical transport to acceptable levels.² At reactor-relevant densities and temperatures, a quasi-symmetric stellarator will have a significant bootstrap current, a current parallel to the magnetic field driven by gradients in the density and temperature. The bootstrap current is a kinetic effect, arising from the different guiding-center trajectories of electrons and ions as well as collisions. The current in the plasma will produce a contribution to the magnetic field that must be added to the field of the external electromagnetic coils. Accurate reckoning of the bootstrap current is, therefore, necessary to know the total magnetic field and, hence, to know almost any property of interest.

In this article, a new efficient method is described for including the bootstrap current in the optimization of quasi-symmetric stellarators. The approach is a combination of three ideas. The first is a new objective function and associated choice of a parameter space for obtaining self-consistency of the current profile. The second is the application of the isomorphism between quasi-symmetric stellarators and tokamaks,^3–5 allowing the bootstrap current to be computed as if the geometry was axisymmetric. The third is the application of a recent analytic formula for the tokamak bootstrap current by Redl et al.⁶ These last two components greatly reduce the computational cost of stellarator bootstrap current calculations and, as we will show, give quite accurate results compared to bootstrap calculations that fully account for three-dimensional geometry. As examples of the new optimization approach, we demonstrate new quasi-symmetric stellarator configurations with self-consistent current profiles, also showing excellent alpha particle confinement and low neoclassical transport.

The bootstrap current on each flux surface depends on the magnetic geometry. A self-consistent state must be determined, a simultaneous solution of the equations of magnetohydrodynamic (MHD) equilibrium and the drift-kinetic equation. In the MHD equilibrium problem, a profile of either toroidal current or rotational transform (along with a pressure profile) is provided as an input, and the magnetic field B is an output. Conversely, in solutions of the drift-kinetic equation, the magnetic field components are taken as input, and the parallel current on each flux surface is an output. Generally, if an equilibrium calculation is followed by a drift-kinetic calculation, there is no reason for the resulting parallel current profile to match the current profile of the MHD equilibrium. However, the drift-kinetic profile of parallel current can be used as input for a new equilibrium calculation, as the basis for a fixed-point iteration. After several iterations between the equilibrium and drift-kinetic calculations, a self-consistent state is obtained.⁷ While this fixed-point iteration works well for a single configuration, it is not clear how to include this iteration within optimization. If one or a few fixed-point iterations were performed for each evaluation of the objective function, the objective would take on slightly different values each time it is evaluated at the same point in the parameter space, potentially confusing the optimization algorithm, and rendering finite difference derivatives very inaccurate. Here, we circumvent these issues by using a different approach in which there is no fixed-point iteration. Instead, consistency of the MHD equilibrium and the drift-kinetic equation is enforced through a penalty in the objective.

Another complication of including bootstrap calculations inside stellarator optimization is their computational cost. Accurate stellarator bootstrap calculations require the solution of the drift-kinetic equation, a high-dimensional advection-dominated integrodifferential equation, which is numerically challenging. For the case of quasi-symmetric stellarators, we can exploit their isomorphism with axisymmetric geometries, identified by Pytte and Boozer.³ In this isomorphism, formulas for neoclassical phenomena in a quasi-symmetric stellarator are the same as the corresponding formulas in axisymmetry, up to a few substitutions. The primary substitution of interest is that the rotational transform ι in axisymmetry is replaced with $ι−N$⁠, where N = 0 for quasi-axisymmetry and N is the number of field periods for quasi-helical symmetry. As a result, the bootstrap current in a quasi-symmetric stellarator can be computed as if the geometry was axisymmetric, at lower computational cost due to the ignorable coordinate.

Moreover, for the bootstrap current in axisymmetry, accurate analytic formulas are available, eliminating the need for any numerical kinetic calculation. Here, we will use a recent formula for the tokamak bootstrap current by Redl et al.,⁶ which improves on an earlier formula by Sauter et al.⁸ The form of this formula is motivated by analytic limits in collisionality and aspect ratio, with coefficients determined by fits to a database of first-principles numerical solutions of the drift-kinetic equation. The resulting formula is applicable at any collisionality and aspect ratio. Here, we show for the first time that these analytic tokamak expressions are accurate also in quasi-symmetric stellarators.

Expressions for the bootstrap current in a general stellarator are available in the limit of long mean free path compared to system size,^9–11 which can be evaluated faster than the general-collisionality drift-kinetic equation. However, to date, these long-mean-free-path expressions tend to produce noisy results,¹² as they depend on the solutions to magnetic differential equations that are singular on any rational surface. Therefore, ad hoc smoothing of the current profile is sometimes applied. This noisiness is unfavorable for optimization, as smooth objective functions are generally easier to optimize. Also, as will be shown below, the long-mean-free-path formula can give significantly different results in practice to the general-collisionality drift-kinetic equation, the latter being more general and accurate. If improved analytic expressions for the bootstrap current in a general stellarator in the long-mean-free-path regime become available, they could be incorporated into stellarator optimization using the same objective function and parameter space described in this paper.

In Secs. II and III, the quasi-symmetry isomorphism is reviewed, and its application to the bootstrap current formula of Redl et al.⁶ is described. In Sec. IV, the new method of calculating the bootstrap current is demonstrated in existing stellarator configurations. Details of the new optimization procedure are given in Sec. V. Optimization results are presented in Sec. VI, at several values of $β=2μ0p/B2$⁠, where p is the plasma pressure and $B=|B|$⁠. Advantages of the isomorphism approach over long-mean-free-path expressions for the stellarator bootstrap current are illustrated in Sec. VII. In Sec. VIII, the excellent confinement of these configurations is demonstrated, and we conclude in Sec. IX.

In this section, we review the isomorphism between quasi-symmetric and axisymmetric configurations,^3–5 focusing on the bootstrap current. For this discussion, we consider the magnetic field represented in terms of the Boozer poloidal and toroidal angles θ and $φ$⁠. We introduce a helical angle $χ=θ−Nφ$ for some integer N, where N is arbitrary for the moment. Quasi-symmetry can be expressed as the condition $B=B(ψ,χ)$⁠, so B depends on θ and $φ$ only through a fixed linear combination of the angles. We will consider two types of quasi-symmetry: quasi-axisymmetry (QA) and quasi-helical (QH) symmetry. QA is the case in which N = 0, so $B=B(ψ,θ)$⁠, while QH symmetry is the case of nonzero N, $B=B(ψ,θ−Nφ)$⁠. Note that axisymmetry is a special case of QA.

Now, without assuming quasi-symmetry, we proceed to determine the bootstrap current in a general stellarator or tokamak. The calculation begins by writing the Boozer coordinate representation of the field

and

Here, ψ is the toroidal flux divided by $2π$⁠, ι is the rotational transform, $ι̃=ι−N$⁠, $G(ψ)$ is $μ0/(2π)$ times the poloidal current outside the flux surface ψ, $I(ψ)$ is $μ0/(2π)$ times the toroidal current inside the flux surface ψ, and $G̃(ψ)=G+NI$⁠. The Jacobian of the $(ψ,χ,φ)$ coordinates is $g=(∇ψ·∇χ×∇φ)−1=(G̃+ι̃I)/B2=(G+ιI)/B2$⁠.

We use the conventional definition of the bootstrap current as $⟨j·B⟩$⁠, where $j=μ0−1∇×B$ is the current density and $⟨…⟩$ is a flux surface average. The flux surface average of any quantity Q is $⟨Q⟩=(1/V′)∫02πdχ∫02πdφ gQ$⁠, where $V′=∫02πdχ∫02πdφ g$⁠. The bootstrap current can be computed from

where α indicates particle species, $Zα$ is the species charge in units of the proton charge e, $v||$ is the velocity component parallel to B, and $fα$ is the distribution function.

We consider the usual case in which the distribution function is approximately a Maxwellian with no mean flow: $fα≈fα,0=nα[mα/(2πTα)]3/2 exp (−mαv2/(2Tα)),where mα, nα(ψ)$⁠, and $Tα(ψ)$ are the species mass, density, and temperature, and $v=|v|$ is the speed. Deviations from the Maxwellian are then accounted for in higher order terms in the distribution function. We can divide the correction into a term that is independent of gyrophase, $fα,1$⁠, and a gyrophase-dependent term with vanishing gyrophase-average, $f̃α,1$⁠, writing $fα≈fα,0+fα,1+f̃α,1$⁠. Since $f̃α,1$ does not contribute to (3), we will not consider this term further. It is convenient to introduce velocity-space variables $λ=v⊥2/(v2B)$ and $σ=sign(v||)$⁠, where $v⊥=v2−v||2$ is the velocity perpendicular to B, so $fα,1=fα,1(ψ,χ,φ,v,λ,σ)$⁠. The distribution $fα,1$ is found by solving the drift-kinetic equation¹³ 

Here, $∇||=b·∇$⁠, where $b=B−1B$⁠, the gradient is performed at fixed v and λ, and the radial guiding center drift is

(Note that for this radial component, no low-β approximation is required to combine the grad-B and curvature drifts.) In (4), $Ĉα,γ=Cα,γ(fα,1,fγ,0)+Cα,γ(fα,0,fγ,1)$ is the linearized collision operator, with $Cα,γ$ being the nonlinear Fokker–Planck collision operator.^14,15 Once the drift-kinetic equation is solved, the bootstrap current can be computed from (3) in the form

We now introduce the assumption of quasi-symmetry: $B=B(ψ,χ)$⁠, so $(∂B/∂φ)χ=0$⁠. In this case, the drift-kinetic equation (4) can be written

The linearized collision operator can be written in terms of the velocity coordinates v and $v||=σv1−λB$⁠, and so it introduces dependence on position only through B. Similarly, (7) depends on position only through ψ and B, and hence only through ψ and χ. It can be seen then that the entire drift-kinetic equation has position dependence only through $(ψ,χ)$⁠, and so its solution $fα,1$ will have this same property. When the moment (6) of the distribution function is formed, position dependence appears only through B and $g$⁠, which depend on position only through $(ψ,χ)$⁠, with no $φ$-dependence introduced. Thus, the equations for computing the bootstrap current in a quasi-symmetric stellarator, (7) and (6), are identical to the equations one can use to compute the bootstrap current in axisymmetry; in both cases, only two spatial dimensions appear in the equations.

Due to these observations, we reach the following conclusion. The bootstrap current on a quasi-symmetric flux surface of a given nonaxisymmetric field is identical to the bootstrap current on an equivalent axisymmetric flux surface. In this case, “equivalent” means that the quantities appearing in (7) and (6) are matched on that flux surface: B, $ι̃, G̃$⁠, I, $nα, Tα, dnα/dψ$⁠, and $dTα/dψ$⁠. (The equivalent axisymmetric surface need not correspond to any global MHD equilibrium.) Note that N = 0 in the equivalent axisymmetric configuration, even if $N≠0$ in the quasi-symmetric configuration. Thus, we arrive at the following rule: to compute the bootstrap current on a quasi-symmetric flux surface with parameters B, ι, G, and I, we can compute instead the bootstrap current in an axisymmetric surface with corresponding parameters B, $ι−N$⁠, G + NI, and I, also matching $nα, Tα, dnα/dψ$⁠, and $dTα/dψ$⁠.

We will use the recent formula for the tokamak bootstrap current by Redl et al.⁶ This formula is based on fitting to a large database of first-principles drift-kinetic calculations generated by the code NEO¹⁶ (which should not be confused with the stellarator neoclassical code of the same name.) The formula is applicable across all tokamak regimes of collisionality and for any aspect ratio.

The expressions for the bootstrap current by Redl et al. are lengthy and so are relegated in the  Appendix. For application in quasi-symmetric stellarators, only a few minor modifications are made, detailed in the remainder of this section.

First, Redl's Eq. (2) contains the flux function $IR=RBϕ$⁠, where R is the major radius and $Bϕ$ is the toroidal field. [Here, we have added the R subscript to I_R to distinguish it from I in (2).] In axisymmetry, I_R coincides with G. Moreover, the I_R factor in Redl's Eq. (2) originates from the fact that the inhomogeneous drive term in the drift-kinetic equation, the right-hand side of (7), is proportional to $G̃$⁠, which equals I_R = G in axisymmetry. Thus, it is appropriate to replace the I_R factor in Ref. 6 with $G̃$⁠.

Next, in the tokamak expressions of Redl et al.,⁶ the following collisionalities

appear. Here, $q=1/ι$ is the safety factor, R is the major radius in meters, the electron density n_e and ion density n_i are measured in m⁻³, and the electron temperature T_e and ion temperature T_i are measured in eV. In Ref. 6, the major radius $R(ψ)$ of a given flux surface was defined by $(Rmax+Rmin)/2$⁠, where $Rmax(ψ)$ and $Rmin(ψ)$ are the maximum and minimum major radii of the surface. The inverse aspect ratio ϵ was computed as $(Rmax−Rmin)/(2R)$⁠. We must consider how to evaluate the factors $qR/ε3/2$ in a stellarator.

To do so, we consider the physical origin of the expressions (8). The expressions (8) arise from the ratio between the collision operator and the parallel streaming term in the drift-kinetic equation, for the trapped part of phase space. The factor $qR/ε3/2$ arises as a product of factors $1/ε$ and $qR/ε$⁠. The factor $1/ε$ comes from the pitch-angle scattering term in the collision operator: each of the two pitch-angle derivatives is magnified by a factor $1/ε$ for trapped particles, since trapped particles have a parallel velocity $v||∼εv$⁠, and so they only need to diffuse by an angle $∼ε$ to become untrapped. The factor $qR/ε$ comes from the inverse of $v||∇||θ$ in the parallel streaming term, again noting $v||∼εv$⁠. In both the collision and streaming terms, the ϵ factors originate from the size of the trapped region of velocity space, which is determined by $Bmax$ and $Bmin$⁠. In tokamaks, $R∼1/B$ to a good approximation, while in stellarators, it is more accurate to evaluate the size of the trapped region in velocity space using B instead of R. Therefore, a reasonable generalization of ϵ in the Redl formula is

where $Bmax(ψ)$ and $Bmin(ψ)$ are the maximum and minimum B on the flux surface. This definition reduces to that of Redl et al. when $R∝1/B$⁠, as is approximately true in axisymmetry, but is a better reflection of the size of the trapped region of velocity space in a quasi-symmetric stellarator.

The qR factors in (8) originate from the inverse of $∇||θ$ in the streaming term of the drift-kinetic equation. In a quasi-symmetric stellarator, the analogous expression is $1/(b·∇χ)=(G+ιI)/[(ι−N)B]$⁠. Some kind of average must be performed so that the overall expression is a flux function. We choose to evaluate (8) using

Other choices of the average could be reasonable and result in only extremely minor differences.

The magnetic geometry also enters Redl's expressions through the effective fraction of trapped particles

a quantity that arises repeatedly in analytic neoclassical calculations for the tokamak banana regime¹⁵ or stellarator $1/ν$ regime.⁹ This expression requires no modification in a stellarator.

Finally, the density and temperature gradients appear in Ref. 6 as $dnα/dψp$ and $dTα/dψp$⁠, where $2πψp$ is the poloidal flux, and $dψp=ιdψ$⁠. Following the isomorphism, we make the substitutions $dnα/dψp→(ι−N)−1dnα/dψ$ and $dTα/dψp→(ι−N)−1dTα/dψ$⁠.

When the field strength is not perfectly quasi-symmetric, there are several reasonable options available for evaluating f_t and ϵ. One option is to evaluate (11) directly, including integration over the toroidal angle in the flux surface averages, and use the actual maximum and minimum B to evaluate ϵ. A second option is to construct Boozer coordinates, remove any Fourier modes of $B(θ,φ)$ that break the symmetry, and evaluate f_t and ϵ for the resulting perfectly symmetric B pattern. The second option has some extra complexity due to the need to construct Boozer coordinates. However, it is also less sensitive to isolated maxima and minima of B on a flux surface when symmetry is imperfect. Both methods have been compared in the optimization method described in Sec. V, and both are effective. Based on experience so far, the Boozer-coordinate approach results in slightly smaller values of the objective function and so is used for the results that follow.

Unfortunately, the Redl formula is not similarly useful for quasi-poloidally symmetric or quasi-isodynamic stellarators. If these optimizations are achieved perfectly, the quasi-symmetry isomorphism and other analytic calculations^17,18 predict $⟨j·B⟩∝I(ψ)$⁠, and so Eq. (C3) in Ref. 18 gives $⟨j·B⟩=0$⁠. In practice, the properties of quasi-isodynamic or quasi-poloidal symmetry are only achieved imperfectly, and one would prefer to have an estimate of the small nonzero values of $⟨j·B⟩$⁠. Obtaining a fast and smooth calculation of the bootstrap current in imperfectly quasi-poloidally symmetric or quasi-isodynamic stellarators is an important question for future research.

Before embarking on new optimizations, one should verify that the Redl formula with the modifications of Sec. III is, indeed, accurate in stellarators. We now demonstrate this correspondence, using the quasi-axisymmetric (QA) and quasi-helically symmetric (QH) configurations from Ref. 2. For this initial test, no attempt will be made to obtain self-consistency between MHD equilibrium and the drift-kinetic equation. Instead, the bootstrap current will be computed in fixed magnetic fields. This calculation can be thought of as one step of the fixed point iteration described in the Introduction.

Both the QA and QH configurations are scaled to a minor radius of 1.70 m and volume-averaged B of 5.86 T, matching the corresponding values of the ARIES-CS reactor study.¹⁹ Here and throughout this paper, density and temperature profiles are specified using

where $nα,0$ and $Tα,0$ are the on-axis values, and $s∈[0,1]$ is the toroidal flux normalized by its value at the plasma boundary. In these first calculations, we consider a hydrogen plasma with $ne,0=nH,0=4.13×1020$ m⁻³ and $Te,0=TH,0=12$ keV. (Numerical calculations of the bootstrap current in a hydrogen plasma are nearly indistinguishable from calculations with a deuterium-tritium mixture.) Here and throughout the remainder of the paper, we assume that there are no impurities, so the effective ion charge is Z_eff = 1, although the bootstrap current formula by Redl et al. can in fact account for impurities. The results of the modified Redl formula for the two configurations are shown in Fig. 1.

To assess the accuracy of the Redl formula, the bootstrap current is also calculated using the code SFINCS.²⁰ This code solves the drift-kinetic Eq. (4) in a stellarator, making no assumption of quasi-symmetry, and, therefore, fully accounting for three-dimensional geometry and imperfections in the quasi-symmetry. The bootstrap current computed by SFINCS using (6) is displayed in Fig. 1.

The drift-kinetic equation solved by SFINCS can also include the effect of a radial electric field through the term $vE·∇fα,1$⁠, where $vE=(dΦ/dψ)B−2B×∇ψ$ is the $E×B$ drift and $Φ(ψ)$ is the electrostatic potential. However, in the SFINCS calculations here, it is very difficult to accurately determine an ambipolar radial electric field, since for any value of the electric field the radial current is found to be within discretization error of zero. This “intrinsic ambipolarity” is a consequence of the precise quasi-symmetry of these configurations.²¹ Another consequence is that the bootstrap current computed by SFINCS is independent of the radial electric field to high precision. Therefore, for all SFINCS calculations in this paper, we set the radial electric field to zero.

It is apparent from Fig. 1 that the generalized Redl formula is reasonably accurate when compared to the first-principles SFINCS calculations. The small disagreement can arise, in principle, both due to imperfections in the quasi-symmetry and due to the fact that even in axisymmetry, the Redl formula is only an approximate fit to solutions of the drift-kinetic equation. Additional SFINCS calculations for the two configurations were carried out based on a Boozer-coordinate representation for the magnetic field, filtering out all Fourier modes of B that break quasi-symmetry. The results are nearly indistinguishable from the SFINCS calculations that include the imperfections in the quasi-symmetry. Therefore, for these geometries, the disagreement in Fig. 1 is dominated by the fact that the Redl formula only approximates true drift-kinetic solutions. This small disagreement is acceptable for the purposes of optimization, and it is within the large uncertainty associated with predicting the pressure profile. Moreover, as will be shown later, this small difference in the bootstrap current can be eliminated by applying a few steps of fixed-point iteration at the end of an optimization. Given the success of this initial test of the Redl formula in stellarators, we now proceed to use it in the optimization of new stellarator configurations.

To optimize for quasi-helically symmetric equilibria, we minimize the objective function

where A is the aspect ratio, a is the minor radius, $B¯$ is the volume-averaged field strength, and quantities with a subscript $*$ are target values. The term f_QS represents the departure from quasi-symmetry and is

The term f_QS or similar expressions have been used previously to optimize for quasi-symmetry.^2,22,23 Finally, the term f_boot penalizes inconsistency in the bootstrap current

Thus, f_boot = 0 if the bootstrap current is self-consistent, whereas if the parallel current in the MHD equilibrium is zero (as it may be in the first iteration), the normalizing denominator in (15) results in f_boot = 1. Overall, the objective (13) promotes quasi-symmetry and bootstrap current self-consistency, while keeping the aspect ratio, device size, and average field strength close to desired target values. For all results shown here, we match the minor radius and field strength of ARIES-CS, $a*=1.70$ m and $B¯*=5.86$ T. We did not find it necessary to introduce weights in front of each term in 13.

For quasi-axisymmetry, one more term is included in the objective

where $ι¯=∫01ds ι$ is the average rotational transform, and $ι¯*$ is a target value. The additional term in quasi-axisymmetry is required to prevent the optimum from being axisymmetric, i.e., a tokamak.

The integrals in the objective are discretized as finite sums, yielding a nonlinear least squares problem. In the quasi-symmetry term, the flux surface averages are discretized using uniform grids in the standard toroidal angle $ϕ$ and in a poloidal angle ϑ, which need not be the Boozer angle. We can then write

where

with $Δϑ$ and $Δϕ$ denoting the grid spacing, $gs=1/(∇s·∇ϑ×∇ϕ)$ is the Jacobian of the ${s,ϑ,ϕ}$ coordinates, and $Vs′=ΔϑΔϕ∑ϑk,ϕℓgs$⁠.

The parameter space for the optimization consists of the shape of the toroidal boundary surface, the enclosed toroidal flux, and a profile of either toroidal current or $ι(s)$⁠. The boundary shape is represented as

and the mode amplitudes $Rm,n$ and $Zm,n$ are included in the parameter space.

Fixed profiles of density and temperature are chosen, used both for the bootstrap current calculations and to set the pressure profile for the MHD equilibrium calculations. For MHD equilibrium, one other profile must be specified, typically either the toroidal current $It(s)$ (or its radial derivative) or ι. We find, in practice, that using the profile of $dIt/ds$ in the parameter space is preferable as it yields slightly lower values of the objective than if $ι(s)$ is used. However, varying $ι(s)$ in the optimization is also effective. The $dIt/ds$ profile is represented as a cubic spline on fixed nodes, with the function values varied during the optimization. Since typical values of $dIt/ds$ in SI units can be $∼106$⁠, the optimizer performs better if the function values are scaled so they are $∼1$⁠.

The optimization problem is solved using the SIMSOPT software,²⁴ with VMEC²⁵ as the equilibrium code. The nonlinear least squares problem is solved using the default algorithm in scipy, “trust region reflective.” Gradients are provided by finite differences, using MPI for concurrent function evaluations. The initial condition is a circular cross-section torus. For quasi-helical symmetry, torsion of the magnetic axis is added to the initial condition via $R1,1$ and $Z1,1$⁠. The parameter space is expanded in a series of steps, with the surface shape and current profile refined at each step. In step $j=1,2,…$⁠, the amplitudes $Rm,n$ and $Zm,n$ with $m≤j$ and $|n|≤j$ are varied, and $2j+3$ spline nodes are used in the current profile. For quasi-axisymmetry, the approach is slightly different and starts from the QA vacuum configuration presented in Ref. 2 to ensure that ι is sufficient to support the pressure and that the local minimum is not a tokamak. As a crude method to address the presence of local minima in the objective, for each case shown in Sec. VI, 12–24 optimizations are run with different finite difference step sizes, and the result with lowest objective function is shown. The software, scripts, input files, and output files associated with the optimizations are available in Ref. 26.

The optimization problem described here is limited to a single density profile and single temperature profile. These profiles would be chosen to be the target operating point of an experiment. However, in a real experiment it would be important to consider how the configuration changes for other profiles, particularly for the low temperature at the start of the discharge. Ideally, one would like to optimize for confinement and ι across a range of profiles. However, it is not clear how the procedure of this section could be generalized to account for multiple sets of possible profiles, since the variation of a configuration with profiles requires a choice of coils and free-boundary calculations, which are not considered here. For experimental scenarios in which the bootstrap current is different from its design value, such as the ramp-up phase at the start of a discharge, associated differences in ι could perhaps be compensated to some extent by other sources of current (driven by radio frequency waves or beams) or control coils. Understanding how to optimize a stellarator design across a range of plasma profiles is an important question for future research.

As a first demonstration of the new optimization method, we consider quasi-helical symmetry in a reactor-scale configuration with four field periods and an aspect ratio $A*=6.5$⁠. We adopt the density and temperature profiles (12) with $ne,0=2.2×1020$ m⁻³ and $Te,0=Ti,0=10$ keV.

The configuration resulting from the optimization has a volume-averaged β of 2.5% and a plasma current of 1.2 MA. The geometry of the configuration is shown in Fig. 2. The geometry is qualitatively similar to previous quasi-helically symmetric stellarators.^2,27–29 Contours of B in Boozer coordinates are shown for several flux surfaces in Fig. 3, showing excellent quasi-symmetry throughout the volume. Indeed, the neoclassical transport coefficient $εeff3/2$ is $<3×10−5$⁠, roughly 30 times smaller than in W7-X. For comparison, a recent estimate³⁰ suggests that $εeff3/2≲10−3$ is sufficiently small for a reactor, and for smaller values, turbulent transport will dominate.

Figure 4 shows the profile of bootstrap current in the optimized configuration, computed in three ways. The solid green curve shows the profile associated with the MHD equilibrium, taken from the VMEC output file. The dashed blue curve shows the result of the Redl formula (as modified in Sec. III). These two curves overlap closely, as expected due to the minimization of f_boot. Finally, the dotted red curve shows a SFINCS calculation for the equilibrium, fully accounting for three-dimensional geometry and imperfections in the quasi-symmetry, and considering the ions to be an equal mixture of deuterium and tritium. This curve also agrees well with the MHD equilibrium, showing that self-consistency of the current profile has been achieved.

A problem arises when this method is applied to larger values of β, which can be understood from Fig. 5. Here, the profile of ι is shown both for the $β=2.5%$ configuration and for a configuration optimized at the same aspect ratio at β = 0. The bootstrap current in quasi-helical stellarators acts to reduce $|ι|$⁠, and for $β=2.5%$⁠, the ι profile is close to crossing ι = 1. This value is the worst possible rational value, in that tiny field errors are likely to introduce large magnetic islands. For larger values of β, a method of controlling ι must be introduced to avoid crossing ι = 1.

To optimize for quasi-helical symmetry at higher β, the following term is added to the objective function:

This term acts as a barrier to prevent ι from dropping below 1.03, providing some margin away from the resonance at ι = 1.

An example of optimization using this barrier term at $β=5%$ is shown in Fig. 6. For this optimization, we again consider an aspect ratio $A*=6.5$ and the profiles (12), now with $ne,0=3×1020$ m⁻³ and $Te,0=Ti,0=15$ keV. Excellent quasi-symmetry throughout the volume is again achieved, shown by the contours of B in Boozer coordinates in Fig. 7.

Figure 8 shows that the rotational transform profile avoids the ι = 1 resonance. For comparison, a similar optimization with the same density and temperature profiles but without the barrier term (20) is also shown. In this case, the ι = 1 resonance is crossed in the middle of the volume. Interestingly, the barrier has the effect of shifting the $ι(s)$ profile by a constant, without affecting the magnetic shear.

Figure 9 shows the current profile for the final configuration, computed in the same three ways as in Fig. 4. Again, it can be seen that the current profile is self-consistent, with only a tiny difference between $⟨j·B⟩$ in the MHD equilibrium compared to the SFINCS calculation.

If there is any concern about the small remaining differences between the MHD and drift-kinetic current profiles, the difference can be made arbitrarily small in a post-processing step. Following the optimization, a small number of fixed-point iterations between MHD equilibrium and the drift-kinetic equation can be applied, as described in the Introduction. The boundary shape is fixed during these iterations. While there is no control on other terms in the objective function such as quasi-symmetry or $fι$ during the fixed-point iterations, the change to the current profile is extremely small, since the current profile was very nearly self-consistent before the fixed-point iterations. Therefore, degradation of the other objectives is negligible.

Fixed-point iterations using SFINCS are applied to the configuration of Fig. 6, and the final current profile is shown in Fig. 10. It can be seen that there is now essentially exact agreement between MHD equilibrium and the drift-kinetic equation. Figure 11 shows B as a function of Boozer angles after the fixed-point iterations. There is no significant degradation in quasi-symmetry compared to Fig. 7, before the fixed-point iterations. The final configuration has $εeff3/2<6×10−5$⁠, so neoclassical transport will be negligible.

For QA configurations, we start the optimization from the “new QA” geometry in Ref. 2. This is an $nfp=2$ vacuum configuration with an average ι of 0.42 and an aspect ratio A = 6.0. Like the other configurations in this paper, we have a = 1.70 m and $B¯=5.86$ T. We use profiles (12) with $ne,0=2.38×1020$ m⁻³ and $Te,0=Ti,0=9.45$ keV. The target average ι is achieved by including a term (16) in the optimization. To prevent ι from crossing 0.5, a barrier term as in (20) is employed, with the barrier at $ι=0.485$ and a maximum employed instead of a minimum to penalize crossing above this ι. No stepping in the number of Fourier-modes is done from this initial condition, as this results in the optimizer finding an axisymmetric tokamak as the local minimum.

The resulting configuration has a volume averaged β of 2.5% and a 2.72 MA plasma current. This value is $2.3×$ larger than the plasma current for the QH configuration at the same β in Sec. VI A, consistent with all neoclassical effects being larger in QA compared to QH symmetry. The neoclassical transport coefficient $εeff3/2$ is $<2×10−5$⁠, which again is more than adequate. The surface shape, bootstrap current profile, the quality of the quasi-axisymmetry, and the ι profile are shown in Figs. 12–15. The figures show that the boundary has modest shaping, the Redl formula is accurate, reasonable quasi-axisymmetry is achieved, and the ι profile avoids 1/2.

The ι profile does, however, cross several other low-order rationals, including 1/4, 2/5, and 1/3, but these crossings are confined to the core region, where the magnetic shear is relatively large, reducing the width of islands even if resonant fields are present. It may be difficult to entirely avoid crossing low-order rational surfaces in QA stellarators at finite β with self-consistent bootstrap current, so having these crossings in the core might be the best achievable outcome. The high shear across these rationals also limits the width of the region with low ι, which is beneficial for fast particle confinement as high ι reduces the orbit width.

The downside of higher ι is that it increases elongation. For vacuum QA configurations, this effect appears to limit ι to values of $ι≲0.9$⁠, based on previous studies exploring the landscape of quasi-symmetric stellarators.³¹ In practice, the space of vacuum QA configurations may be yet more limited, as the plasma boundary tends to become unreasonably elongated already for mean ι of about 0.6.

Adding finite β to these configurations, in the same manner as was done here with using the vacuum configuration as an initial guess in the optimization, does not make the plasma boundaries less elongated in our optimizations.

In general, trying to increase β much beyond 2.5% for quasi-axisymmetric configurations, trying various choices of $ι¯*$⁠, initial conditions, and optimization sequences, we have so far encountered three different possible modes of failure. (1) The optimal configuration can have an unreasonably elongated boundary, e.g., the ratio of boundary surface height to width can exceed 10 in the bean-shaped cross section. This is probably unacceptable since the alpha particle Larmor radius can become a sizeable fraction of the shorter dimension, and the very large pressure gradient in real space could drive instabilities. (2) The optimal configuration can become close to a tokamak, with most of ι being created by the bootstrap current. (3) The configuration may not be numerically converged: high-resolution VMEC calculations may show interior flux surface shapes that appear unphysical, being far from elliptical in the core. In our optimizations, failure state (1) is the most common. Failure state (2) is less common for larger values of $ι¯*$⁠.

It remains possible that for some other choices of parameters and objective function, configurations with self-consistent bootstrap current and a similar degree of quasi-axisymmetry may exist at higher β and/or at higher ι. In the future, other functions in the objective could be explored for avoiding axisymmetric solutions and solutions with large elongation.

Previously, a different approach for computing the bootstrap current in stellarators has been to use a formula derived by expanding in low collisionality, i.e., long mean-free-path compared to system size.^32,33 Variations of this formula have been derived by several authors.^9–11 Due to the low-collisionality expansion, these formulas allow the bootstrap current to be computed faster than if the finite-collisionality 3D drift-kinetic equation is solved (as in the codes SFINCS and DKES³⁴). However, for quasi-symmetric stellarators, the Redl formula is superior to these low-collisionality formulas for several reasons: the Redl formula is faster to evaluate, it is more accurate, and the low-collisionality formulas are plagued by noise.

These issues are shown in Fig. 16. Here, three methods for computing the bootstrap current are compared for the quasi-axisymmetric configuration of Sec. VI C with $β=2.5%$⁠. The red dotted curve shows a calculation of the full 3D drift-kinetic equation with SFINCS, which is the most complete physics model. The Redl formula is shown in the blue dashed curve, displaying close agreement with SFINCS. The orange curve shows the low-collisionality bootstrap formula computed with the code BOOTSJ.^10,35,36 This curve is extremely noisy, looking nothing like the Redl or SFINCS results. BOOTSJ includes an ad hoc smoothing parameter named “damp_bs,” denoted here by d. This parameter is defined as follows: resonances $1/x$ in the bootstrap formula are replaced with $x/(x2+m2d2)$⁠, where $x=m−n/ι$⁠, and m and n are the poloidal and toroidal mode numbers in the solution. The green and pink curves in Fig. 16 show BOOTSJ results for d = 0.01 and d = 1. It can be seen that even though the spikes of the orange curve can be smoothed over with a sufficiently large d, the current profile still does not match the SFINCS or Redl results. The reason may be in part that the bootstrap current approaches the low-collisionality asymptote extremely slowly, as shown for example in Fig. 3(b) of Ref. 37. In contrast, the Redl formula is valid for arbitrary collisionality. Another contributing factor to the difference may be finite aspect ratio effects not included in BOOTSJ.

In summary, for computing the bootstrap current in quasi-symmetric stellarators, the Redl formula has many advantages over low-collisionality asymptotic formula.

Confinement of fusion-produced alpha particles has long been a serious challenge for the stellarator concept. Several QA and QH optimizations have been recently demonstrated to achieve extremely good confinement of fast ions^2,38–40 in vacuum fields, but a natural question is whether similarly good confinement is possible at finite beta. We can now answer this question in the affirmative: the stellarator configurations obtained in the preceding sections have excellent alpha confinement, as a result of the high degree of quasi-symmetry. Both the problem and the recent progress toward its resolution are shown in Fig. 17. This chart displays the fraction of alpha particle energy that is lost for a variety of stellarator configurations, using calculations similar to the recent study in Ref. 41. Each configuration has been scaled to the minor radius and averaged field strength of ARIES-CS. In each configuration, 20 000 alpha particles are initialized isotropically in velocity space at 3.5 MeV, with a spatial density proportional to the local fusion reaction rate. The latter is determined using the same profiles for the density and temperature of deuterium and tritium in each configuration, $nD=nT=(2×1020 m−3)(1−s5)$ and $TD=TT=(12 keV)(1−s)$⁠. The alpha guiding centers are followed with collisions until a particle either thermalizes with the bulk plasma or moves outside the s = 1 surface, whereupon it is considered lost (a conservative assumption). The calculations are performed with the ANTS code.⁴² ANTS uses the relativistic guiding center equations on pages 42–43 of Ref. 43. Test particle collisions with the bulk plasma are computed using the same profiles as for the alpha birth distribution, along with an electron species with $ne=nD+nT$ and T_e = T_D. These profiles are not generally consistent with the pressure profile used to compute each MHD equilibrium. However, to fairly compare alpha confinement across configurations with varied β, it is more important that the birth and collision profiles be matched.

The set of magnetic configurations includes W7-X,⁴⁴ NCSX,⁴⁵ ARIES-CS,¹⁹ CFQS,⁴⁶ HSX,²⁷ and LHD.⁴⁷ For LHD, since it is known that confinement is sensitive to shifts in the major radius, configurations with two values of major radius are included. The configurations from Refs. 2, 29, 38, 40, 41, and 48–50 are also included.

Figure 17 shows that for the majority of the stellarators designed before the year 2021, alpha energy losses of tens of percent are typical. These losses would likely cause significant damage to the plasma facing materials. For comparison, a perfectly axisymmetric ITER hybrid configuration is included at the bottom of the figure, with losses of only $∼1%$⁠. (Although this ITER configuration is perfectly symmetric, losses are still possible due to direct orbit loss and collisional diffusion.) It is interesting that the LHD configuration with smaller major radius has lower losses than many stellarators that are nominally optimized for reduced neoclassical transport. However, the losses are significantly lower, $∼2%$ or smaller, for the configurations from Sec. VI of this paper, along with the other configurations reported in the past year. It is interesting that the QH configurations can have losses lower than the ITER configuration. This is due to the fact that the width of banana orbits in an axisymmetric or quasi-symmetric configuration scales as $∝1/|ι−N|$⁠. Therefore, the banana orbits are significantly thinner in the QH configurations, leading to reduced alpha loss. Note that the finite width of banana orbits is independent of the bounce-averaged drifts, so confinement measures based on bounce averages such as Γ_c in Ref. 51 will not account for this difference between QA and QH configurations.

Comparing alpha particle confinement between stellarators and tokamaks is complicated, so the ITER bar in Fig. 17 must be interpreted with caution. On the one hand, higher field strengths may be possible in tokamaks, since greater structural support would be needed in a stellarator to prevent bending deformations of non-planar coils under electromagnetic loads. This difference favors better alpha confinement in tokamaks. On the other hand, stellarators can operate at higher collisionality, due to the absence of the Greenwald density limit, and the shorter slowing-down time would favor confinement in stellarators.

In this work, we have demonstrated a new method for computing the bootstrap current in quasi-symmetric stellarators, leveraging the isomorphism with tokamaks. The method enables the optimization of quasi-symmetric stellarators such that the current profile in the MHD equilibrium is consistent with the drift-kinetic equation. Since the bootstrap current is computed in this method using explicit analytic formulas, the computational cost of the bootstrap calculations is negligible compared to the MHD equilibrium calculations. The new approach has the additional advantage that the current profile is smooth, in contrast to the noisy low-collisionality asymptotic formulas for the bootstrap current in general, stellarators.

For equal minor radius, field strength, and profiles of density and temperature, the new configurations have alpha particle confinement comparable to or even better than a tokamak. Similarly, excellent confinement of alpha particles is seen for other β = 0 stellarators reported in recent months.^2,38,40

Although quasi-symmetry imperfections in the configurations here are consistently smaller than in pre-2021 configurations, the imperfections at finite β are larger than in vacuum configurations at the same aspect ratio. In future work, it would be valuable to understand this phenomenon. One explanation could be the following. As discussed in Ref. 52, in an expansion in small distance from the magnetic axis compared to the scale length of the axis shape, quasi-symmetry is likely to be possible to two orders but not three. Therefore, quasi-symmetry is likely best when gradients in the minor radius direction are small, so the third-order terms are as small as possible. The bootstrap current profile introduces relatively short-scale variation in the minor radius direction, potentially increasing the size of the symmetry-breaking third-order terms. Supporting this theory, preliminary numerical experiments show that for fixed total toroidal current, quasi-symmetry errors are reduced if an s-independent current profile is prescribed, in contrast to the self-consistent bootstrap current profiles used here. In the future, further investigations of the β-dependence of quasi-symmetry errors would be worthwhile.

Significant additional work would be needed before the configurations in this paper could be considered complete conceptual designs. First, MHD stability would need to be assessed and possibly optimized. Second, flux surface quality would need to be checked using an MHD or extended-MHD code that does not assume the existence of flux surfaces. While both of the QH configurations avoid the low-order rationals at ι = 1 and 4/3, they both cross $ι=8/7≈1.14$ and $12/11≈1.09$⁠, and the $β=5%$ case crosses $ι=12/10=1.2$⁠. Islands could exist at any of these resonances. Third, coil feasibility would need to be studied.

The most uncertain aspects of these configurations are the density and temperature profiles. In this work, we have considered only prescribed profiles. In a real experiment, these profiles would largely be determined by the sources of heat and particles together with turbulent transport, which is challenging to model. Due to this significant uncertainty, the sensitivity of the confinement to different density and temperature profiles should be studied.

Another question for future work is the extent to which, within the constraint of QA or QH symmetry, the bootstrap current magnitude can be minimized. Minimization of the current would help to preserve the other optimized properties over a wider range of density and temperature profiles, and reduce a potential source of MHD instability. While the quasi-symmetry isomorphism implies that the bootstrap current in QA or QH symmetry cannot be eliminated, some freedom does exist in its magnitude. For given density and temperature profiles, the bootstrap current can be varied through ι and f_t, the latter of which can, in principle, be controlled through the ratio $Bmax/Bmin$ on the surfaces.

Magnetic configurations in Fig. 17 were provided by Andrew Giuliani, Sophia Henneberg, Geoffrey McFadden, Shoichi Okamura, Neil Pomphrey, John Schmitt, Don Spong, Yasuhiro Suzuki, and Florian Wechsung. Thanks also to Andreas Redl and Aaron Bader for providing clarifying input related to this research. Support from the SIMSOPT team is gratefully acknowledged. This work was supported by the U.S. Department of Energy under Contract Nos. DE-AC02-09CH11466 and DE-FG02-93ER54197. This work was also supported by a grant from the Simons Foundation (No. 560651, M.L.).

The authors have no conflicts to disclose.

Matt Landreman: Conceptualization (lead); Investigation (lead); Writing – original draft (lead). Stefan Buller: Investigation (supporting); Writing – original draft (supporting). Michael Drevlak: Software (supporting).

The data that support the findings of this study are openly available in Zenodo at https://doi.org/10.5281/zenodo.6520103, Ref. 26.

In this section, to avoid confusion with the conventional tokamak notation for the poloidal flux, we denote the toroidal flux as ψ_t. Assuming zero induced electric field (⁠$E||=0$⁠) and applying the substitutions in Sec. III, the Redl formula for a quasi-symmetric stellarator is, from Eq. (2) in Ref. 6,

The coefficients are from Eqs. (10)–(16) and (19)–(21) in Ref. 6 

The trapped particle fraction, denoted $ftrap$ in Ref. 6, is

The effective collisionalities are

The Coulomb logarithms are⁸ 

with the temperatures in eV and the densities in m⁻³.

In (A1), if light impurities are included in the calculation, the terms should be substituted as

although this is only an approximation, as it assumes that the impurity and bulk ion transport coefficients (described by $L31$ and α) are the same. This assumption is always somewhat violated in practice, if the impurities have different collisionalities than the bulk ions, but this can be expected to be partly compensated for by the fact that all the transport coefficients have been fitted to simulations with carbon impurities at different $Zeff$⁠.