2.1. Neoclassical and turbulent transport

In stellarators, particle and energy transport are caused by both neoclassical transport and turbulence. The diffusion coefficient associated with the latter is, according to local gyrokinetic theory, of the order of the gyro-Bohm value (Hagan & Frieman Reference Hagan and Frieman1986; Connor Reference Connor1988):

(2.1)\begin{equation} D_{\rm gB} \sim \rho_{{\ast} i}^2 v_{\rm Ti} a, \end{equation}

where $v_{{\rm Ti}} = \sqrt {2 T_i/m_i}$ is the ion thermal velocity and $\rho _{\ast i} = \rho _i/a$ the gyroradius normalised to the minor radius $a$.Footnote ² The magnitude of the neoclassical transport depends sensitively on the geometry of the magnetic field and the collisionality. The electrons are usually in the ‘$1/\nu$ collisionality regime’, in which the diffusion coefficient scales as (Galeev et al. Reference Galeev, Sagdeev, Furth and Rosenbluth1969; Galeev & Sagdeev Reference Galeev and Sagdeev1979; Ho & Kulsrud Reference Ho and Kulsrud1987; Beidler et al. Reference Beidler2011)

(2.2)\begin{equation} D^{1/\nu}_e \sim \frac{\epsilon_{\rm eff}^{3/2} v_{\rm de}^2}{\nu_e}, \end{equation}

where $v_{\rm de} \sim \tau \rho _{\ast i} v_{\rm Ti}$ denotes the electron drift velocity and $\nu _e$ the electron collision frequency, which is related to that of the ions by

(2.3)\begin{equation} \nu_e \sim \frac{M^{1/2} \nu_i}{\tau^{3/2}}, \end{equation}

with $\tau = T_e/T_i$ and $M = m_i / m_e$ (Helander & Sigmar Reference Helander and Sigmar2005). The so-called effective magnetic ripple $\epsilon _{\rm eff}$ depends on the geometry of the magnetic field and measures the net effect of the radial magnetic drift, integrated over all trapped orbits, for this type of transport (Nemov et al. Reference Nemov, Kasilov, Kernbichler and Heyn1999). It vanishes in an exactly omnigenous field, where, by definition, drift surfaces and flux surfaces coincide. An important goal of stellarator optimisation is to make the ratio of neoclassical to turbulent transport smaller than unity:

(2.4)\begin{equation} \frac{D^{1/\nu}_e}{D_{\rm gB}} \sim \frac{\epsilon_{\rm eff}^{3/2}\tau^{7/2}}{M^{1/2} \nu_{{\ast} i}} \lesssim 1, \end{equation}

where $\nu _\ast = \nu _i a / v_{\rm Ti}$. In practice, $\epsilon _{\rm eff}$ therefore needs to be a few per cent or smaller for relevant values of the plasma parameters. That this goal is achievable in a large, low-collisionality stellarator has recently been demonstrated by experiments in W7-X (Beidler et al. Reference Beidler, Smith, Alonso, Andreeva, Baldzuhn, Beurskens, Borchardt, Bozhenkov, Brunner and Damm2021).

The ions are usually in a different collisionality regime, the so-called $\sqrt {\nu }$ regime, where the diffusion coefficient scales as (Galeev et al. Reference Galeev, Sagdeev, Furth and Rosenbluth1969; Ho & Kulsrud Reference Ho and Kulsrud1987; Calvo et al. Reference Calvo, Parra, Velasco and Alonso2017)

(2.5)\begin{equation} D_i^{\sqrt{\nu}} \sim c_{\sqrt{\nu}} \frac{\nu_i^{1/2} \rho_{{\ast} i}^2 v_{\rm Ti}^2}{\omega_E^{3/2}} \sqrt{\ln \left( \frac{\omega_E}{\nu_i} \right)}. \end{equation}

Here, the coefficient $c_{\sqrt {\nu }}$ encapsulates the effect of the magnetic field geometry on this type of transport, and is given by $c_{\sqrt {\nu }} = (b_{10}/\epsilon _t)^2$ in a classical stellarator, where $b_{10}$ measures the lowest-order poloidal harmonic in the variation of the magnetic field strength and $\epsilon _t$ is the so-called toroidal ripple. We refer the reader to Beidler et al. (Reference Beidler2011) for more details. The denominator contains the ${\boldsymbol {E}} \times {\boldsymbol {B}}$ rotation frequency $\omega _E = |E_r |/ (rB)$, where $E_r$ denotes the radial electric field, $r$ the minor radius and $B$ the magnetic field. The logarithmic factor, which was relatively recently found by Calvo et al. (Reference Calvo, Parra, Velasco and Alonso2017), is ignored in the following since it is rarely very different from unity.

