Optimization of quasi-axisymmetric stellarators with varied elongation

Quasi symmetry has been proven to be a powerful concept for improving neoclassical transport in stellarators. The idea came from the original work of Boozer¹ who showed that the particle orbits in 3D stellarator geometry only depend on the magnitude of the magnetic field in Boozer coordinates. Thus, the particle orbits in stellarators are equivalent to those of tokamaks if the magnetic field strength is symmetric in Boozer coordinates even though the vector magnetic field and geometry of magnetic surfaces are 3D. Using this idea, Nührenberg and Zille² demonstrated numerically the existence of quasi-helically symmetric (QH) stellarators. QH means that in the Fourier expansion of the magnitude of magnetic field on a magnetic flux, $B=∑m,nBm,n cos (mθ−nϕ)$⁠, the dominant Fourier components $Bm,n$ are of only one m/n ratio, where θ and $ϕ$ are poloidal and toroidal Boozer coordinates, respectively. This led to the Helical Symmetric eXperiment (HSX) which validated the theory of QH experimentally.³ Furthermore, the existence of Quasi Axisymmetric (QA) stellarators was also shown,^4,5 where QA means that the dominant Fourier components $Bm,n$ are only of n = 0. Since then, a number of QA stellarators have been designed including NCSX,⁶ CHS-qa,⁷ ESTELL,⁸ and more recently, CFQS⁹ and the new QA stellarator by Henneberg et al.¹⁰ In particular, CFQS is currently being built in China. A reactor concept based on an optimized QA configuration has also been designed.¹¹ Compared to tokamak reactors, QA stellarators have advantages of steady state operation without an external current drive and potentially disruption free operation. The freedom of 3D geometry also allows more external control over important properties such as MHD stability and plasma transport while maintaining quasisymmetry.¹² 

In this work, we have carried out an optimization study of QA stellarators with varied elongation. It is known that, for tokamaks, the elongation has strong effects on the MHD stability¹³ and plasma transport.¹⁴ Thus, it is important to investigate the dependence of QA stellarators on elongation. The present work starts from a three field period aspect ratio = 6 QA configuration obtained by Ku and Boozer¹⁵ who performed a systematic study of the dependence of modular coil geometry on the aspect ratio, number of field periods, external rotational transform, and plasma-coil spacing. Our work here, on the other hand, is focused on the effects of elongation. This work extends the previous work to a range of elongation at both zero beta and finite beta with self-consistent bootstrap current. The elongation scan is carried out both at a fixed plasma volume and with volume proportional to elongation. The latter scan is motivated by the fact that for axisymmetric tokamaks, the plasma volume is proportional to elongation at a fixed major radius and minor radius. Our results demonstrate that good QA stellarators exist with a range of elongations that are stable to external kink modes.

This paper is organized as follows. Sections II and III present the results of an elongation scan at a fixed volume and at a volume proportional to elongation. Section IV describes the calculation of the self-consistent bootstrap current. Section V presents the results of MHD stability optimization. Finally, conclusions are given in Sec. VI.

Here, we carry out an elongation scan starting from the three field period aspect ratio = 6 QA configuration obtained in the work by Ku and Boozer.¹⁵ Elongation was not a target function in the original work, but rather an output quantity with plasma volume as a constraint. We use the STELLOPT code¹⁶ to do the optimization using the genetic algorithm with differential evolution (GADE)^17,18 at each value of elongation in the zero beta case. In the optimization, the sum of $χi2=(fitarget−fi)2/σi2$ is minimized, where f_i is the ith variable chosen to be targeted, $fitarget$ is its targeted value, and σ_i is a weight parameter. We consider several critical physical properties listed in Table I.

Here, the elongation κ is defined by the shape of the plasma boundary at the $ϕ=0$ plane. Specifically κ is the total height divided by the width at the mid-plane of the last closed flux surface; $ϵeff3/2$ is the effective helical ripple calculated by the NEO code;¹⁹ s is the normalized toroidal flux; QA is enforced (helicity constraint) by assigning a target value of zero to each calculated non-axisymmetric component of field strength in Boozer coordinates, where $Bm,n$ harmonics with $n=0−16$ and $|m|=0−32$ are included, and identical weights are assumed for each targeted harmonic. The weights are chosen by trial and error to give the desirable results.

The following figures show the evolution of one optimization calculation targeting elongation to be $κ=3.64$ (⁠$κ=3.84$ is the initial configuration, from Ref. 15).

