Investigating nonlinear magnetohydrodynamics in an optimized, reactor-scale quasi-axisymmetric stellarator

The success of Wendelstein 7-X,¹ together with advances in optimization,² has driven a surge of interest in stellarators as a fusion power plant concept. Like tokamaks, stellarators rely on externally generated magnetic fields to confine plasma into topologically toroidal configurations.³ A defining characteristic of stellarators is the intentional breaking of the continuous toroidal symmetry, nominally present in tokamaks, i.e., $∂φ≠0$⁠, where $φ$ is the toroidal angle. This means that, unlike tokamaks, stellarators are not reliant on currents driven in the plasma to generate a confining magnetic field. Plasma current is a key driver of magnetohydrodynamic (MHD) instabilities, which can be detrimental to confinement. From this perspective, stellarators offer advantages with respect to macroscopic stability.

Magnetic field lines lying on a magnetic surface are confined to it. To the lowest order, charged particles in a plasma are confined to a magnetic field line. This can make continuously nested magnetic surfaces a desirable characteristic for global confinement in stellarator plasmas. Since magnetic field lines are described by a Hamiltonian system, when axisymmetry is broken, magnetic field structure is no longer guaranteed to consist of continuously nested magnetic surfaces. In other words, three-dimensional (3D) magnetic fields can be non-integrable.^3,4 Consequently, stellarator magnetic fields are notoriously difficult to design^5–7 precisely because the magnetic field structure is not known a priori. Despite this, there exist approximate symmetries of the magnetic field, known as quasi-symmetry (QS), that can reduce neoclassical particle transport and, therefore, improve confinement in 3D magnetic fields.⁸ This motivates, at least in part, the recent explosion of interest in quasi-symmetric configurations.^9,10

It is widely held that stellarators are more robust to macroscopic instabilities than tokamaks.² This is supported by experimental observations on the Large Helical Device (LHD), which have shown that linear MHD stability thresholds can be exceeded without significant degradation of energy confinement.¹¹ Earlier stellarator optimization efforts sought absolute stability against linear ideal MHD modes as part of the optimization criteria.^12,13 More recent studies, however, suggest that such a stringent stability constraint trades off against coil complexity.¹⁴ For the development of stellarators as a viable fusion power plant concept, this is problematic since coil complexity is a major driver of cost.^15,16 In contemporary optimization, a concerted effort has been made by the stellarator community to focus on QS, energetic particle confinement, and turbulent transport.² In this work, we analyze the nonlinear MHD properties of an optimized, finite-beta, reactor-scale stellarator equilibrium, designed in accordance with these principles.¹⁷ Our study demonstrates that, when optimizing for finite-β quasi-axisymmetric equilibria, MHD stability is at least as important a target for optimization. This motivates the development of more sophisticated MHD stability metrics than what is currently available.

The remaining sections are organized as follows: in Sec. II, we describe the methods used for the numerical simulations. In Sec. III, we present the results and analysis. Discussion and conclusions are in Sec. IV.

Our simulations use the high-fidelity, extended-MHD code, M3D-C1,¹⁸ that now accommodates strongly shaped, non-axisymmetric computational domains.¹⁹ To date, in stellarator geometry, the code has been compared to the PIES equilibrium code¹⁹ as well as to experiments on W7-X for simulating sawteeth-like phenomena.²⁰ We use a single-fluid model that includes heat transport, resistivity, and viscosity, with conducting wall boundary conditions. The full set of equations can be found in Jardin et al.¹⁸ For convenience, the single-fluid equations are repeated in  Appendix A. For the simulations presented, there is no vacuum region between the initial plasma boundary and vacuum vessel, i.e., we consider a limited plasma. Both the temperature and current are assumed to vanish at the wall. The density at the wall is fixed at its initial value and, therefore, acts as an implicit source. The shape of the computation domain is the same as the shape of the equilibrium plasma. The equilibrium under consideration was developed using the SIMSOPT optimization framework²¹ and designed for energetic particle confinement, quasi-symmetry, and self-consistent bootstrap current, assuming a magnetic field with continuously nested magnetic surfaces.¹⁷ According to theory, these three properties are critical for a range of high-performance stellarators; energetic particle confinement is critical for burning fusion plasmas,²² quasi-symmetry is a common strategy for reducing neoclassical transport of particles,⁸ while, in stellarators, bootstrap current can be stabilizing against magnetic islands.^23,24 On the other hand, MHD stability was not considered in this optimization. A widely used metric for linear ideal MHD stability in stellarators is the so-called magnetic well, which is the average curvature on a magnetic surface.²⁵ If a surface is linearly stable to perturbations of the local pressure gradient, then the surface is said to possess a magnetic well. For the equilibrium considered here, a posteriori analysis shows that there is a magnetic well for all surfaces.