Depending on the magnetic-field geometry and the collisionality, regimes other than (2.5) are also possible (Mynick Reference Mynick2006). In practice, these are not well separated from one another, and an accurate calculation of the transport must therefore be done numerically. A lower bound on the transport is given by the banana-regime coefficient associated with an orbit width proportional to the gyroradius:

(2.6)\begin{equation} D_i^{\rm ban} \sim c_{\rm ban} \nu_i \rho_i^2, \end{equation}

which has a scaling similar to that of classical (Braginskii) transport, tokamak banana transport and Pfirsch–Schlüter transport (Helander & Sigmar Reference Helander and Sigmar2005). (A similar bound applies to the electrons, which is, however, extremely small.) An exactly omnigenous stellarator would exhibit such transport (Boozer Reference Boozer1983; Helander & Nührenberg Reference Helander and Nührenberg2009) and is indeed seen numerically in extremely well-optimised configurations (Goodman et al. Reference Goodman, Camacho Mata, Henneberg, Jorge, Landreman, Plunk, Smith, Mackenbach, Beidler and Helander2023).

2.2. Ambipolarity

An important qualitative difference between neoclassical and turbulent transport is the issue of intrinsic ambipolarity. Particle transport is said to be intrinsically (or automatically) ambipolar if the radial electric current, i.e. the sum over all species of the charge-weighted particle fluxes, vanishes for any value of the radial electric field. Stellarator neoclassical transport is never intrinsically ambipolar unless the magnetic field is quasi-symmetric (Helander & Simakov Reference Helander and Simakov2008; Helander Reference Helander2014). Only for certain specific values of the radial electric field does the radial current vanish.

In contrast, turbulent transport is intrinsically ambipolar to leading order in the gyroradius expansion according to standard gyrokinetic theory (Sugama et al. Reference Sugama, Okamoto, Horton and Wakatani1996; Parra & Catto Reference Parra and Catto2009; Helander Reference Helander2014). Any non-ambipolar contribution to the turbulent transport is thus expected to be associated with a diffusion coefficient of at most

(2.7)\begin{equation} D_{\rm turb}^{\rm na} \lesssim \rho_{{\ast} i} D_{\rm gB}. \end{equation}

As a result, even if a stellarator has been optimised well enough that (2.4) holds, the non-ambipolar part of the transport will usually be dominated by the neoclassical contribution. Only if the geometry of the magnetic field has been optimised to the extent that

(2.8)\begin{equation} {D_e^{1/\nu}}\lesssim \rho_\ast {D_{\rm gB}} \end{equation}

does the turbulent transport significantly affect the ambipolarity condition. Such extreme optimisation is, however, unnecessary for energy confinement and is rarely undertaken in practice. As a result, although most of the energy transport in large stellarators such as LHD and W7-X (Beidler et al. Reference Beidler, Smith, Alonso, Andreeva, Baldzuhn, Beurskens, Borchardt, Bozhenkov, Brunner and Damm2021) is caused by plasma turbulence, the radial electric field is nevertheless expected to be determined by neoclassical transport to leading order in the smallness of the gyroradius (Sugama et al. Reference Sugama, Okamoto, Horton and Wakatani1996; Helander & Simakov Reference Helander and Simakov2008; Helander Reference Helander2014). An exception occurs on small radial scales, where short-wavelength zonal flows can be excited by turbulence (Helander & Simakov Reference Helander and Simakov2008).

2.3. Electron and ion roots

In non-quasi-symmetric stellarators, then, ambipolarity is not automatic but is realised only for one or several particular values of the radial electric field. The dependence of the radial current on $E_r$ is nonlinear, and it has long been known that it can have up to three roots (Mynick & Hitchon Reference Mynick and Hitchon1983; Shaing Reference Shaing1984; Hastings, Houlberg & Shaing Reference Hastings, Houlberg and Shaing1985; Hastings Reference Hastings1986). If the ion and electron temperatures are comparable, there is usually only the so-called ‘ion root’ associated with an inward-pointing electric field, $E_r < 0$. However, if the electrons are relatively hot, $\tau > 1$, an ‘electron root’ with $E_r >0$ is possible. A third, intermediate, root is unstable and therefore expected to be irrelevant.Footnote ³

These theoretical expectations are borne out in practice. Stellarator experiments demonstrate that the radial electric field is usually negative but switches sign in the centre of the plasma if strong heating is selectively applied to the electrons. In a range of several very different stellarators (CHS, LHD, TJ-II, W7-AS), localised electron cyclotron resonance heating (ECRH) resulted in a regime of improved confinement, the so-called ‘core electron-root confinement regime’, associated with postive $E_r$ (Yokoyama et al. Reference Yokoyama2006). Above a certain ECRH power, the electron temperature profile became very steep, enabling a relatively high electron temperature to be reached in the core. Similar plasmas have more recently been created in W7-X (Geiger et al. Reference Geiger2019). In HSX, however, only the ion root has been observed although the electrons are usually considerably hotter than the ions, but this stellarator is quasi-helically symmetric and should therefore behave differently due to intrinsic ambipolarity of the neoclassical transport.

