From tokamaks to stellarators: understanding the role of 3D shaping

TL;DR

This paper investigates how three-dimensional shaping influences rotational transform in toroidal magnetic confinement devices, comparing theoretical susceptance matrices with simulations, and exploring the effects of different shaping modes on plasma stability.

Contribution

It introduces a comparison between the susceptance matrix and 3D MHD equilibrium simulations, revealing how specific shaping modes influence rotational transform in stellarators.

Findings

High poloidal mode number shapes affect edge transform minimally.

Modes with m < 2 influence both edge and core rotational transforms.

3D shape non-linearity impacts magnetic symmetry and plasma behavior.

Abstract

In this work, the role which three-dimensional shaping plays in the generation of rotational transform in toroidal magnetically confinement devices is explored. The susceptance matrix as defined by Strand & Houlberg (2001) is presented and compared to simulations of three dimensional MHD equilibria. The dependence of the edge rotational transform on axisymmetric shaping is briefly explored. It is then shown how simple three dimensional shaping can drive rotational transform in a simply shaped stellarator. It is found that modes with high poloidal mode number drive edge transform with little effect on the core, while modes with poloidal mode m < 2 drive both edge and core rotational transforms. This work concludes with a discussion of the non-linearity of the three dimensional shape space as defined in terms of magnetic symmetry.

Equations

$$\vec{F}=q\left(\vec{E}+\vec{v}\times \vec{B}\right)$$

$$\iota =\frac{d\Psi }{d\Phi }.$$

$$\iota =\frac{{\mu }_{0}I}{S_{11}{\Phi }^{\prime }}-\frac{S_{12}}{S_{11}}$$

$${\mu }_0\left(\begin{array}{c}I\\ F\end{array}\right)=\left(\begin{array}{cc}{S}_{11}& S_{12}\\ S_{21}& S_{22}\end{array}\right)\left(\begin{array}{c}d\Psi /d\rho \\ d\Phi /d\rho \end{array}\right)$$

$$\vec{B}=\frac{1}{2\pi \sqrt{g}}\left[\left(\frac{d\Psi }{d\rho }-\frac{d\lambda }{d\zeta }\frac{d\Phi }{d\rho }\right){\vec{e}}_{\theta }+\left(1+\frac{d\lambda }{d\theta }\right)\frac{d\Phi }{d\rho }{\vec{e}}_{\zeta }\right].$$

$$S_{11}=\frac{dV/d\rho }{4{\pi }^2}⟨\frac{g_{\theta \theta }}{g}⟩$$

$$S_{12}=\frac{dV/d\rho }{4{\pi }^2}⟨\frac{g_{\theta \zeta }\left(1+d\lambda /d\theta \right)-g_{\theta \theta }d\lambda /d\zeta }{g}⟩$$

$$S_{21}=\frac{dV/d\rho }{4{\pi }^2}⟨\frac{g_{\zeta \zeta }\left(1+d\lambda /d\theta \right)-g_{\zeta \theta }d\lambda /d\zeta }{g}⟩$$

$$S_{22}=\frac{dV/d\rho }{4{\pi }^2}⟨\frac{g_{\zeta \theta }}{g}⟩$$

$$R\left(\rho ,{\theta }^{\ast },\zeta \right)=\sum _{n=-N}^N\sum _{m=0}^{M}R\left(\rho \right)\cos \left(m{\theta }^{\ast }-n{N}_{FP}\zeta \right)$$

$$Z\left(\rho ,{\theta }^{\ast },\zeta \right)=\sum _{n=-N}^N\sum _{m=0}^{M}Z\left(\rho \right)\sin \left(m{\theta }^{\ast }-n{N}_{FP}\zeta \right)$$

$$R=\frac{\delta P}{\delta I}=\frac{{\mu }_0/{\Phi }^{\prime }-\partial S_{12}/\partial I}{{\mu }_{0}I/{\Phi }^{\prime }-S_{12}}$$

Full text

1 Introduction

The goal of magnetic confinement nuclear fusion is to use magnetic fields to confine a hot plasma long enough so that a burning plasma state can be reached. While many magnetic configurations have been considered, the topological torus shows the greatest maturity. In this class of device two subclasses exist, those which are axisymmetric and those which are non-axisymmetric. In the axisymmetric branch, exist the tokamak and field reversed configuration. While in the non-axisymmetric branch one finds stellarators and heliotrons. Spanning the space are the reverse field pinch devices which can exhibit an axisymmetric nature, or bifurcate into a state with non-axisymmetric plasma cores. In this paper, the role of three-dimensional shaping in confinement is explored.

2 Theory

Consider a charged particle in a uniform magnetic field. The relevant force acting on the particle is the Lorentz force