The entire plasma is contained within the computational domain, and no symmetries (discrete or continuous) are imposed on the solutions. The simulations presented are run on a semi-structured 3D mesh with 72 toroidal planes, corresponding to $8.64×105$ 3D elements. The compute time for each simulation is approximately $4×105$ CPU-hours, corresponding to a wall-clock time of 14 days. Convergence studies were performed to identify optimal runtime parameters. No hyper-diffusion parameters are applied.

We start the M3D-C1 simulations from an initial state that has nested magnetic surfaces and follow the subsequent nonlinear evolution. In Fig. 1, we present Poincaré sections of the magnetic field [1(a)–1(c)] at $t=0τA$ and $t=825τA$ for both the constant and $η0=1$ Spitzer resistivity models. The full sequence of Poincaré sections for the constant resistivity model is shown in  Appendix B (see Fig. 6). The corresponding electron temperature profiles are shown in Figs. 1(d)–1(f). Overall, the equilibrium is highly unstable to multiple MHD modes that grow on Alfvénic timescales. We also find the plasma evolution to be qualitatively identical, regardless of resistivity model. This is indicative of an ideally unstable plasma. Modes with high poloidal and toroidal Fourier mode numbers, m and n, respectively, appear almost immediately near the plasma edge. Within approximately $500τA$⁠, a chaotic region with a thickness that is around 5%–10% of the minor radius is formed near the plasma boundary. Over the next $300τA$⁠, the chaotic layer grows considerably as additional moderate to high-(m, n) modes begin to develop on the remaining outer flux surfaces. By around $800τA$⁠, the chaotic layer has expanded to nearly 25% of the minor radius. The chaotic layer continues to expand until any confinement is effectively lost as the core is enveloped by a sea of chaotic magnetic fields. In Fig. 2, we show a Poincaré section for the constant resistivity case, as a representative example. We also observe the development of two core instabilities, with m = 5 and m = 3, that have much lower growth rates. Consequently, they do not reach a significant amplitude in the time simulated here. The simulations are stopped once a large volume of chaotic field has formed, as the plasma is no longer confined (discussed below). This occurs at different rates for each of the profiles. Hence, the simulations are stopped at different times. Since the unstable modes are internal in character, they are unaffected by the boundary conditions chosen for this study.

To explain the appearance of temperature gradients despite the presence of a significant volume of chaotic magnetic fields [see Figs. 1(b), 1(c), 1(e), and 1(f)], we can consider the pressure relaxation time (⁠ $τrelaxation$⁠). For the equilibrium under consideration, which has a minor radius $a=1.704 m$⁠, we note that that $τrelaxation∼a2/χ⊥=1.33 s$⁠. By contrast, $1500τA$ (c.f. Fig. 2) corresponds to $∼700 μs$⁠. Thus, temperature gradients remain because the magnetic surfaces break up so quickly that the temperature profile has not had time to relax. In the single-fluid MHD model, pressure is the product of the temperature and density. Hence, our observation is consistent with previous work,²⁸ which demonstrated, both numerically and experimentally, that pressure gradients can persist in regions where the magnetic field is chaotic. In our case, the observed breakup of magnetic surfaces is due to changes in the magnetic field that are produced by rapidly growing MHD instabilities. This is mechanistically distinct from the break up of magnetic surfaces due to currents induced by equilibrium pressure gradients.²⁸ 

To quantify the observed plasma dynamics, we Fourier decompose the kinetic energy and compute approximate growth rates for each toroidal Fourier mode by taking the derivative of the logarithm of the kinetic energy with respect to time. (M3D-C1 itself is not a spectral code.) Because the absolute values may be small, we apply a Gaussian kernel with standard deviation $±5 τA$ to get a relatively smooth estimate of the growth rates as a function of time. In Fig. 4, for each resistivity model, we plot growth rates as a function of time for even-numbered toroidal Fourier modes with $6≤n≤30$⁠. We find that at $∼200τA$⁠, a large number of these modes begin to develop. The largest growth rates correspond to the toroidal Fourier modes with largest n. To confirm this, in Fig. 3, we plot the maximum growth rate of even-numbered toroidal Fourier modes with $6≤n≤30$ for each resistivity model considered. For all resistivity models, the even-numbered toroidal Fourier modes with n > 6 have positive maximum growth rates, indicating that these modes are unstable. In each case, modes with the largest overall growth rates are the toroidal Fourier modes where $20≤n≤28$⁠.

The substantial number of toroidal Fourier modes that were destabilized can be explained by considering the N-fold rotational symmetry of the equilibrium. In this case, N = 2. This discrete symmetry leads to coupling of Fourier modes in the toroidal direction. Thus, an instability with toroidal Fourier mode number n₀ can couple to other toroidal Fourier modes with $nl=ln0±N$ for $l∈ℤ$⁠. For the N = 2 equilibrium, this implies that a single unstable eigenmode will be composed of numerous even-numbered toroidal Fourier modes, which is consistent with what is observed.