Theoretically, the ion root appears because the electrons are the ‘rate-controlling’ species (Ho & Kulsrud Reference Ho and Kulsrud1987). In the absence of an electric field, both the electrons and the ions would be in their respective $1/\nu$ regimes, but the ion transport would then be much larger than the electron transport, violating ambipolarity. A radial electric field must therefore arise to reduce the ion transport to the electron level, which thus sets the overall transport. The electric field makes the ion transport follow the $\sqrt {\nu }$ scaling (2.5) and adjusts to a value where the latter becomes comparable to (2.2).

For a hydrogen plasma with the electrons in the $1/\nu$ regime and the ions in the $\sqrt {\nu }$ regime, the ambipolarity condition is

(2.9)\begin{equation} D_e^{1/\nu} \left( \frac{{\rm d}\ln n}{{\rm d}r} + \frac{e E_r}{T_e} + \delta_e \frac{{\rm d} \ln T_e}{{\rm d}r} \right) = D_i^{\sqrt{\nu}} \left( \frac{{\rm d}\ln n}{{\rm d}r} - \frac{e E_r}{T_i} + \delta_i \frac{{\rm d} \ln T_i}{{\rm d}r} \right), \end{equation}

with $\delta _e = 7/2$ and $\delta _i = 5/4$ (Beidler et al. Reference Beidler, Smith, Alonso, Andreeva, Baldzuhn, Beurskens, Borchardt, Bozhenkov, Brunner and Damm2021), and implies a radial electric field given by

(2.10)\begin{equation} \frac{e E_r}{T_e} = \frac{1 + \eta_i \delta_i - (1 + \eta_e \delta_e) x}{x + \tau} \cdot \frac{{\rm d}\ln n}{{\rm d}r} , \end{equation}

where $x =D_e^{1/\nu } / D_i^{\sqrt {\nu }}$ and $\eta _\sigma = {\rm d} \ln T_\sigma / {\rm d} \ln n$. For a conventional density profile with $n'(r) < 0$, a positive radial electric field is only realised if

(2.11)\begin{equation} x > \frac{1 + \eta_i \delta_i }{1 + \eta_e \delta_e}. \end{equation}

In order to derive a crude scaling law for plasma parameters admitting an electron root, we take the right-hand side of this equation to be of order unity and radial scale lengths to be comparable to $a$. The electric field then becomes of order

(2.12)\begin{equation} E_r \sim \frac{T_i r}{e a^2}, \end{equation}

so that $\omega _E \sim \rho _{\ast i} v_{{\rm Ti}}/a$. If the logarithmic factor in (2.5) is neglected and $r \sim a$, the condition $x \gtrsim 1$ leads to the criterion

(2.13)\begin{equation} \tau \gtrsim M^{1/7} c_{\sqrt{\nu}}^{2/7} \left( \frac{\nu_{{\ast} i}}{\epsilon_{\rm eff} \rho_{{\ast} i}} \right)^{3/7} \equiv \tau_\ast, \end{equation}

which should approximately predict the appearance of an electron root. The scalings in this criterion have recently been found to agree with global gyrokinetic simulations (Kuczynski et al. Reference Kuczynski, Kleiber, Smith, Beidler, Borchardt, Geiger and Helander2024). Several interesting conclusions can be drawn from this result.

Firstly, it explains why an electron root is at all possible despite the large mass difference between electrons and ions, $M \sim 10^3$. The latter is responsible for the fact that ions have larger neoclassical transport than electrons, and one may be forgiven for wondering whether the reverse can ever be possible. However, the ratio of electron to ion temperature at which this happens is, according to (2.13), only proportional to $M^{1/7}$, which is about 3 in a hydrogen plasma. The electrons therefore need not be very much hotter than the ions to enable the appearance of an electron root.

Secondly, it is interesting to note that, although the $\rho _{\ast i}$ scaling in (2.13) is unfavourable for large stellarators, such devices have relatively good confinement, enabling higher temperatures and lower collisionalities. Thanks to the favourable scaling of (2.13) with $\nu _{\ast i}$, which enters as $\nu _{\ast i} / \rho _{\ast i} \propto a^2 n /T_i^{5/2}$, the ratio of electron to ion temperature $\tau _\ast$ necessary for an electron root is therefore not necessarily much higher in a stellarator reactor than in present experiments. However, the ion and electron temperatures will almost be equal in a reactor, since the energy confinement time is likely to exceed the Braginskii temperature-equilibration time, making it difficult to attain $\tau > 1$.