[TABLE]

where $q$ is the particles charge, $\vec{E}$ the electric field, $\vec{v}$ its velocity, and $\vec{B}$ is the magnetic field. If the particle possesses any component of velocity perpendicular to the magnetic field it will begin to exhibit cyclotron motion. One may generate such a field using a solenoid, however collisions between particles will tend to equalize the components of the velocity parallel and perpendicular to the magnetic field. Thus particles will tend to stream along the magnetic field line at finite temperatures (ignoring an additional hierarchy of complicated physics). This implies that we must in some way ‘cap’ the ends of this magnetic bottle.

One concept for ‘capping’ the bottle is to use a magnetic mirror force to reflect the particles back along the field line. However, this is only a partial solution as the fastest particle will still tend to escape. Ultimately cooling the plasma before being able to reach fusion relevant temperatures. The more agreed upon solution is to deform the straight field into a topological torus. While this avoids the problems of end losses, the magnetic field now becomes curved (and a field gradient forms). This curvature of the magnetic field causes particles to drift across magnetic field lines. For particles on the inboard side of the torus, such a drift brings them toward the center of the device. While for particles on the outboard side, the motion is away from the center.

If the magnetic field lines could be curved in such a way that particles spend more of their time in the ‘good’ curvature region than in the ‘bad,’ the net effect of the magnetic field would be to ‘confine’ particles to the hot core region. The ‘Hairy Ball Theorem’ of algebraic topology (Eisenberg & Guy, 1979) tells us that for a torus this is possible. From this point forward we will restrict ourselves to a discussion where nested surfaces of toroidal flux exist everywhere in our magnetic configuration up to some enclosing flux surface. This enclosing surface will be called the last closed magnetic surface (LCMS). At the center of the torus we thus have a single magnetic field line. While this assumption may in fact be artificial, the ramifications of such an assumption will be discussed later in this work.

The ‘rate’ at which magnetic field lines rotate about the magnetic axis is defined by a quantity called the rotational transform, which is defined in terms of poloidal ($\Psi$) and toroidal ($\Phi$) magnetic fluxes (Bateman, 1973; Strand & Houlberg, 2001)

[TABLE]

So, put another way, the act of confining charged particles with magnetic fields in a toroidal device is the act of generating rotational transform.

Rotational transform can be generated by non-axisymmetric shaping or toroidal current as defined by

[TABLE]

where $\mu_{0}$ is the permeability of free space, $I$ the toroidal current, $\Phi^{\prime}=d\Phi/d\rho$, and $S$ are the components of the susceptance matrix. In general the suseptance matrix is defined as:

[TABLE]

where $F$ is the poloidal current and $\rho$ is a radial coordinate. Assuming a nested flux surface representation and a magnetic stream function $\lambda$ we can redefine our poloidal angle ($\theta^{*}=\theta+\lambda\left(\rho,\theta,\zeta\right)$), then allowing one to write the magnetic field in the form

[TABLE]

where $\sqrt{g}$ is the coordinate Jacobian. Here we note that our choice of magnetic field representation is exactly that of the VMEC equilibrium code (Hirshman & Whitson, 1983), allowing us to write the components of the susceptance matrix in terms of geometric quantities

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $g_{kk}$ are the components to the coordinate metric tensor and $V$ is the volume enclosed in a flux surface. Clearly the $S_{11}$ and $S_{22}$ coefficients depend purely on the geometry of the flux surfaces. Figure 1 depicts the values of the $S_{11}$ and $S_{12}$ coefficients for a few common MHD equilibria. Note that at the edge ($x=1$) the W7-X device has a rather small value of $S_{11}$. This suggests a strong sensitivity of the edge rotational transform to enclosed toroidal current (an already well documented issue) (Bozhenkov et al., 2013). The relatively large values of $S_{11}$ in the NCSX and QHS configurations suggest a robustness of the configurations to current (an already documented feature) (Bader et al., 2017).

3 The dependence of rotational transform on flux surface geometry

Equation 3 tells us that there are two ways of generating transform: toroidal current and three-dimensional shaping. To explore this statement we employ the VMEC equilibrium code to calculate various equilibria and assess the impact of magnetic geometry on rotational transform. To do this we first examine the effect of the $g_{\theta\theta}$ term in axisymmetry then move onto non-axisymmetric perturbations.

As VMEC will be our tool it is worthwhile to review the nature of the flux surface representation used in the code. For stellarator symmetric equilibrium the flux surfaces are parametrized by functions which trace out the magnetic surfaces, namely

[TABLE]

[TABLE]

where $N$ is the toroidal truncation of modes, $M$ the poloidal truncation of modes, and $N_{FP}$ the magnetic field periodicity. The $\lambda\left(\rho,\theta^{*},\zeta\right)$ function has odd parity, thus a similar form to equation 11.