Finally, we remark on some differences in the nonlinear behavior between the constant and Spitzer resistivity profiles, which are temperature dependent. While evident in Fig. 4 for toroidal Fourier modes with $20≤n≤30$⁠, in Fig. 5, we plot the growth rate for just n = 30, for clarity. In Fig. 5, at first, the growth rate increases linearly for both the temperature dependent and independent resistivity profiles. At $t≈500τA$⁠, the growth rate for the constant resistivity profile plateaus, which is characteristic of the linear growth phase of MHD instabilities. At $t≈700τA$⁠, the growth rate decreases toward zero, which suggests the mode is tending toward nonlinear saturation. By contrast, for the temperature-dependent models, the growth rate does not plateau after the initial phase. Instead, we see a decrease rate of change of the growth rate followed shortly thereafter by a significant increase, before it peaks and tends toward zero. The secondary increase in the growth rate of the cases with temperature dependent resistivity is suggestive of strong nonlinear interactions, which are accelerating the growth of those modes. Further work is warranted to investigate the details of these interactions, which could be understood, for example, by undertaking an energy transfer analysis. Evidently, the Spitzer resistivity profiles are more realistic. So, although beyond the scope of the present work, these observations do highlight the potential need to proceed with caution when using linear approximations to predict nonlinear MHD behavior in stellarators.

The extended-MHD code, M3D-C1, has been applied to calculate the nonlinear evolution of a finite-beta, reactor-scale stellarator equilibrium that was optimized for quasi-axisymmetry, energetic particle confinement, and with self-consistent bootstrap current.

We performed a series of high-fidelity simulations to model the nonlinear MHD evolution of this initial state. We use a single-fluid MHD model that evolves the magnetic field, temperature, density, and flow. Importantly for stellarators, we do not impose constraints on the structure of the magnetic field. Thus, it (along with the other fluid variables) evolves self-consistently with the plasma dynamics that include effects due to resistivity, viscosity, and anisotropic heat conductivity. We use a reduced model for heat transport that captures the most crucial physics property, namely, strong anisotropy in the directions perpendicular and parallel to the magnetic field which has a ratio of $10−6$⁠. We consider four different resistivity profiles: one that is constant, a temperature dependent Spitzer model, and two others that are scaled versions of the latter.

The simulations find that the equilibrium is highly unstable to pressure-driven MHD modes. The behavior is qualitatively unchanged for the different resistivity profiles. This, together with the fast onset, suggests that the equilibrium is ideally MHD unstable. The magnetic field rapidly develops a chaotic layer near the plasma edge that continues to grow in time, until the magnetic field ceases to be confining.

The destabilization of a large number of toroidal Fourier modes is related to the N-fold rotational symmetry of the equilibrium (N = 2). This discrete symmetry leads to coupling of Fourier modes in the toroidal direction. For an equilibrium with N = 2, an unstable even-numbered toroidal Fourier mode couples, in principle, to all other Fourier modes with even toroidal mode number. The same is true of an odd-numbered toroidal Fourier mode. In this case, the single unstable eigenmode consisted of many even-numbered toroidal Fourier modes.

Finally, we observed some notable differences in the nonlinear behavior of the high-n Fourier modes between the temperature dependent and independent resistivity profiles. In the case with constant resistivity, the growth rates plateau before decreasing and tending to zero. This is evidence of the linear growth phase of an MHD instability, followed by nonlinear saturation. On the other hand, the temperature dependent profiles experience a second phase where the growth rate accelerates. This suggests there are strong nonlinear interactions that provide additional drive to the mode.

The observations of this work offer new insight into the criteria that are important for stellarator optimization. Foremost, field-periodicity plays an important role in nonlinear MHD stability, especially if there is a linearly unstable mode. Given the interest in relaxing linear ideal MHD stability as an optimization constraint, this trade-off should be understood. Second, simulations with the (more realistic) Spitzer resistivity profile observe nonlinear interactions that accelerate the growth of unstable toroidal Fourier modes. By contrast, the simulation with a constant resistivity profile predicts nonlinear saturation. This difference is important to consider when developing reduced models for nonlinear MHD in stellarators, especially for application to optimization.

The authors gratefully acknowledge valuable discussions with M. Landreman. This work was supported by Department of Energy Contract No. DE-AC02-09-CH11466.

The authors have no conflicts to disclose.

A. M. Wright: Conceptualization (lead); Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (lead); Project administration (lead); Resources (lead); Software (lead); Validation (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (lead). N. M. Ferraro: Software (supporting); Supervision (lead); Writing – review & editing (supporting).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

The full sequence of Poincaré sections for the case with constant resistivity, shown as a function of time in Fig. 6.