Finally, and most importantly, although electron roots are usually only observed when the electrons are hotter than the ions, (2.13) suggests that an electron root should be possible even if $\tau = 1$. If a stellarator could be optimised to have a sufficiently small ratio $c_{\sqrt {\nu }} / \epsilon _{\rm eff}^{3/2}$, an electron root should be attainable in temperature-equilibrated plasmas.

These conclusions do not change much if the ion transport is instead given by the lower bound of (2.6). Then $D_e > D_i$ when

(2.14)\begin{equation} \tau > \frac{c_{\rm ban}^{2/7}}{\epsilon_{\rm eff}^{3/7}} \, M^{1/7} \nu_{{\ast} i}^{4/7}, \end{equation}

which shares many features with (2.13). Only the scaling with $\rho _{\ast i}$ is significantly different.

3.1. Complete omnigenity

Hall & McNamara (Reference Hall and McNamara1975) introduced the notion of magnetic fields in which the net radial drift of all trapped particles vanishes, and called such fields omnigenous. Cary & Shasharina (Reference Cary and Shasharina1997) derived a basic criterion for exact omnigenity, which we now re-derive. If the magnetic field is written as $\boldsymbol {B} = \boldsymbol {\nabla } \psi \times \boldsymbol {\nabla } \alpha$, with $\psi$ denoting the toroidal magnetic flux divided by $2 {\rm \pi}$ and $\alpha = \theta - \iota \varphi$ in Boozer coordinates, the net drift of a magnetically trapped particle in the $\psi$ direction is equal to (Helander Reference Helander2014)

(3.1)\begin{equation} \Delta \psi = \frac{1}{q} \frac{\partial \mathcal{J}}{\partial \alpha} , \end{equation}

where

(3.2)\begin{equation} \mathcal{J}(\psi,\alpha,y) \equiv \int mv_\parallel \,{\rm d}l = \sqrt{2m\mathcal{E}} \int\sqrt{1-B(\psi,\alpha,l)/y}\,{\rm d}l \end{equation}

is the so-called second (or parallel) adiabatic invariant. Here $m$ denotes the mass, $q$ the charge and $v_\| = \pm v \sqrt {1- B/y}$ the parallel velocity of the particle, where $y = \mathcal {E} / \mu$ is the ratio of kinetic energy ${\mathcal {E}} = mv^2/2$ to magnetic moment $\mu = m v_\perp ^2/2B$, both of which are constant for a particle moving within an equipotential magnetic surface. The integral is taken along the magnetic field between two consecutive bounce points of equal field strength $y$, with the arc length along the field line between these points denoted by $l$.

It follows that omnigenity is equivalent to the condition

(3.3)\begin{equation} \frac{\partial}{\partial \alpha} \int \sqrt{y-B(\psi,\alpha,l)} \, {\rm d}l = 0, \end{equation}

for all values of $y \in [B_{\rm min}(\psi ),B_{\rm max}(\psi )]$ between the minimum and maximum field strengths on the magnetic surface in question. Differentiating this equation with respect to $y$ and changing the integration variable from $l$ to $B$ gives

(3.4)\begin{equation} \frac{\partial}{\partial \alpha} \int_{B_{\rm min}}^y F(\psi,\alpha,B) \frac{{\rm d}B}{\sqrt{y-B}} = 0, \end{equation}

where the function $F$ is defined by

(3.5)\begin{equation} F(\psi,\alpha,B) = \sum \left| \frac{\partial B}{\partial l} \right|^{{-}1}, \end{equation}

and the sum is taken over the two bounce points. Equation (3.4) is an Abel integral equation,

(3.6)\begin{equation} \int_a^y f(x) \frac{{\rm d}\kern0.7pt x}{\sqrt{y-x}} = g(y), \end{equation}

whose general solution is

(3.7)\begin{equation} f(y) = \frac{1}{\rm \pi} \frac{{\rm d}}{{\rm d}y} \int_a^y g(x) \frac{{\rm d}\kern0.7pt x}{\sqrt{y-x}}. \end{equation}

Since the right-hand side of (3.3) vanishes, we must therefore have

(3.8)\begin{equation} \frac{\partial F}{\partial \alpha} = 0, \end{equation}

which is the omnigenity criterion derived by Cary & Shasharina (Reference Cary and Shasharina1997). As they note, it follows that, for any function $h(\psi,B)$,

(3.9)\begin{equation} \frac{\partial}{\partial \alpha} \int_{B(\psi,\alpha,l) \leqslant y} h[\psi,B(\psi,\alpha,l)] \,{\rm d}l = 0, \end{equation}