3.1 Axisymmetric equilibria

For axisymmetric equilibria the $S_{12}$ component of the susceptance matrix vanishes, leaving only the $S_{11}$ component (and toroidal current) to generate rotational transform. An aspect ratio $10$ circular cross section equilibrium (at zero beta but finite toroidal current) serves as a basis for examining the relationship between $S_{11}$ and rotational transform. It is then a simple matter to elliptically deform the equilibria by adjusting the $n=0,m=1$ mode of the $Z$ harmonics. Figure 2 depicts the dependence of edge transform on the $S_{11}$ component of the susceptance matrix as the equilibrium is deformed. In this plot the equilibria with larger deformation are located on the left of the graph. It is important to remember that rotational transform is the inverse of the safety factor ($\iota=1/q$). The result tells us that as we add ellipticity to the axisymmetric equilibrium the value of $q$ increases (holding all else constant).

If we consider the effect on the iota profile, this result holds for other modes which increase the $g_{\theta\theta}$ metric. The toroidal current to zeroth order acts as a rigid shift to the rotational transform (ignoring profile effects). While the boundary perturbations have more nuanced effects on the rotational transform profile. The $m=1$ modes (which generate ellipticity) tend to penetrate to the core, meanwhile the higher order modes tend to decay into the core with higher orders falling off faster. This makes it possible to tailor the iota profile through shaping and current. This implies that multiple combinations of current profile and axisymmetric shaping may be employed to generate a given transform profile.

3.2 Non-axisymmetric equilibria

When the toroidal symmetry of the device is broken the $S_{12}$ term becomes non-zero allowing rotational transform to be produced even when the enclosed toroidal current becomes zero. Such is the design concept surrounding stellarators and heliotrons. Here the ratio $\iota_{vac}=-S_{12}/S_{11}$ defines the ‘vacuum’ transform. It should now be clear that by maximizing the $S_{11}$ component of the susceptance, a robustness of the iota profile to currents (be they ohmic, bootstrap, beam driven or electron-cyclotron drive) can be built into an equilibrium. However, this also suggests that the $S_{12}$ component be maximized as well to generate the vacuum transform.

A better understanding of how boundary shapes drive vacuum transform in stellarators can be seen in figure 3. We begin by perturbing the axisymmetric equilibria (toroidal current set to zero) with $n=1,m=1$ (upper plots) and $n=2,m=2$ perturbations (lower plots) and examine the resulting rotational transform profiles. The $n=1,m=1$ perturbation deforms the equilibria into a rotating ellipse. The resulting rotational transform profile posses reversed shear and modifies both the edge and axis values of the profile. The $n=2,m=2$ perturbation creates a rotating triangular feature to the equilibrium, becoming slightly scalloped at the largest amplitude explored. It is worthwhile to note that the $n=2$ nature implies that the half-field period plane has the same cross section as the other symmetry planes. The $m=2$ nature of this perturbation results in a modification to the rotational transform at the edge but not in the core. We note that the sign of these transform profiles is a function of the choice of sign in the toroidal modes.

While these two examples appear trivial, they already show that three dimensional shaping can be used to tailor the vacuum rotational transform of a stellarator. Specifically the low order modes ($m<2$) should provide an ability to tailor the transform on axis. Meanwhile the higher order modes allow tailoring of the edge transform. Thus the design of a stellarator becomes a the task of combining these modes to produce an equilibrium with the desired properties.

Stellarators usually use a multitude of modes to drive rotational transform. In figure 4, a set of boundary deformations is explored where the relative amplitudes of the modes are equal. We see that modes with $m=1$ nature penetrate to the core region while modes of higher m are damped. By changing the sign of the modes the sign of their contribution to the rotational transform can also be changed. We see that the $n=1,m=2$ mode drive a similar level of core transform as that of the $n=1,m=1$ mode but significantly more edge transform. Simulations indicate that it is possible to combine the $n=1,m=1$ boundary modes and $n=2,m=2$ boundary modes and recompute the equilibria. The resulting equilibria has a transform almost identically matching that found though summation of the two individual transforms. This result may not hold should the magnetic axis become sufficiently non-axisymmetric.

4 Discussion

This work has begun to explore the role of three dimensional shaping on the generation (and modification) of rotational transform in ideal magnetohydrodynamic toroidal equilibria. While not exhaustive, this work shows that modes with poloidal modes less than $2$ are the primary means by which transform on axis can be driven. The dependence of core transform on edge poloidal perturbations appears weak at best, with the $n=1,m=1$ mode driving nearly the same amount of transform as the $n=2,m=1$ mode for similar perturbation amplitudes. At the edge of the equilibrium modes with higher poloidal mode number appear to drive larger amounts of transform for the same level of boundary displacement. However, this suggests that the same (or similar) rotational transform profiles can be achieved through a multitude of boundary modes. Numerical simulations combining different modes supports the idea of linear combinations of modes to produce a profile which is the the sum of the profiles driven by different modes.