Figure 1(a) shows the evolution of the total $χ2$ and its components. We observe that the total $χ2$ decreases to a constant value after approximately 100 iterations with $χi2$ of the neoclassical transport being the dominant component. Figures 1(b)–1(d) show the evolution of the effective ripple (⁠$ϵeff3/2$⁠) profile, the iota profile, and the flux surface shape. We observe that the effective ripple decreases by a factor of two after optimization. The ι profile is almost identical to the initial one. The plasma boundary shape at $ϕ=0$ shows a reduction of elongation to its target value of 3.64. Figure 1(e) plots the dominant Boozer $|Bnm|$ components as a function of normalized toroidal flux. It indicates that the optimized configuration is indeed a QA stellarator since the axi-symmetric components are much larger than the non-symmetric components. Figures 1(f) and 1(g) show the $|B|$ patterns on the last closed flux surface for the initial and final cases. The results indicate that the magnetic field strength of the final case is nearly axi-symmetric and its symmetry is a little improved over that of the initial case.

The above optimization process is repeated for each target elongation value. This is done by decreasing the target value of elongation in each run by small increments. The results of this elongation scan are shown in Fig. 2 which plots the flux surfaces for several representative values of elongation. Figure 3 shows the corresponding effective ripple (⁠$ϵeff3/2$⁠) profile for the different elongation cases. The results suggest that the neoclassical transport can be kept at low levels when the elongation is reduced to lower values. It should be noted that the effective ripple level of the original configuration is significantly higher because it was optimized at finite plasma beta whereas here the zero beta case is considered. Table II shows the data of major physical quantities for different elongation cases after optimization, where R_major is the averaged major radius calculated from the volume and cross section, a_minor is the averaged minor radius calculated from the cross section area, $⟨B⟩$ is the volume averaged magnetic field strength, and B₀ is the averaged magnetic field strength on axis. These quantities are kept almost the same for different elongation cases. This is achieved by targeting the aspect ratio, volume, and the rotational transform to be the same as in the original case. The edge toroidal flux is held fixed during the optimization, leading to a slight variation in the magnetic field strength.

Here, we carry out a different elongation scan at zero beta by targeting the plasma volume proportional to κ while other targets are kept the same as before. This is motivated by the fact that for an axisymmetric tokamak geometry, the volume is proportional to κ.

Figures 4 and 5 show the optimized outer boundary flux surfaces and the corresponding effective ripple (⁠$ϵeff3/2$⁠) profiles for several values of elongation. We see that the neoclassical transport levels are still kept at low levels.

There is a question on how to define κ values for 3D stellarators. In this work, we use κ at $ϕ=0$ plane as a measure of elongation. On the other hand, we can also use the averaged elongation, i.e., the elongation of the corresponding tokamak geometry defined by the n = 0 components of the 3D stellarator geometry. Figure 6 plots the volume vs $R0a02⟨κ⟩$ for the stellarator and its corresponding tokamak.

R₀, a₀, $⟨κ⟩$ is the major radius, minor radius, and elongation of the corresponding tokamak, respectively. We observe that the volume of the corresponding tokamak is proportional to the product of $R0a02⟨κ⟩$⁠. The volume of the stellarator is a little smaller than that of the corresponding tokamak, but it is indeed proportional to $R0a02⟨κ⟩$ approximately.

In quasisymmetric stellarators, the bootstrap current²⁰ is important to consider, as in axisymmetric tokamaks, since the physics of bootstrap current is the same in the Boozer coordinates. It can significantly modify the equilibrium properties through change in the rotational transform. Here, we use the drift-kinetic code SFINCS²¹ to calculate the bootstrap current. The SFINCS code solves the steady-state drift-kinetic equation for multiple species, allowing arbitrary collisionality and magnetic geometry (subject to the assumption that magnetic surfaces exist), and using the full linearized Fokker–Planck–Landau collision operator. For calculations shown here, the variation of the electrostatic potential on magnetic surfaces is neglected. The SFINCS code is state of art for calculation of the bootstrap current, however it is relatively expensive computationally and cannot be used directly in the optimization process. Thus, here we use it iteratively in our optimization. Specifically, we first calculate the bootstrap current at the start of the optimization and then keep it fixed during the optimization process using STELLOPT. After each optimization, we calculate bootstrap again using the optimized configuration and then repeat the optimization process using the new bootstrap current profile. This iterative process is repeated until the final configuration converges. We will show that this iterative method works for the cases considered in this work.

In our bootstrap current calculation, the pressure profile is given by $P0(1−s2)3$⁠, where P₀ is a constant adjusted by STELLOPT to achieve the target value of volume averaged beta and s is the normalized toroidal flux. The ion and electron density profiles are uniform at $4×1020 m−3$⁠, and the temperature profiles of the two species have the same shape as the pressure profile. Figure 7 plots the current profile and total bootstrap current calculated by SFINCS for different β values for the lowest elongation $κ=2.57$ case in Fig. 2. We observe that the dependence of the bootstrap current profile on the β value and pressure profile is consistent with theoretical expectation.²² It should be noted that the central iota decreases as beta increases. This is due to the effect of beta on the equilibrium.