for all $y \in [B_{\rm min}(\psi ), B_{\rm max}(\psi )]$. For instance, the choice $h=1$ implies that the distance along the field between two consecutive bounce points with field strength $y$ is independent of $\alpha$. In other words, the distance between bounce points is the same for all field lines on the same flux surface. This result is a powerful tool for proving theorems about plasmas confined by omnigenous magnetic fields (Helander & Nührenberg Reference Helander and Nührenberg2009; Helander, Geiger & Maaßberg Reference Helander, Geiger and Maaßberg2011; Landreman & Catto Reference Landreman and Catto2012).

3.2. Partial omnigenity

In a perfectly omnigenous field, there is neither $1/\nu$ transport nor $\sqrt {\nu }$ transport, i.e. $\epsilon _{\rm eff} = c_{\sqrt {\nu }} = 0$, and the neoclassical transport will therefore be small and tokamak-like. In any actual stellarator, however, the net radial drift never vanishes for all orbits, and there will be some enhancement of the neoclassical transport above the axisymmetric base level. It is therefore of interest to consider how the Cary–Shasharina theorem (3.8) is modified if the net radial drift (3.1) vanishes for some, but not all, trapped orbits.

The simplest case occurs if deeply trapped particles are omnigenous but less deeply trapped ones are not.Footnote ⁴ The integral equation (3.3) is thus taken to hold for all values of $y$ in the range $B_{\rm min} \leqslant y \leqslant B_1$, where $B_{\rm min} \leqslant B_1 < B_{\rm max}$ and $B_{\rm min}$ is independent of $\alpha$ due to omnigenity (Cary & Shasharina Reference Cary and Shasharina1997; Landreman & Catto Reference Landreman and Catto2012; Helander Reference Helander2014). The solution is then the same as before, (3.8), but it is only valid for $y \in [B_{\rm min},B_1]$ and not for $y > B_1$. The reason is obvious: the deeply trapped particles (i.e. those with $y < B_1$) only sample the magnetic field in the omnigenous part of the trapping well, and are thus unaffected by the lack of omnigenity for more shallowly trapped particles (i.e. those with $y > B_1$).

More interesting, and much more relevant to our present considerations, is the opposite case, in which shallowly trapped particles are taken to be omnigenous but deeply trapped ones are not.Footnote ⁵ The question is whether particles that are shallowly trapped can be omnigenous although they must pass through regions where more deeply trapped ones, which are not omnigenous, reside. We thus examine whether (3.4) can hold for $y \in [B_1,B_{\rm max}]$ but not for $y < B_1$. In the former region, the following integral equation must then be satisfied:

(3.10)\begin{equation} \int_{B_1}^y \frac{\partial F(\psi,\alpha,B)}{\partial \alpha} \frac{{\rm d}B}{\sqrt{y-B}} ={-} \frac{\partial}{\partial \alpha} \int_{B_{\rm min}(\psi,\alpha)}^{B_1} F(\psi,\alpha,B) \frac{{\rm d}B}{\sqrt{y-B}}, \end{equation}

where we note that $B_{\rm min}$ may depend on $\alpha$ since the deeply trapped particles are not omnigenous. According to (3.7), the solution is

(3.11)\begin{equation} \frac{\partial F(\psi,\alpha,y)}{\partial \alpha} ={-} \frac{1}{\rm \pi} \frac{\partial^2}{ \partial y \partial \alpha} \int_{B_1}^y \frac{{\rm d}\kern0.7pt x}{\sqrt{y-x}} \int_{B_{\rm min}(\psi,\alpha)}^{B_1} F(\psi,\alpha,B) \frac{{\rm d}B}{\sqrt{x-B}}, \end{equation}

where the two integrals can be interchanged and that over $x$ performed:

(3.12)\begin{equation} \int_{B_1}^y \frac{{\rm d}\kern0.7pt x}{\sqrt{(y-x)(x-B)}} = {\rm \pi}-2 \arctan \sqrt{\frac{B_1 - B}{y-B}}, \end{equation}

giving

(3.13)\begin{equation} \frac{\partial F(\psi,\alpha,y)}{\partial \alpha} ={-} \frac{1}{{\rm \pi}\sqrt{y-B_1}} \frac{\partial}{\partial \alpha} \int_{B_{\rm min}(\psi,\alpha)}^{B_1} F(\psi,\alpha,B) \frac{\sqrt{B_1 - B}}{y-B} \,{\rm d}B. \end{equation}

It is straightforward, though slightly tedious, to verify that this function has the property

(3.14)\begin{equation} \int_{B_1}^y \frac{\partial F(\psi,\alpha,B)}{\partial \alpha} \sqrt{y-B} \,{\rm d}B ={-} \frac{\partial }{\partial \alpha} \int_{B_{\rm min}}^{B_1} F(\psi,\alpha,x) \left( \sqrt{y-x} - \sqrt{B_1 - x} \right) {{\rm d}\kern0.7pt x}, \end{equation}