This work suggests a new target for stellarator optimization. Specifically if a stellarator was optimized in such a way to maximize the $S_{11}$ coefficient, the resulting configuration would become robust to current changes. This would provide rotational transform robustness to bootstrap currents. This is desirable in stellarators which employ island divertors (such as Wendelstein 7-X). These devices rely on an edge island chain of known position and width. Thus these devices would benefit from edge iota profile rigidity. One metric could be defined by the transform profile rigidity as

[TABLE]

where we have defined low values of rigidity ($R$) to indicate a lack of sensitivity to toroidal currents. The variable $P$ is used to represent the normalized change in rotational transform for a given unit current. It should be noted that the partial derivative ($\partial S_{12}/\partial I$) can become quite complex and as such we leave leave further analysis to future work.

It is important to note that in this work, we have ignored the possibility of islands and stochastic regions at rational surfaces. Moreover, we have completely neglected the recent results by Loizu et al. (2015) regarding the necessity of non-smooth rotational transform profiles across rational surfaces. If we imagine a plasma model which allows the islands to open, then we expect the island width to scale as the resonant harmonic. While modifying the profile in the neighborhood of the island, overall the equilibrium profile is mostly unperturbed (excluding pathological examples). We note that there clearly exists a freedom in approximating a rotational transform profile using many different modes. This suggest that one could design a stellarator by first choosing a rotational transform profile with known resonances, then explicitly setting the resonant component to zero and using the non-resonant boundary deformations to create the transform.

If we consider inclusion of shielding currents, then our equilibrium rotational transform profile becomes discontinuous at each resonant rational surface. These discontinuities scale as the amplitude of the resonant harmonic and should be small in most cases. The resulting equilibria are only modified in the vicinity of these surfaces and should not impact the generic features of the equilibrium (Lazerson et al., 2016). One should note that such shielding currents require toroidal sheet currents and as such are inconsistent with the zero beta, zero current equilibrium limit.

Before concluding this discussion the authors would like to take the time to point out that generation of transform alone is not enough to confine the particles in a nuclear fusion device. To this end a great deal of work (mostly associated with stellarators) has gone into the development of optimized three-dimensional equilibria (Mynick et al., 1982; Spong et al., 1998; Mynick, 2006). Here the optimization usually focuses on minimizing the radial transport of particles through magnetic symmetries. These symmetries include helically symmetric (the HSX device) (Boozer, 1995), quasi-oniginatiy (W7-X) (Hall & McNamara, 1975), and quasi-axisymmetry (NCSX) (Garabedian, 1996). NCSX being relevant to this paper as the tokamak falls into the quasi-axisymmetric (QAS) category (axisymmetry being the exact realization of quasi-axisymmetry). Of course any perturbation to a tokamak which breaks toroidal symmetry takes the equilibrium away from it’s ideal axisymmetric state. However, the situation in practice is quite more complicated. In figure 5, we have mapped the n-dimensional space around the NCSX equilibrium by defining two vectors and evaluating equilibria in the subspace defined by these two vectors. In one direction (along x-axis) we smoothly deform the NCSX equilibria into an axisymmetric tokamak equilibria. In the other direction, we smoothly deform the boundary toward a turbulent transport optimized version of NCSX (Mynick et al., 2010). The amplitude of the color map indicates how non-axisymmetric the magnetic equilibria are, a sum over $|B|$ modes with toroidal mode numbers other than zero (normalized to the amplitude of those modes with $n=0$.) This database was generated from $\sim 10000$ VMEC runs produced by the STELLOPT code.

The purpose of showing this seemingly separate work is to demonstrate that in moving from an axisymmetric tokamak to the quasi-axisymmetric NCSX equilibrium the non-axisymmetry increases and then decreases. And while NCSX is clearly not axisymmetric, it is predicted to have transport levels approaching that of a tokamak. It should now be clear that the selection of modes to generate transform in stellarators is much more complicated than the analysis presented here. Still this work acts as a starting point, and future works will explore the dependency of radial transport, stability, and bootstrap current in future works.

5 Citations and references

The authors would like to acknowledge fruitful discussions with J. Loizu, S. Hudson, A. Boozer, H. Mynick, N. Pomphrey, S. Prager, and D. Gates. Furthermore, the authors would like to thank S. Hirshman, S. Seal, and M. Cianciosa for access to the VMEC and PARVMEC codes.