Figure 8 shows the total bootstrap current as a function of κ, holding the plasma $β=0.05$⁠. The curve for the constant volume series indicates that the total bootstrap current increases a little when the elongation decreases. The other curve shows that the bootstrap current is approximately proportional to elongation or volume. In tokamaks, the bootstrap current density is given by $⟨J·B⟩=−cRBϕdP/dψ$ or $⟨J·B⟩=−c(Bϕ/Bp)dP/dr$⁠, where c is a constant coefficient, ψ is the poloidal flux, and B_p is the poloidal magnetic field.²² Thus, the elongation affects the total bootstrap current through the cross section area and the poloidal field. The poloidal B field in the stellarator is mainly supplied by external coils when the bootstrap current is small. In this limit, the bootstrap current is proportional to elongation since the cross section area or plasma volume is proportional to elongation. This explains why in Fig. 8, the total bootstrap current is roughly independent of kappa for the constant volume case, while for the varying volume case, the total bootstrap current is almost proportional to kappa (volume).

MHD stability is an important property for all magnetic confinement devices. Strong global MHD instabilities can seriously degrade the confinement of plasmas. Therefore, it is important to consider MHD stability in designing optimized stellarators.

We first consider the stability of MHD modes for the $κ=2.57$ case with a plasma beta value of $β=0.05$⁠. The equilibrium profiles of pressure, current, and rotational transform are shown in Fig. 9.

We see that the rotational transform profile is very different from the zero beta case [see Fig. 1(c)]. Iota increases radially from about 0.1 at the magnetic axis to 0.75 at s = 0.75 and then decreases to 0.6 at the edge. Figure 10 shows the five dominant Fourier harmonics for the n = 0 and n = 1 families calculated by the TERPSICHORE code.^23,24 We note that, for the n = 0 mode family, the largest component [(m, n) = (4, 3)] peaks near the $ι=0.75$ rational surface as expected. On the other hand, for the n = 1 mode family, the mode is clearly external with the dominant (2, 1) component peaks at the edge. This is consistent with the fact that iota = 0.5 is the nearest major rational surface in the vacuum. The normalized eigenvalues for the two families (n = 0 and n = 1) are $ω2=−2.269×10−2$ and $−4.988×10−2$⁠, respectively. The magnitudes of the eigenvalue indicate that the modes are strongly unstable. Thus, these modes need to be stabilized via optimization.

As done above, we keep the bootstrap current profile fixed as given in Fig. 9(b) although the plasma boundary will be changed as the optimization progresses due to the computational expense of calculating the bootstrap current by SFINCS. We use a Levenberg–Marquardt (L–M) algorithm^25,26 with finite differences for derivative evaluation in our optimization when targeting MHD stability. The total number of function evaluations required to achieve a satisfactory optimization is significantly less for L–M than for GADE, and trapping in local minima was not found to be a significant problem.

Figure 11(a) shows the evolution of key targets including kink stability in the STELLOPT calculation for the above case. We observe that the unstable kink growth rate squared value is reduced over four orders of magnitude, but the neoclassical transport level is increased somewhat. The later iterations are mainly optimizing the neoclassical transport. Figure 11(d) shows clearly that the negative curvature at the inner side of the flux surface at the $ϕ=0$ plane is increased, i.e., higher triangularity. This is consistent with our expectation that the negative curvature improves the MHD stability.²⁷ Figure 12 shows that it requires a larger negative curvature for the higher plasma β cases.

Since the flux surfaces are changed after optimization, we need to re-calculate the bootstrap current again to make it consistent with the new configuration, and check the MHD stability again. If it is kink unstable, we need to do the optimization again. The iterative process needs to be repeated a few times until the results are converged. Figure 13 shows the evolution of key targets during the iteration. Using the final case of Fig. 11, we calculate its self-consistent bootstrap current $jbs1$ (⁠$β=0.05$⁠) and find that it is kink unstable. After the 2nd optimization, the final case with its self-consistent bootstrap current $jbs2$ is much more stable (⁠$ω2=1.29×10−5$ and $1.64×10−5$ for n = 0 and n = 1 family, respectively), and the neoclassical transport level is also relatively low [Fig. 13(b)].

In this work, a systematic optimization study has been carried out starting from a three field period quasi-symmetric stellarator with an aspect ratio of 6. A series of zero beta QA configurations with varied elongation are generated using the stellarator optimization code STELLOPT. For finite beta stellarator plasmas, the bootstrap current is calculated self-consistently using the SFINCS code. The bootstrap current increases approximately linearly with elongation for QA stellarators, if the volume is proportional to elongation. Good QA stellarators are obtained with varied plasma beta, self-consistent bootstrap current, low neoclassical transport, and good kink stability. This study demonstrates that good QA stellarators with varied elongation and plasma beta exist. Future work will explore larger parameter space of quality QA stellarator configurations.

This work was supported by Zhejiang University's startup funding for one of the authors (Guoyong Fu) and by the U.S. Department of Energy under Contract No. DE-AC02-09CH11466.