and thus indeed satisfies the integral equation (3.4), but it does not satisfy (3.3) since with this solution

(3.15)\begin{equation} \frac{\partial }{\partial \alpha} \int_{B_{\rm min}}^y F(\psi,\alpha,B) \sqrt{y-B} \,{\rm d}B = \frac{\partial }{\partial \alpha} \int_{B_{\rm min}}^{B_1} F(\psi,\alpha,x) \sqrt{B_1 - x} \, {{\rm d}\kern0.7pt x}. \end{equation}

The right-hand side of this equation vanishes only if particles with $y = B_1$ are omnigenous, which, however, will not generally be the case in a field that is non-omnigenous for $B(l) < B_1$. Unsurprisingly, if orbits with $y < B_1$ are not omnigenous, those with $y = B_1$ will in general not be so either.

However, it is possible that orbits with slightly larger values of the bounce-point field strength are omnigenous. Indeed, given the dependence of the magnetic field strength on $(\psi,\alpha,l)$ for values up to $B_1$, it is possible to extend $B(\psi,\alpha,l)$ to values $[B_1,B_{\rm max}(\psi )]$ in such a way that all particle orbits with $y \geqslant B_2$ are omnigenous, where $B_2$ is an arbitrarily chosen field strength in the interval $(B_1,B_{\rm max})$. To do so, we merely need to define $F(\psi,\alpha,y)$ for $y \in [B_1,B_2]$ in such a way that

(3.16)\begin{equation} \frac{\partial }{\partial \alpha} \int_{B_{\rm min}}^{B_2} F(\psi,\alpha,x) \sqrt{B_2 - x} \, {{\rm d}\kern0.7pt x} =0, \end{equation}

which can, for instance, be achieved by choosing

(3.17)\begin{equation} \frac{\partial F(\psi,\alpha,y)}{\partial \alpha} ={-} \frac{3}{2} \left(B_2-B_1\right)^{{-}3/2} \frac{\partial}{\partial \alpha} \int_{B_{\rm min}}^{B_1} F(\psi,\alpha,x) \sqrt{B_2 - x} \, {{\rm d}\kern0.7pt x} \end{equation}

independently of $y$ throughout the region $y \in [B_1,B_2]$.

We are thus led to the conclusion that it should be possible to devise magnetic fields in which deeply trapped particles are less well confined than shallowly trapped ones. This should improve the neoclassical confinement of ions insofar as they are in the $\sqrt {\nu }$ regime since this type of transport vanishes if orbits close to the trapped–passing boundary are perfectly omnigenous. According to the estimate (2.13), the ratio of electron to ion temperature necessary for an electron root would then be reduced, and there is no reason why it could not be as low as unity or smaller.

5.1. Magnetic configuration

We now proceed to show a stellarator plasma which has been designed to have an electron root. As we have seen, it is possible, in principle, to optimise for easy access to the electron root, but such a device will only be attractive as a power plant if it satisfies a range of other requirements too. For this reason, we adopt the optimisation strategy outlined by Goodman et al. (Reference Goodman, Xanthopoulos, Plunk, Smith, Nührenberg, Beidler, Henneberg, Roberg-Clark, Drevlak and Helander2024), who devised a method by which a stellarator magnetic field can be numerically optimised to be QI while simultaneously satisfying several other essential physics criteria. Here, we make two amendments to this method. First, we add an additional target function, suggested by Kappel, Landreman & Malhotra (Reference Kappel, Landreman and Malhotra2024), as an inequality constraint, to improve the configuration's compatibility with a set of practically feasible filamentary coils. This additional target had little-to-no effect on the physics performance of the resulting configuration.

Second, and more importantly, we make several modifications to the implementation of the QI target function employed by Goodman et al. (Reference Goodman, Xanthopoulos, Plunk, Smith, Nührenberg, Beidler, Henneberg, Roberg-Clark, Drevlak and Helander2024), which, in essence, penalises differences between the adiabatic invariant (3.2) corresponding to some input field, $\mathcal {J}_I$, and that in a constructed, perfectly QI field, $\mathcal {J}_C$. These differences are measured along various field lines on a flux surface, as well as multiple values of $y$. The form of this target is thus given by

(5.1)\begin{equation} f_\textrm{QI} \propto \sum_{\tilde{y}} \sum_i \sum_j^{j_\textrm{max}} \left[ \mathcal{J}_I(\psi,\alpha_i,\tilde{y})^{p_\mathcal{J}} - \mathcal{J}_C(\psi, \alpha_j, \tilde{y})^{p_\mathcal{J}}\right]^2 , \end{equation}

where $p_\mathcal {J} = 1$ and the various values of $y$, $\alpha _i$ and $\alpha _j$ are measured on uniform grids. Through numerical optimisation, $f_\textrm {QI}$ can be minimised on various flux surfaces in an equilibrium, resulting in QI fields. In this work, we make three minor modifications to this target function.

First, since we need not necessarily sample $y$ on a uniform grid, we biased our optimisation in favour of more shallowly trapped particles by sampling more values of $y$ near $B_\textrm {max}$ than $B_\textrm {min}$. The exact spacing of these points was varied throughout the optimisation to ensure that the confinement of deeply trapped particles was not spoiled too significantly. We parametrise this spacing with the term $p_y$, and define the $i$th value of our new $y$ grid as

(5.2)\begin{equation} \tilde{y}_i = B_\textrm{min} + \left(\frac{i}{N_y}(B_\textrm{max} - B_\textrm{min})^{1/p_y}\right)^{p_y}. \end{equation}

By increasing $p_y$, we create more points near $B_\textrm {min}$ than $B_\textrm {max}$. We varied $p_y$ in the range $[0.5, 2]$ throughout this optimisation, as we felt needed.

Second, we need not necessarily compare $\mathcal {J}$ along every field line $\alpha _i$ to every other field line $\alpha _j$ on the same flux surface. For instance, one could consider only comparing neighbouring field lines, thus more closely approximating an expression for $\partial \mathcal {J}/\partial \alpha$. We found that, by increasing the range of field lines along which we compared $\mathcal {J}$, shallowly trapped particles were better confined. This range of field lines is parametrised in $j_\textrm {max}$ in (5.1), and was also varied throughout the optimisation process. At times, $j_\textrm {max}$ was large enough to encompass the entire flux surface, while at other times, it was small enough to see only neighbouring field lines.

Finally, we note that $\mathcal {J}$ is always larger for shallowly trapped particles than for deeply trapped ones, since both the integrand and the range of integration in (3.2) are larger. Hence, another tool to bias the optimisation towards better confined shallowly trapped particles is to increase the parameter $p_\mathcal {J}$ in (5.1), which increases the relative importance of these particles in the optimisation. As with the other new parameters, $p_\mathcal {J}\in [0.5, 2]$ was varied throughout the optimisation, which was carried out with the simsopt optimisation framework (Landreman et al. Reference Landreman, Medasani, Wechsung, Giuliani, Jorge and Zhu2021) using the VMEC magnetohydrodynamic equilibrium code (Hirshman, van RIJ & Merkel Reference Hirshman, van RIJ and Merkel1986).

The result of this optimisation is a QI stellarator configuration, shown in figure 1, with many attractive properties, the majority of which will be discussed in a subsequent publication. For now, we focus on its properties relevant to this work. The configuration has visually excellent QI quality, as evidenced by its poloidally closed $B$ contours shown in figure 1. The coefficient $\epsilon _{\rm eff}$ varies from 0.004 on the magnetic axis to 0.025 at the plasma edge. We also note that, while the quality of the omnigenity of its deeply trapped particles is diminished, this configuration is able to confine the high-energy ions generated from fusion reactions as well as can other highly optimised configurations (Landreman & Paul Reference Landreman and Paul2022; Goodman et al. Reference Goodman, Camacho Mata, Henneberg, Jorge, Landreman, Plunk, Smith, Mackenbach, Beidler and Helander2023, Reference Goodman, Xanthopoulos, Plunk, Smith, Nührenberg, Beidler, Henneberg, Roberg-Clark, Drevlak and Helander2024). In collisionless orbit-following simulations for 0.2 s (which corresponds to a typical slowing-down time), no particles launched at half-radius are lost, if the magnetic field is scaled to a mean on-axis field strength of $5.4\ \textrm {T}$ and the minor radius of the plasma exceeds 60 cm. This is due to the relatively large radial gradient of $\mathcal {J}$, which is known to improve the confinement of these ions (Sánchez et al. Reference Sánchez, Velasco, Calvo and Mulas2023; Velasco et al. Reference Velasco, Calvo, Sánchez and Parra2023).

Figure 1.

(a) The boundary of the configuration found in this work, with field lines shown in black. (b) This configuration's $B$ contours, in Boozer coordinates, on various flux surfaces within this configuration.

For the purposes of the present paper, the most important property of the configuration is the possibility of an electron root for reactor-relevant plasma parameters, to which we now turn our attention.

5.2. Performance under reactor conditions

To model reactor plasmas, the optimised equilibrium was scaled to a volume of $1450\ {\rm m}^3$, yielding a major radius of $R=20.13$ m and a minor radius of $a=1.91$ m. The average magnetic field strength at the plasma axis was taken to be $B_0=7.5$ T. Simulations were performed for this equilibrium with the predictive version of the one-dimensional transport code NTSS (Turkin et al. Reference Turkin, Beidler, Maaßberg, Murakami, Tribaldos and Wakasa2011). A density scan was performed, varying the central densities of deuterium and tritium, but without changing the profile shape. NTSS then solves the coupled energy balance equations for electrons ($\sigma ={\rm e}$), deuterium ($\sigma ={\rm D}$), tritium ($\sigma ={\rm T}$) and helium ash ($\sigma ={\rm He}$), modelling the energy flux densities as the sum of neoclassical and turbulent contributions, $Q_\sigma = Q_\sigma ^{\rm neo} + Q_\sigma ^{\rm turb}$. The neoclassical portion is determined using tabulated values of the mono-energetic transport transport coefficients calculated by the Drift Kinetic Equation Solver (DKES) (van Rij & Hirshman Reference van Rij and Hirshman1989), while the part from turbulence is modelled by the simple formula $Q_\sigma ^{\rm turb} = -n_\sigma \chi ^{\rm turb} ({\rm d}T_\sigma /{\rm d}r)$, with the magnitude of $\chi ^{\rm turb}$ being adjusted during the course of the density scan to yield fusion plasmas producing $P_\alpha = 600$ MW of alpha-particle heating power. The profile of the radial electric field is simultaneously calculated, with the particle flux densities also being modelled as a sum of neoclassical and turbulent contributions: $\varGamma _\sigma = \varGamma _\sigma ^{\rm neo} + \varGamma _\sigma ^{\rm turb}$. The neoclassical portion can again be calculated from the DKES dataset, and the turbulent portion is obtained from $\varGamma _\sigma ^{\rm turb}= - D^{\rm turb} ({\rm d}n_\sigma /{\rm d}r)$, with the additional assumption that $\chi ^{\rm turb}/D^{\rm turb} = 7.5$. A differential equation is solved for the radial electric field to ensure a unique value of $E_r$ in parts of the plasma where multiple roots of the ambipolarity constraint exist. This equation serves to minimise the generalised heat production rate, thereby introducing thermodynamic considerations into the determination of the radius at which root transitions take place (Shaing Reference Shaing1984; Turkin et al. Reference Turkin, Beidler, Maaßberg, Murakami, Tribaldos and Wakasa2011).

An example from this density scan is illustrated by the results depicted in figure 2. The profiles of density, temperature and radial electric field are shown for the case with central deuterium/tritium densities of $0.9 \times 10^{20}\ {\rm m}^{-3}$. These profiles for the bulk-ion species are shown in the $n_\sigma$ plot by the blue dotted curve, and their sum is given by the blue solid curve. The density of helium ash is shown by the black dotted curve and the electron density profile is plotted in red. The same colour scheme is used to depict the temperature profiles, although the differences in $T_\sigma$ for the the ion species are too small to be discernable on the plot. The profile of the radial electric field obtained from the differential equation is shown by the black solid curve, and the electron-, unstable- and ion-root solutions of the ambipolarity constraint are given by the black data points. The volume average of the normalised plasma pressure in this simulation has the value of $\langle \beta \rangle = 2.04\,\%$.

Figure 2.

(a) Assumed density profiles of electrons (red solid), ${\rm D}+{\rm T}$ fuel ions (blue solid), D and T individually (blue dotted) and He (black dotted), as functions of normalised radius $\rho = r/a$. (b) Temperature profiles for electrons (red solid) and all ions (blue dotted). (c) Radial electric field calculated through diffusive model (black solid) and the three roots of the ambipolarity equation (dotted). Note the rapid transition from the electron root for $\rho < 0.45$ to the ion root for $\rho > 0.5$.

For this example, the electron root extends nearly to a normalised plasma radius of $\rho =r/a=0.5$. This extent is increased (decreased) if the simulation is performed at lower (higher) density and higher (lower) temperature, as expected theoretically. In this example, the energy confinement time is $\tau _E=1.69$ s, which is somewhat below the $\tau _E=1.99$ s predicted by the International Stellarator Scaling ISS04 (Yamada et al. Reference Yamada2005) for the results of this simulation. The range of densities for which an electron root is possible will increase for smaller $B$, but the optimisation would then need to remain effective at larger values of $\langle \beta \rangle$.

It is interesting to note that, although the electron root is realised for $0 \leqslant \rho \leqslant 0.5$ in this calculation, two other roots also exist as solutions to the ambipolarity equation in this region. These are indicated by dotted curves in figure 2. In the centre of the plasma, the electron root is much stronger (has larger amplitude) but the weaker ion root may also be realised. The diffusive model used to calculate $E_r$ exhibits hysteresis, and the outcome thus depends on the pre-history of the system (Hastings Reference Hastings1986). The electron root is selected if the temperature of the electrons is initially chosen to be higher than that of the ions, as would be the case in a device with dominant electron heating.