Electron root optimisation for stellarator reactor designs

E. Lascas Neto; R. Jorge; C.D. Beidler; J. Lion

doi:10.1017/S0022377824001466

Electron root optimisation for stellarator reactor designs

Published online by Cambridge University Press: 22 January 2025

and

E. Lascas Neto*: Affiliation:
Instituto de Plasmas e Fusão Nuclear, Instituto Superior Técnico, Universidade de Lisboa, 1049-001 Lisboa, Portugal Department of Physics, University of Wisconsin-Madison, Madison, WI 53706, USA
R. Jorge: Affiliation:
Department of Physics, University of Wisconsin-Madison, Madison, WI 53706, USA
C.D. Beidler: Affiliation:
Max-Planck-Institut für Plasmaphysik, D-17491 Greifswald, Germany
J. Lion: Affiliation:
Proxima Fusion GmbH, 81671 Munich, Germany
*: †Email address for correspondence: [email protected]

Article contents

Abstract
Introduction
Theoretical background
Optimisation method
Validation method
Results
Conclusion
Funding
Declaration of interests
References

Rights & Permissions

Abstract

In this work, we propose a method for optimising stellarator devices to favour the presence of an electron root solution of the radial electric field. Such a solution can help avoid heavy impurity accumulation, improve neoclassical thermal ion confinement and helium ash exhaust, and possibly reduce turbulence. This study shows that an optimisation for such a root is possible in quasi-isodynamic stellarators. Examples are shown for both vacuum and finite plasma pressure configurations.

Keywords

fusion plasma plasma confinement

Type: Research Article
Information: Journal of Plasma Physics , Volume 91 , Issue 1 , February 2025 , E24

DOI: https://doi.org/10.1017/S0022377824001466 [Opens in a new window]

NASA ADS Abstract Service [Opens in a new window]
Creative Commons: This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted re-use, distribution and reproduction, provided the original article is properly cited.
Copyright: Copyright © The Author(s), 2025. Published by Cambridge University Press

1 Introduction

Stellarators are a strong candidate for a fusion power plant reactor as they are not prone to current-driven disruptions, have lower recirculating power fraction and have an inherent flexibility of their plasma shape (Helander Reference Helander2014). Such flexibility allows for optimisation of the different plasma and magnetic configurations to achieve the desired reactor properties. This flexibility also comes with the caveat of the complexity of the optimisation problem and the three-dimensional geometry physics. In particular, coil complexity, increased neoclassical transport and turbulence have been seen as obstacles to the design of a reactor-size device of this type. However, recent work has shown that these difficulties can be overcome. Improved coil optimisation techniques have shown that simple coil systems can be achieved (Wechsung et al. Reference Wechsung, Landreman, Giuliani, Cerfon and Stadler2022; Jorge et al. Reference Jorge, Goodman, Landreman, Rodrigues and Wechsung2023; Kappel, Landreman & Malhotra Reference Kappel, Landreman and Malhotra2023). Experimental results in W7-X (Wendelstein 7-X) have shown reduced levels of neoclassical transport (Beidler, Smith & Alonso Reference Beidler, Smith, Alonso, Andreeva, Baldzuhn, Beurskens, Borchardt, Bozhenkov, Brunner and Damm2021). New optimisation techniques have made precise quasisymmetric (QS) (Landreman & Paul Reference Landreman and Paul2022) and quasi-isodynamic (QI) (Goodman et al. Reference Goodman, Camacho Mata, Henneberg, Jorge, Landreman, Plunk, Smith, Mackenbach, Beidler and Helander2023) stellarator configurations numerically achievable. Advances in turbulence optimisation have also been made (Roberg-Clark et al. Reference Roberg-Clark, Plunk, Xanthopoulos, Nührenberg, Henneberg and Smith2023; Jorge et al. Reference Jorge, Dorland, Kim, Landreman, Mandell, Merlo and Qian2024; Kim et al. Reference Kim, Buller, Conlin, Dorland, Dudt, Gaur, Jorge, Kolemen, Landreman and Mandell2024). With this work, we intend to show that stellarators can also be optimised for properties that entail less heavy impurity accumulation, improved helium ash exhaust, good energy confinement with less restrictive optimisation criteria, and possibly, reduced turbulence.

In a fusion reactor, the ion temperature will need to be large in order to achieve sufficient deuterium–tritium reactivity for the fusion reactions to be sustainable. In most reactor-relevant designs, the large heat load reaching the wall is redirected to the divertor. A carbon divertor in a stellarator has been studied in the W7-X experiment (Pedersen et al. Reference Pedersen, König, Krychowiak, Jakubowski, Baldzuhn, Bozhenkov, Fuchert, Langenberg, Niemann and Zhang2018), and partial tungsten coating of the divertor tiles has only been initiated recently in the LHD (Large Helical Device) (Tokitani, Masuzaki & Murase Reference Tokitani, Masuzaki and Murase2019) and in W7-X (Naujoks et al. Reference Naujoks, Dhard, Feng, Gao, Stange, Buttenschön, Bozhenkov, Brezinsek, Brunner and Cseh2023). Tungsten is the first wall material of choice for high-temperature reactor-like designs as it allows for good power exhaust while allowing for lower tritium retention. Experimentally, tungsten divertors have been studied in large tokamak experiments such as ASDEX-U (axially symmetric divertor experiment upgrade) and JET (Joint European Torus) (Pütterich et al. Reference Pütterich, Dux, Neu, Bernert, Beurskens, Bobkov, Brezinsek, Challis, Coenen and Coffey2013). In these experiments, it is seen that the sputtering of the heavy tungsten ions followed by their accumulation in the plasma core, can degrade and even terminate operation due to radiative collapse (Pütterich et al. Reference Pütterich, Dux, Neu, Bernert, Beurskens, Bobkov, Brezinsek, Challis, Coenen and Coffey2013). Heavy impurity transport is thus bound to be an issue in the path to a future stellarator reactor design.

The importance of the radial electric field and the impurity density asymmetries as the main drive mechanisms of a strong heavy impurity core accumulation has been known and considered extensively (Maassberg et al. Reference Maassberg, Brakel, Burhenn, Gasparino, Grigull, Kick, Kuhner, Ringler, Sardei and Stroth1993; Helander et al. Reference Helander, Newton, Mollén and Smith2017; Buller et al. Reference Buller, Smith, Helander, Mollén, Newton and Pusztai2018; Calvo et al. Reference Calvo, Parra, Velasco, Alonso and García-Regaña2018). Optimisation of the heavy impurity density asymmetries through control of the parallel variation of the electrostatic potential has been considered in QI stellarators (Buller, Smith & Mollén Reference Buller, Smith and Mollén2021). However, such an optimisation technique requires a considerable number of Fourier components of the electrostatic potential (Buller et al. Reference Buller, Smith and Mollén2021). Moreover, the available degrees of freedom to change the electrostatic potential are limited in QI designs, as strong flows are likely not present (Helander & Simakov Reference Helander and Simakov2008). Thus, ion cyclotron resonance heating and fast particle asymmetries are the few options to control the heavy impurity density asymmetries (Buller et al. Reference Buller, Smith and Mollén2021).

Another approach to optimise for low impurity confinement is to control the ambipolarity constraint so that the plasma will operate with a core positive radial electric field, often called the electron root. In general, for heavy impurities, due to their high charge, the main drive of neoclassical radial transport is the electric field. Thus, the electron root can help maintain low core heavy impurity accumulation for most situations, as the neoclassical pinch of impurities can point radially outwards (Tamura et al. Reference Tamura, Sudo, Suzuki, Funaba, Nakamura, Tanaka, Yoshinuma and Ida2016), thus helping flush heavy tungsten impurities out of the plasma core. Such an electric field solution may also improve helium ash exhaust, which, if not performed effectively may lead to intolerable fuel dilution (Beidler et al. Reference Beidler, Drevlak, Geiger, Helander, Smith and Turkin2024). Additionally, the transition from a positive to negative electric field often observed when an electron root is present in the plasma core, may generate large enough sheared flows with the potential to suppress turbulence. The effect of an electron root on ion-temperature-gradient turbulence was studied in Fu et al. (Reference Fu, Nicolau, Liu, Wei, Xiao and Lin2021). Electron root plasma operation has been considered experimentally in Pablant et al. (Reference Pablant, Langenberg, Alonso, Beidler, Bitter, Bozhenkov, Burhenn, Beurskens, Delgado-Aparicio and Dinklage2018), Fujiwara et al. (Reference Fujiwara, Kawahata, Ohyabu, Kaneko, Komori, Yamada, Ashikawa, Baylor, Combs and deVries2001) and Yokoyama et al. (Reference Yokoyama, Maaßberg, Beidler, Tribaldos, Ida, Estrada, Castejon, Fujisawa, Minami and Shimozuma2007) and was seen to have beneficial effects on the electron heat confinement (Yokoyama et al. Reference Yokoyama, Maaßberg, Beidler, Tribaldos, Ida, Estrada, Castejon, Fujisawa, Minami and Shimozuma2007). However, such electron root scenarios were mostly obtained on neutral-beam-injection-heated plasmas with low density or on plasmas with high electron temperatures provided by electron-cyclotron resonance heating (ECRH). Optimising for an electron root operation in new optimised stellarator reactor designs without the need for such auxiliary mechanisms is thus an attractive solution to tackle the heavy impurity accumulation problem that is foreseen in future devices.

Finding the electric field solution to the ambipolarity constraint in three-dimensional magnetic configurations is a complex task. This is because, at the values of temperature and density relevant to the reactor plasma core, the ions and electrons can be at different collisionality regimes, leading to a nonlinear equation in the electric field which can have positive (electron root) and negative (ion root) solutions, or both, throughout the radial extent of the device. This equation depends on the radial transport coefficients of the electrons and ions, which in turn depend on the magnetic configuration and plasma parameters. Theoretical work to unravel the relation between the solutions of this equation with the magnetic configuration and plasma parameters is scarce. Efforts have been made to obtain the ion and electron monoenergetic transport coefficients to help identify the space of solutions of the ambipolarity in various stellarator devices (Beidler et al. Reference Beidler, Allmaier, Isaev, Kasilov, Kernbichler, Leitold, Maaßberg, Mikkelsen, Murakami and Schmidt2011) and the study of helium ash transport has been recently done for different W7-X-like configurations in Beidler et al. (Reference Beidler, Drevlak, Geiger, Helander, Smith and Turkin2024). In this work, we propose a new method to optimise for such radial electric field solutions based on maximising the plasma parameter space of electron root solutions. We note that two other studies targeting the optimisation of an electron root (Lee et al. Reference Lee, Lazerson, Smith, Beidler and Pablant2024; Helander et al. Reference Helander, Goodman, Beidler, Kuczynski and Smith2024) were made in parallel to the work presented here. Nevertheless, the three works present substantial differences in the methods used. In Lee et al. (Reference Lee, Lazerson, Smith, Beidler and Pablant2024) the electric field is optimised using the monoenergetic ion transport coefficient obtained with the code DKES (drift kinetic equation solver) and a measure of the effective ripple. The work presented here uses directly the ratio between the ion and electron thermal transport coefficients, and thus a direct calculation of the effective ripple is not necessary. In Helander et al. (Reference Helander, Goodman, Beidler, Kuczynski and Smith2024) an electron root is optimised by tailoring the magnetic field and not by a direct calculation of the neoclassical radial transport coefficients. We start by reviewing the physics behind the cost function and the numerical methods used. Then, we describe the optimisation method used as well as the method used to obtain the electric field solution in the optimised configurations. Finally, we discuss our results by analysing different optimised configurations, comparing the respective electric field solutions and discussing how they affect several properties of interest for a fusion reactor design.

2 Theoretical background

In the standard approach, stellarator optimisation is performed by finding a magnetic configuration that conforms to different aspects of interest defined by a given cost function. Such magnetic configurations are usually bound to obey the magnetohydrodynamics (MHD) equilibrium equation (Freidberg Reference Freidberg2014)

(2.1)\begin{equation} \boldsymbol{J}\times \boldsymbol{B}=\boldsymbol{\nabla} P, \end{equation}

where $\boldsymbol {J}=\boldsymbol {\nabla }\times \boldsymbol {B}/\mu _0$ is the plasma current density, $\boldsymbol {B}$ is the magnetic field and $P$ is the plasma pressure. Equation (2.1) governs the plasma equilibrium of fusion reactors in the absence of strong flows or strong pressure anisotropies. If the equilibrium is assumed to be described by nested closed surfaces of constant toroidal magnetic flux $\psi$, the plasma pressure $P=P(\psi )$ is a flux function. The plasma pressure is the sum of both electron and ion contributions $P=n_eT_e+n_iT_i$, where $n_a(\psi )$ and $T_a(\psi )$ are the density and temperature of species $a$, respectively. In the set of Boozer flux coordinates $(\psi,\theta,\phi )$ (Boozer Reference Boozer1981) the magnetic field $\boldsymbol {B}$ can be written in its contravariant and covariant forms, respectively, as

(2.2)\begin{equation} \boldsymbol{B}=\boldsymbol{\nabla} \psi\times \boldsymbol{\nabla} \alpha =B_{\psi}(\psi,\theta,\phi)\boldsymbol{\nabla} \psi +I(\psi)\boldsymbol{\nabla} \theta +G(\psi) \boldsymbol{\nabla} \phi, \end{equation}

where $\alpha =\theta -\iota (\psi ) \phi$ is the field line label, $\iota$ is the rotational transform, $G$ is $\mu _0/(2{\rm \pi} )$ times the poloidal current outside the surface, $I$ is $\mu _0/(2{\rm \pi} )$ times the toroidal current inside the surface and $B_{\psi }$ is associated with the Pfirsch–Schlütter current.

In fusion reactors, it is of interest to have a magnetic configuration in which the plasma is well confined in order to be able to sustain fusion reactions. In toroidal magnetic confinement devices such as stellarators, particles can drift radially, leading to unwanted transport. This transport is usually dictated by the particles trapped in the magnetic field wells (neoclassical effects) and the anomalous transport due to turbulence. State-of-the-art stellarator configurations are optimised to reduce neoclassical transport across flux surfaces (Beidler et al. Reference Beidler, Smith, Alonso, Andreeva, Baldzuhn, Beurskens, Borchardt, Bozhenkov, Brunner and Damm2021; Landreman & Paul Reference Landreman and Paul2022; Goodman et al. Reference Goodman, Camacho Mata, Henneberg, Jorge, Landreman, Plunk, Smith, Mackenbach, Beidler and Helander2023). Magnetic configurations for which the neoclassical transport across flux surfaces vanishes are called omnigenous. Such magnetic configurations are obtained by tailoring the magnetic field strength properties such that the bounce average of the radial drift of the trapped particles along their bounce orbit is zero. Two particular cases of interest are QS (Landreman & Paul Reference Landreman and Paul2022) and QI (Camacho Mata, Plunk & Jorge Reference Camacho Mata, Plunk and Jorge2022; Jorge et al. Reference Jorge, Plunk, Drevlak, Landreman, Lobsien, Camacho Mata and Helander2022; Goodman et al. Reference Goodman, Camacho Mata, Henneberg, Jorge, Landreman, Plunk, Smith, Mackenbach, Beidler and Helander2023) configurations. In QS devices, the magnetic field strength obeys $B=B(M\theta -N\phi )$, with $M$ and $N$ unique integers. Such a symmetry allows for the drift across flux surfaces to vanish along the bounce orbit of the trapped particles as the magnetic well is symmetric about its minimum value, similar to tokamaks (Helander & Sigmar Reference Helander and Sigmar2002). On the other hand, QI configurations are less restrictive and omnigenity is achieved by imposing

(2.3)\begin{gather} \left. \frac{\partial B_\text{min}(\psi,\alpha)}{\partial \alpha}\right|_{\psi}=\left.\frac{\partial B_\text{max}(\psi,\alpha)}{\partial \alpha}\right|_{\psi} = 0 \ \forall\ \psi, \end{gather}

(2.4)\begin{gather}\left.\frac{\partial l_b(\psi,B,\alpha)}{\partial \alpha}\right|_{\psi,B} = 0 \ \forall\ \psi,B , \end{gather}

in which $B_\text {min}(\psi,\alpha )$ and $B_\text {max}(\psi,\alpha )$ are, respectively, the minimum and the maximum of the magnetic field strength in a particular flux surface and magnetic field line, and $l_b(\psi,B,\alpha )$ is the distance along the particle orbit between two consecutive bounce points. These constraints impose certain properties in the magnetic wells (Cary & Shasharina Reference Cary and Shasharina1997), which are less restrictive than the symmetry imposed by QS. Quasi-isodynamic devices also impose an extra constraint that the contours of $B$ close poloidally to target the lowest possible bootstrap current in an omnigenous configuration (Helander Reference Helander2014). While perfect omnigenity was shown to be impossible to achieve (Cary & Shasharina Reference Cary and Shasharina1997), precise QS and QI magnetic fields can be achieved numerically (Landreman & Paul Reference Landreman and Paul2022; Goodman et al. Reference Goodman, Camacho Mata, Henneberg, Jorge, Landreman, Plunk, Smith, Mackenbach, Beidler and Helander2023) with error fields as small as the geomagnetic field.

Ambipolarity in a fusion plasma is achieved when the net radial currents in the plasma vanish. Such a condition is linked to the temporal evolution of the radial electric field, which is governed by the electric field transport equation. Such an equation is obtained by flux averaging Ampère's law

(2.5)\begin{equation} \varepsilon\frac{\partial E^{\psi}}{\partial t} +J^{\psi}=0, \end{equation}

in which $E^{\psi }=\langle \boldsymbol {E}\boldsymbol {\cdot }\boldsymbol {\nabla } \psi \rangle$ and $J^{\psi }=\langle \boldsymbol {J}\boldsymbol {\cdot }\boldsymbol {\nabla } \psi \rangle$ are the flux-averaged contravariant components of the electric field $\boldsymbol {E}$ and the plasma current density across flux surfaces, and $\varepsilon$ is the plasma dielectric constant. In (2.5), the fact that the plasma obeys (2.1) is used. Thus, the currents generated by the magnetic field have no component across the flux surface, i.e. $(\boldsymbol {\nabla } \times \boldsymbol {B})\boldsymbol {\cdot }\boldsymbol {\nabla }\psi =0$. Therefore, in (2.5), $J^{\psi }$ is only generated by the plasma (neglecting external sources) and can be obtained from the distribution function $f^a$ of species $a$. Expanding the distribution function in $\varDelta _{\psi }/L$, with $\varDelta _{\psi }$ the orbit width and $L$ the typical macroscopic length of the plasma gradients, we can write $f^a=f_M^a+f_1^a+f_h^a$, with $f_M^a$ a zeroth-order Maxwellian, $f^a_1$ the first-order part of $f_a$ and $f_h^a$ the part with higher-order terms (Hastings Reference Hastings1984; Shaing Reference Shaing1984). The flux-averaged radial current density due to the plasma can then be written as

(2.6)\begin{align} J^{\psi} & = \sum_a \left \langle\int Z_a e \boldsymbol{v_d}\boldsymbol{\cdot}\boldsymbol{\nabla} \psi f^a \,{\rm d}^3 v \right \rangle \nonumber\\ & =\sum_a \left\langle\int Z_a e \boldsymbol{v_d}\boldsymbol{\cdot}\boldsymbol{\nabla} \psi f^a_1 \,{\rm d}^3 v \right\rangle+\sum_a \left\langle\int Z_a e \boldsymbol{v_d}\boldsymbol{\cdot}\boldsymbol{\nabla} \psi f^a_h \,{\rm d}^3 v \right\rangle \nonumber\\ & = \sum_a e Z_a \varGamma_a + \varepsilon \left(\frac{{\rm d}V}{{\rm d}\psi}\right)^{{-}1}\frac{\partial}{\partial \psi}\left[\frac{{\rm d}V}{{\rm d}\psi}\sum_a (e Z_aD_{E^{\psi}}^a) \frac{\partial E_{\psi}}{\partial\psi}\right]. \end{align}

In (2.6), $\varGamma _a=\langle \boldsymbol {\varGamma } \boldsymbol {\cdot }\boldsymbol {\nabla } \psi \rangle$ is the usual flux-averaged particle-flux-density across flux surfaces obtained at lowest order in $\varDelta _{\psi }/L$, and $D_{E_{\psi }}^a$ is the electric field diffusivity associated with species $a$ which comes from considering higher-order orbit width effects (Hastings Reference Hastings1984; Shaing Reference Shaing1984), $e$ is the electron charge, $Z_a$ is the charge number of species $a$, $\boldsymbol {v_d}$ is the drift velocity of a particle and $V(\psi )$ is the volume enclosed by the flux surface $\psi$. If the higher-order orbit width effects from $f_h^a$ are neglected, the steady-state radial electric field is obtained from the ambipolarity constraint

(2.7)\begin{equation} \sum_a e Z_a \varGamma_a =0, \end{equation}

which can have multiple solutions due to the different collisionality regimes of the electrons and ions. The choice between these multiple solutions depends on thermodynamical considerations (Shaing Reference Shaing1984; Hastings, Houlberg & Shaing Reference Hastings, Houlberg and Shaing1985) and large transitions of the electric field sign between neighbouring flux surfaces can occur. In these regions, the gradient of the electric field gets so large that higher-order orbit width effects from $f_h^a$ become important. When considering this correction, one obtains the electric field diffusion term in (2.6) and information about the electric field gradient is retained (Hastings Reference Hastings1984). This extra term is related to the minimisation of the heat transfer rate (Shaing Reference Shaing1984). Taking into account this extra term, the electric field transport equation can be written as

(2.8)\begin{equation} \varepsilon\frac{\partial E^{\psi}}{\partial t} - \varepsilon\left(\frac{{\rm d}V}{{\rm d}\psi}\right)^{{-}1}\frac{\partial}{\partial \psi}\left[\frac{{\rm d}V}{{\rm d}\psi}\sum_a \left(e Z_aD_{E^{\psi}}^a\right) \frac{\partial E^{\psi}}{\partial \psi}\right] = \sum_a e Z_a \varGamma_a, \end{equation}

yielding a single solution for the radial electric field. However, the electric field affects the power balance in a fusion reactor and vice versa. Therefore, to obtain a consistent solution of the electric field, (2.8) has to be solved together with a set of energy transport equations. A common form of the energy transport equation for species $a$ is (Turkin et al. Reference Turkin, Beidler, Maaberg, Murakami, Tribaldos and Wakasa2011)

(2.9)\begin{equation} \frac{3}{2}\frac{\partial n_{a}T_{a}}{\partial t}+\left(\frac{{\rm d}V}{{\rm d}\psi}\right)^{{-}1}\frac{\partial}{\partial \psi}\left[\frac{{\rm d}V}{{\rm d}\psi}\left(Q_a+\varGamma_aT_a-\chi_a n_a\frac{\partial T_a}{\partial \psi}\right)\right]=\mathcal{P}_a+e Z_a \varGamma_a E_{\psi}, \end{equation}

where $Q_a=\langle Q \boldsymbol {\cdot }\boldsymbol {\nabla } \psi \rangle$ is the flux-averaged energy-flux-density of species $a$ across flux surfaces, $\chi _a$ is the anomalous heat diffusivity coefficient of species $a$ and $E_{\psi }=\boldsymbol {E}\boldsymbol {\cdot }\boldsymbol {e}_{\psi }$. On the right-hand side of (2.9), the first term, $\mathcal {P}_a$, is the total power balance of external energy sources and sinks related to species $a$ and the second term is the work done by the electric field in the fluid of species $a$. If the density of the different species is allowed to vary, then a set of continuity equations also needs to be solved. A common form of the continuity equation of species $a$ is (Turkin et al. Reference Turkin, Beidler, Maaberg, Murakami, Tribaldos and Wakasa2011)

(2.10)\begin{equation} \frac{\partial n_{a}}{\partial t} +\left(\frac{{\rm d}V}{{\rm d}\psi}\right)^{{-}1}\frac{\partial}{\partial \psi}\left[\frac{{\rm d}V}{{\rm d}\psi}\left(\varGamma_a-\mathcal{D}_a \frac{\partial n_a}{\partial \psi}\right)\right]=S_a, \end{equation}

in which $\mathcal {D}_a$ is the particle diffusivity of species $a$ associated with the anomalous transport and $S_a$ is the source term of particles of species $a$. To close the system of equations (2.8)–(2.10) we use the quasineutrality constraint $\sum _a Z_a n_a=0$. The particle and energy-flux-densities can be prescribed using a reduced model or obtained from solving the drift kinetic equation, while the anomalous diffusivities are obtained either from a reduced model or from solving the gyrokinetic equation. We note that in the previous discussion, we have assumed that the fluxes contributing to the ambipolarity term of the radial current density in (2.6) are the neoclassical fluxes. Such an assumption comes from the fact that in stellarators, at the usual gyrokinetic orderings, the anomalous transport is intrinsically ambipolar (Sugama et al. Reference Sugama, Okamoto, Horton and Wakatani1996; Helander Reference Helander2014) and will thus have no influence on solutions of (2.7). Thus, unless a perfect omnigenous magnetic configuration is achieved to reduce neoclassical fluxes to extremely low levels, the neoclassical fluxes are the most important players in the ambipolarity constraint. Therefore, we now show how the neoclassical theory can help obtain the distribution function $f_a$, and find $\varGamma _a$ and $Q_a$.

The usual local neoclassical theory allows us to obtain the steady-state solution of the first-order distribution function $f_1^a$ which is dictated by the drift kinetic equation (Landreman et al. Reference Landreman, Smith, Mollén and Helander2014)

(2.11)\begin{align} & \dot{\psi}\frac{\partial f_1^a}{\partial \psi}+\dot{\theta}\frac{\partial f_1^a}{\partial \theta}+\dot{\phi}\frac{\partial f_1^a}{\partial \phi}+\dot{x}\frac{\partial f_1^a}{\partial x}+\dot{\xi}\frac{\partial f_1^a}{\partial \xi}-C(f^a _1) \nonumber\\ & \quad=\frac{T_a}{ZeB^3\sqrt{g}}x^2(1+\xi^2)\left(G\frac{\partial B}{\partial \theta}-I\frac{\partial B}{\partial \phi}\right)f_M^a (A_1^a+x^2A_2^a)-\frac{v_{{\parallel}}}{B_0} A_3^a f_M^a, \end{align}

where the velocity variables are the normalised velocity $x=v/v_{th}^a$ and the pitch angle $\xi =v_{\parallel }/v$, with $v_{th}^a=\sqrt {2T_a/m_a}$ the thermal velocity of species $a$ and $m_a$ the mass of species $a$. We also used the fact that the magnetic drift in Boozer coordinates can be written as $\boldsymbol {v_m}\boldsymbol {\cdot }\boldsymbol {\nabla }\psi =-({T_a}/{ZeB^3\sqrt {g}})x^2(1+\xi ^2)(G({\partial B}/{\partial \theta })-I({\partial B}/{\partial \phi }))$. The Jacobian determinant $\sqrt {g}$ in (2.11) is written as $\sqrt {g}=(G+\iota I)/B^2$ in Boozer coordinates. In (2.11), the thermodynamical forces are given by

(2.12)\begin{gather} A_1^a=\frac{1}{n_a}\frac{{\rm d} n_a}{{\rm d} \psi}-\frac{Z_a e E_{\psi}}{T_a}-\frac{3}{2}\frac{1}{T_a}\frac{{\rm d} T_a}{{\rm d}\psi}, \end{gather}

(2.13)\begin{gather}A_2^a=\frac{1}{T_a}\frac{{\rm d} T_a}{{\rm d} \psi}, \end{gather}

(2.14)\begin{gather}A_3^a={-}\frac{Z_a e B_0}{T_a} \frac{\langle E_{{\parallel}}B\rangle}{\langle B^2\rangle}, \end{gather}

with the flux surface average of a quantity $X$ defined as $\langle X \rangle = \int _0^{2{\rm \pi} } \int _0^{2{\rm \pi} } {\rm d}\theta \,{\rm d}\phi \,X \sqrt {g} / ({\rm d}V/ {\rm d}\psi )$. In the local approximation, $\dot {\psi }$ is usually neglected, and at a large aspect ratio we can write the poloidal velocity of the particles as

(2.15)\begin{equation} \dot{\theta}=v_{{\parallel}}\frac{\boldsymbol{B}}{B} \boldsymbol{\cdot}\boldsymbol{\nabla} \theta + \boldsymbol{v}_{E}\boldsymbol{\cdot}\boldsymbol{\nabla}\theta=\frac{xv_{th}^a\xi}{B\sqrt{g}}\iota -\frac{G}{B^2 \sqrt{g}}E_{\psi}=\left(x\xi-E^*_{\psi}\right)\frac{\iota v_{th}^a}{B\sqrt{g}}, \end{equation}

and the toroidal velocity as

(2.16)\begin{equation} \dot{\phi}=v_{{\parallel}}\frac{\boldsymbol{B}}{B}\boldsymbol{\cdot}\boldsymbol{\nabla} \phi+ \boldsymbol{v}_{E}\boldsymbol{\cdot}\boldsymbol{\nabla}\phi=\frac{xv_{th}^a\xi}{B\sqrt{g}} +\frac{1}{B^2 \sqrt{g}}E_{\psi}=\left(x\xi+E^*_{\psi}\frac{I \iota}{G}\right)\frac{ v_{th}^a}{B\sqrt{g}}, \end{equation}

where the magnetic drift tangent to the flux surfaces has been neglected, and thus the poloidal and toroidal velocity of the particles combine the motion along field lines and the $E\times B$ drift defined as $\boldsymbol {v_{E}}$. In (2.15) and (2.16), $E^*_{\psi }=E_{\psi }G/(\iota v_{th}^a B)$ is a measure of the strength of the $E\times B$ drift (drift frequency) relative to the parallel orbit motion of the particle (transit frequency). The parallel acceleration of the particles is given by

(2.17)\begin{align} \dot{\xi}& ={-}\frac{(1-\xi^2)}{\xi 2 B^2}\boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla} B +\frac{Z_a e}{2 T_a x^2 \xi}(1-\xi^2)E_{\psi} \boldsymbol{v_m} \boldsymbol{\cdot}\boldsymbol{\nabla}\psi \nonumber\\ & ={-}\frac{x(1-\xi^2)v_{th}^a}{2B^2\sqrt{g}}\left(\iota\frac{\partial B}{\partial \theta}+\frac{\partial B}{\partial \phi}\right)-\frac{\iota v_{th}^a}{2 \sqrt{g}B^2}\xi(1-\xi^2)\left(\frac{\partial B}{\partial \theta}-\frac{I}{G}\frac{\partial B}{\partial \phi}\right)E^*_{\psi}, \end{align}

and the kinetic energy temporal variation by

(2.18)\begin{equation} \dot{x}=\frac{Z_a e}{2 T_a x}E_{\psi} \boldsymbol{v_m} \boldsymbol{\cdot}\boldsymbol{\nabla} \psi={-}\frac{\iota v_{th}^a}{2 \sqrt{g}B^2}x(1+\xi^2)\left(\frac{\partial B}{\partial \theta}-\frac{I}{G}\frac{\partial B}{\partial \phi}\right)E^*_{\psi}. \end{equation}

These two equations show that trapped/passing domains are influenced by the parallel mirror force as well as by the effect of a strong radial electric field. The set of (2.15)–(2.18) for the orbits of the particles in the phase space are usually called ‘full trajectories’ (Landreman et al. Reference Landreman, Smith, Mollén and Helander2014). If we neglect the electric field corrections in (2.17) and (2.18), we obtain the so-called ‘partial trajectories’ (Landreman et al. Reference Landreman, Smith, Mollén and Helander2014).

It is useful to linearise the drift kinetic equation, (2.11), to obtain an approximate, but more expeditiously and less numerically expensive result. Such a task is accomplished by looking at the right-hand side term of (2.11), which is associated with the thermodynamical forces in (2.12)–(2.14). In (2.11), it can be observed that each thermodynamical force term is associated with a particular dependency on the velocity space variables $(x,\xi )$. The first-order correction to the distribution function can thus be linearised as

(2.19)\begin{equation} f_1^a=\frac{n_a}{{\rm \pi}^{{3}/{2}}{v^a_{th}}^{3}}\sum_{p=1}^3 f_{1p}^a A_p^a \end{equation}

which leads to the set of equations

(2.20)\begin{align} & (\xi x-E^*_{\psi})\iota \frac{\partial f_{1p}^a}{\partial \theta}+\left(\xi x+\frac{\iota I}{G}E^*_{\psi}\right)\frac{\partial f_{1p}^a}{\partial \phi}-\frac{x(1-\xi^2)}{2}\left(\frac{\iota}{B}\frac{\partial B}{\partial \theta}+\frac{1}{B}\frac{\partial B}{\partial \phi}\right)\frac{\partial f_{1p}^a}{\partial \xi}\nonumber\\ & \quad -\frac{(1-\xi^2)}{2}\left(\frac{\iota}{B}\frac{\partial B}{\partial \theta}-\frac{\iota I}{G B}\frac{\partial B}{\partial \phi}\right)\left[\xi\frac{f_{1p}^a}{\partial \xi}+x\frac{f_{1p}^a}{\partial x}\right]E_{\psi}(\psi)=\frac{B_0}{B}\nu'C(f_{1p}^a)+\mathcal{S}^a_p, \end{align}

where

(2.21)\begin{gather} \mathcal{S}_1^a=\frac{T_a G}{Z_aeB_0}x^2(1+\xi^2)e^{{-}x^2}\frac{B_0}{B}\left(\frac{1}{B}\frac{\partial B}{\partial \theta}-\frac{ I}{G B}\frac{\partial B}{\partial \phi}\right), \end{gather}

(2.22)\begin{gather}\mathcal{S}_2^a=\frac{T_a G}{Z_aeB_0}x^4(1+\xi^2)e^{{-}x^2}\frac{B_0}{B}\left(\frac{1}{B}\frac{\partial B}{\partial \theta}-\frac{ I}{G B}\frac{\partial B}{\partial \phi}\right), \end{gather}

(2.23)\begin{gather}\mathcal{S}_3^a={-}x\xi e^{{-}x^2}\frac{(G+\iota I)}{B_0}. \end{gather}

In (2.20), $\nu ^{\prime }=\nu _{aa}(G+\iota I)/v_{th}^a$ and $\nu _{aa}=n_a(Z_a e)^4\ln \varLambda /(4{\rm \pi} ^4\epsilon _0^2m_a^2{v_{th}^a}^3)$ is the typical self-collisions frequency of species $a$ (Helander & Sigmar Reference Helander and Sigmar2002), $\ln \varLambda$ is the Coulomb logarithm and $B_0$ is the magnetic field strength on-axis. The set of (2.20)–(2.23) describes the linear approximation of the 4-D (four-dimensional) drift kinetic equation, which, throughout this work, we will refer to as the ‘4-D approach’ method to solve the drift kinetic equation. The radial particle-flux-density of species $a$ can then be calculated as

(2.24)\begin{equation} \varGamma_a={-}n_a\left[L_{11}^a A_1^a+L_{12}^a A_2^a+L_{13}^a A_3^a\right], \end{equation}

the radial energy-flux-density as

(2.25)\begin{equation} Q_a={-}n_a T_a\left[L_{21}^a A_1^a+L_{22}^a A_2^a+L_{23}^a A_3^a\right] \end{equation}

and the parallel flow as

(2.26)\begin{equation} \frac{\langle U_{a_{{\parallel}}}B\rangle}{B_0}={-}\left[L_{31}^a A_1^a+L_{32}^a A_2^a+L_{33}^a A_3^a\right]. \end{equation}

The thermal transport coefficients are thus defined as

(2.27)\begin{gather} L_{1p}^a={-}\left\langle\int {{\rm d}x}^3 {\rm \pi}^{-{3}/{2}}f_{1p}^a \boldsymbol{v_d}\boldsymbol{\cdot}\boldsymbol{\nabla} \psi \right\rangle, \end{gather}

(2.28)\begin{gather}L_{2p}^a={-}\left\langle\int {{\rm d}x}^3 {\rm \pi}^{-{3}/{2}} f_{1p}^a x^2 \boldsymbol{v_d}\boldsymbol{\cdot}\boldsymbol{\nabla} \psi \right\rangle, \end{gather}

(2.29)\begin{gather}L_{3p}^a={-}\left\langle\int {{\rm d}x}^3 {\rm \pi}^{-{3}/{2}} \frac{B}{B_0} v_{th}^a f_{1p}^a \xi x \right\rangle. \end{gather}

These coefficients contain the pieces of the distribution function related to the neoclassical radial particle-flux-density, radial energy-flux-density and parallel flow, respectively. We note that these coefficients depend on the flux coordinate $\psi$, the magnetic configuration, the inverse mean free path $\nu _{aa}/v_{th}^a$ and the normalised $E\times B$ velocity, $E_{\psi }/(v_{th}^aB)$, or equivalent parameters like $\nu '$ or collisionality $\nu ^*=\nu _{aa} R_0/(\iota v_{th}^a)$, and $E^*_{\psi }$, where $R_0$ is the plasma major radius. The monoenergetic approximation further simplifies (2.20)–(2.23) by assuming the particle velocity to be a parameter. This assumption implies the energy $x$ to be a parameter rather than a variable. With this assumption, the set of equations (2.20)–(2.23) only needs to be solved for the $p=1$ and $p=3$ components, because the components $p=1$ and $p=2$ carry the same information as they have identical dependencies in $\xi$. Additionally, because this monoenergetic description assumes $\dot {x}=0$, the electric field terms in (2.18) and (2.17) have to be dropped for consistency, and thus only the ‘partial trajectories’ can be considered in such a description. The monoenergetic transport coefficients can then be defined as

(2.30)\begin{equation} \left.\begin{gathered} D_{1p}={-}\frac{1}{2}\left\langle\int^{1}_{{-}1} {\rm d}\xi \,f_{1p}^a e^{x^2} \boldsymbol{v_d}\boldsymbol{\cdot}\boldsymbol{\nabla} \psi \right\rangle, \\ D_{3p}={-}\frac{1}{2}\left\langle \int^{1}_{{-}1} {\rm d}\xi \,v_{th}^a\frac{B}{B_0} f_{1p}^a e^{x^2} \xi x \right\rangle,\quad p=1,3, \end{gathered}\right\}\end{equation}

with $D_{22}=D_{12}=D_{21}=D_{11}$, $D_{23}=D_{13}$ and $D_{32}=D_{31}$. The thermal coefficients can then be evaluated by energy convolution of the monoenergetic coefficients

(2.31)\begin{equation} L_{ij}^a=\frac{2}{\sqrt{\rm \pi}}\int_0^{\infty} {{\rm d}x}^2 \,D_{ij}(x) xe^{{-}x^2}h_ih_j, \end{equation}

with $h_1=h_3=1$ and $h_2=x^2$.

The DKES code solves the drift kinetic equation by solving the monoenergetic linearised set of equations using a variational principle (van Rij & Hirshman Reference van Rij and Hirshman1989). Therefore, the set of the two monoenergetic equations in such a code has to be cast in a conservative form. To achieve this, the DKES code makes use of the ‘partial trajectories’ with the extra assumption that the $E\times B$ drift is incompressible such that

(2.32)\begin{equation} v_{E}\approx\frac{\boldsymbol{B}\times \boldsymbol{E}}{\langle B^2 \rangle}. \end{equation}

With this extra approximation, one has the so-called ‘DKES trajectories’. The SFINCS code (Landreman et al. Reference Landreman, Smith, Mollén and Helander2014) uses a different numerical scheme and is able to solve (2.11), or both the 4-D and monoenergetic approaches to the set of (2.20)–(2.23), using the various options which SFINCS provides for describing the trajectories.

A moment should be taken here to note the difference between the thermal and monoenergetic transport coefficients, which will become important when understanding choices made in the optimisation method. Notice that the thermal coefficients are written with a superscript associated with the species $a$ while the monoenergetic ones are not. This is due to the thermal velocity dependence of the $\nu _{aa}/v_{th}^a$ and $E_{\psi }/(B v_{th}^a)$ variables. We thus need to specify the species temperature and density to obtain the thermal transport coefficients. In the monoenergetic case, the variables of interest become $\nu _{ii}(v=v_{th}^i)/v$ (in which the ion–ion collision frequency of ions is chosen as the reference value) and $E_{\psi }/(B v)$, where the $v$ is a parameter and appears only for normalisation purposes. Therefore, the monoenergetic coefficients are species-independent. This makes the 4-D approach quite useful when we seek to compute values of the transport coefficients at specific collisionalities and electric field values for different species easily. Yet, the monoenergetic approach is more suitable in cases where we intend to create a full species-independent scan of the transport coefficients as a basis to analyse a specific magnetic configuration.

The different terms in the drift kinetic equation will lead to different collisionality regimes in which the transport coefficients scale differently with the collisionality. While these regimes have been extensively analysed (Shaing, Ida & Sabbagh Reference Shaing, Ida and Sabbagh2015), we point here to some physical aspects of the main regimes. At high collisionality, the most important mechanism for radial transport is collisions and the radial transport scales linearly with collisionality. When collisionality is close to unity, a plateau regime may form due to the resonance between the orbit movement along the magnetic field and collisions. In this case, radial transport is independent of collisionality and proportional to the transit frequency $\omega _T=\iota v_{th}^a/R_0$. At lower collisionalities, several regimes may exist. If the $E\times B$ drift is less important than collisions, the latter will limit the radial drift of particles. The stronger the collisions the more the radial excursions of trapped particles are reduced and thus the transport coefficients scale with the inverse of collisionality. This is the $1/\nu ^*$ regime that is problematic in stellarators as the heat fluxes increase strongly with temperature in this regime (Beidler et al. Reference Beidler, Allmaier, Isaev, Kasilov, Kernbichler, Leitold, Maaßberg, Mikkelsen, Murakami and Schmidt2011, Reference Beidler, Drevlak, Geiger, Helander, Smith and Turkin2024). When the electric field magnitude is strong enough for the $E\times B$ poloidal drift to be more important than collisions, the radial excursion of the particles will be limited by this poloidal drift. Two asymptotic scalings can be identified in such a case, a scaling with $\sqrt {\nu ^*}$ and a linear scaling with $\nu ^*$ depending on which detrapping mechanism is more relevant, collisions or drift-associated collisionless trapping detrapping processes. In both of these regimes, the radial excursions of trapped particles are also reduced with increasing $E\times B$ drift frequency. These two asymptotic regimes may not be well separated and in such cases, other intermediate scalings may coexist. The fact that different species may be at different regimes in low collisionalities is the reason why the ambipolarity constraint can have multiple solutions.

Finally, we discuss the collision operator $C(f_1)$. There are three widely used forms of the linearised collision operator considered in drift kinetic calculations. The first is the linearised Fokker–Planck collision operator (Landreman et al. Reference Landreman, Smith, Mollén and Helander2014), which considers both energy and momentum collisional transfer between the first-order distribution and the zeroth-order Maxwellian distribution, as well as the respective correction terms that ensure energy and momentum conservation. If the transfer of energy is neglected, one has the second form which considers only a Lorentz scattering operator and a momentum correction term. The third form drops the momentum correction term and considers only a Lorentz scattering operator. In this work, we will consider the third and simplest form of the collision operator unless stated otherwise. The Lorentz scattering operator can be written as

(2.33)\begin{equation} C(f_1)=\frac{\nu_D(x)}{2\nu_{aa}}\frac{\partial}{\partial \xi}\left[(1-\xi^2)\frac{\partial}{\partial \xi}\right] \end{equation}

where

(2.34)\begin{equation} \nu_D(x)=\nu_{aa}\left(\text{erf}(x)-\frac{\text{erf}(x)}{2x^2}+\frac{1}{2x}({\rm d} [\text{erf}(x)]/{{\rm d}x})\right). \end{equation}

We note that such a choice is expected to be enough to treat the main ion and electron species at the low collisionality regime important for the optimisation step as will be discussed in the next section.

We can now determine the ambipolarity constraint for a simple deuterium (D) and tritium (T) plasma. In this case, quasineutrality implies that $n_e=n_i=n_D+n_T$. In general, the tritium concentration will be a fraction of deuterium, $n_T=\gamma n_D$, with a usual choice of equal concentrations $\gamma =1$. Assuming ion species of equal temperature such that $T_D=T_T=T_i$ we thus can write the ambipolarity constraint as $\varGamma _e=\varGamma _D+\varGamma _T$ which, as shown in Beidler et al. (Reference Beidler, Smith, Alonso, Andreeva, Baldzuhn, Beurskens, Borchardt, Bozhenkov, Brunner and Damm2021, Reference Beidler, Drevlak, Geiger, Helander, Smith and Turkin2024), has the solution

(2.35)\begin{equation} \frac{e E_{\psi}}{T_i}=\frac{1}{\mathcal{A}}\left[\left(1-\frac{L_{11}^e}{L_{11}^i}\right)\frac{1}{n_i}\frac{{\rm d}n_i}{{\rm d} \psi} +\delta_{12}^i\frac{1}{T_i}\frac{{\rm d}T_i}{{\rm d} \psi}-\frac{L_{11}^e}{L_{11}^i}\delta_{12}^e \frac{1}{T_e}\frac{{\rm d}T_e}{{\rm d}\psi}\right], \end{equation}

with

(2.36)\begin{gather} \mathcal{A}=1+\frac{L_{11}^e}{L_{11}^i}\frac{T_i}{T_e}, \end{gather}

(2.37)\begin{gather}L_{11}^i=\frac{L_{11}^D+\gamma L_{11}^T}{1+\gamma}, \end{gather}

(2.38)\begin{gather}\delta_{12}^a=\frac{L_{12}^a}{L_{11}^a} \end{gather}

and

(2.39)\begin{equation} \delta_{12}^i=\frac{L_{11}^D\delta_{12}^D +\gamma L_{11}^T\delta_{12}^T }{L_{11}^D +\gamma L_{11}^T}. \end{equation}

As stated in Beidler et al. (Reference Beidler, Drevlak, Geiger, Helander, Smith and Turkin2024), for relevant fusion plasmas, it is usually verified that $\delta _{12}^a > 0$, $\textrm {d}n_a/\textrm {d}\psi n_a^{-1}<0$ and $\textrm {d}T_a/\textrm {d}\psi T_a^{-1}<0$. Thus, a condition that maximises the parameter space of positive electric field solutions (electron root) is to maximise the ratio $L_{11}^e/L_{11}^i$. Indeed, experimental data was seen to link large $L^e_{11}$ values to the occurrence of a positive electric field in LHD plasmas (Velasco et al. Reference Velasco, Calvo, Satake, Alonso, Nunami, Yokoyama, Sato, Estrada, Fontdecaba and Liniers2016). In the following, we consider the ambipolar electric field to be set by (2.35) and, therefore seek stellarator solutions that maximise this ratio. We note that if $L^e_{11}/L^i_{11}$ is not sufficiently large, then although an electron root could be found, it could only be done at a restricted space of temperatures and density profiles. Therefore, maximising the ratio $L^e_{11}/L^i_{11}$ as much as possible while achieving low heat transport is important to allow for a large space of plasma parameters to be explored for other optimisation metrics, while still achieving plasmas which enable electron root operation.

3 Optimisation method

Our main goal is to obtain a magnetic field equilibrium that allows for an electron root solution at a large space of plasma parameters. For this purpose, we use the MHD VMEC (variational moments equilibrium code) equilibrium code (Hirshman & Whitson Reference Hirshman and Whitson1983), together with the optimisation framework SIMSOPT (Landreman et al. Reference Landreman, Medasani, Wechsung, Giuliani, Jorge and Zhu2021). The VMEC code uses the set of flux coordinates $(s,\theta _V,\phi _V)$ where $s=\psi /\psi _b$ is the flux label with $\psi _b$ the toroidal flux at the plasma boundary, $\phi _V$ is the toroidal angle in cylindrical coordinates and $\theta _V$ is a poloidal angle. The VMEC discretises the flux surface label $s$ with a uniform radial grid and treats the $(\theta _V,\phi _V)$ in the Fourier space. A flux surface $S$ is then parameterised using the usual cylindrical coordinate system as $S=[R(\theta _V,\phi _V)\cos (\phi _V),R(\theta _V,\phi _V)\sin (\phi _V),Z(\theta _V,\phi _V)]$ with

(3.1)\begin{gather} R(\theta_V,\phi_V)=\sum_{m=0}^{M_\text{pol}}\sum_{n={-}N_{\text{tor}}}^{N_{\text{tor}}}\text{RBC}_{m,n}\cos(m\theta_V-n\phi_V), \end{gather}

(3.2)\begin{gather}Z(\theta_V,\phi_V)=\sum_{m=0}^{M_\text{pol}}\sum_{n={-}N_{\text{tor}}}^{N_{\text{tor}}}\text{ZBS}_{m,n}\sin(m\theta_V-n\phi_V). \end{gather}

The integers $M_\textrm {pol}$ and $N_\textrm {tor}$ refer to the number of poloidal and toroidal mode numbers used to describe the Fourier decomposition, respectively. The decomposition in (3.1) and (3.2) set the sine terms of $R$ and the cosine terms of $Z$ to zero to ensure stellarator-symmetry (Dewar & Hudson Reference Dewar and Hudson1998). Such a symmetry simplifies the problem by reducing the degrees of freedom, while still achieving good optimised solutions. In this work, we use the fixed boundary mode of VMEC where, given the boundary of the plasma, the plasma pressure $P(\psi )$ and net toroidal current $I(\psi )$ profiles, VMEC will converge to a solution in which at each flux label $s$ the MHD force balance equation is satisfied. In stellarator optimisation, we are interested in finding solutions with net zero toroidal current and so we set $I(\psi )=0$. Optimisations are often done for the vacuum solution of the MHD force balance in which $P=0$. Solutions which include the correction effects of the plasma pressure are called finite-$\beta$ solutions with $\beta =P/(2\mu _0 B^2)$ the ratio between the plasma and magnetic energies. The magnetic configuration in VMEC coordinates is transformed to Boozer coordinates using the BOOZ_XFORM code (Sanchez et al. Reference Sanchez, Hirshman, Ware, Berry and Spong2000).

To target a precise QI magnetic geometry, we use the cost function

(3.3)\begin{equation} J_\text{QI}=\frac{n_\text{fp}}{4{\rm \pi}^2 (B_\text{max}-B_\text{min})^2}\int^{2{\rm \pi}}_0 {\rm d}\alpha \int^{2{\rm \pi} n_\text{fp}}_0 {\rm d}\phi\left[B( s,\alpha,\phi)-B_\text{QI}( s,\alpha,\phi)\right]^2 \end{equation}

proposed in Goodman et al. (Reference Goodman, Camacho Mata, Henneberg, Jorge, Landreman, Plunk, Smith, Mackenbach, Beidler and Helander2023). In this function, $n_\textrm {fp}$ is the number of field periods and $B_\textrm {QI}$ is the target magnetic field constructed from $B(s,\alpha,\phi )$ in such a way that its magnetic wells satisfy (2.3) and (2.4). We follow Goodman et al. (Reference Goodman, Camacho Mata, Henneberg, Jorge, Landreman, Plunk, Smith, Mackenbach, Beidler and Helander2023) and avoid large stresses in the coils by constraining the overall mirror of the magnetic field

(3.4)\begin{equation} \varDelta_M=\frac{B_\text{max}-B_\text{min}}{B_\text{max}+B_\text{min}}, \end{equation}

to be below a certain threshold. To avoid small minor radius solutions to be favoured during optimisation, we constrain the aspect ratio $A$ of the configuration to be close to the aspect ratio of the initial configuration. We add the additional constraint to the mean of the rotational transform $\bar {\iota }$ in order to avoid low-order rational surfaces. Thus, the total cost function that directly targets the magnetic geometry properties can be written as

(3.5)\begin{equation} J_1 = w_\text{QI}J_\text{QI}+w_{\iota}(\bar{\iota}-\bar{\iota}_\text{target})^2 +w_{A}(A-A_\text{target})^2+w_{M}\max [0,(\varDelta_M-\varDelta_{M_\text{target}})]^2, \end{equation}

where $w_\textrm {QI}$, $w_{\iota }$, $w_{A}$ and $w_{M}$ are the weights that multiply the different cost functions and can be adjusted in order to give more focus to one or more of the targets during optimisation.

As seen in the previous section, maximising the plasma parameter space of the electron root solution can be achieved by selecting a magnetic configuration that maximises the ratio of the radial diffusion coefficients $L_{11}^e/L_{11}^i$. Thus, the choice here for a cost function to target magnetic configurations prone to have an electron root is

(3.6)\begin{equation} J_{E_{\psi}} = \left(\frac{L_{11}^i}{L_{11}^e} \right)^2, \end{equation}

where we emphasise again the dependence of these coefficients in the magnetic configuration and the parameters $\nu '$ and $E_{\psi }^*$. There are several codes which are suitable to calculate the transport coefficients $L_{11}^i$ and $L_{11}^e$. Two of these codes have been used consistently for stellarator transport calculations: DKES (van Rij & Hirshman Reference van Rij and Hirshman1989) and SFINCS (Landreman et al. Reference Landreman, Smith, Mollén and Helander2014). We choose to use SFINCS for the calculation of our electron root cost function $J_{E_{\psi }}$. This is because DKES is only able to calculate the transport matrix in the monoenergetic approximation, while SFINCS is capable of solving the linearised drift kinetic equation in both monoenergetic and 4-D approaches (Landreman et al. Reference Landreman, Smith, Mollén and Helander2014). We preferably want to use the 4-D approach to calculate the cost function $J_{E_{\psi }}$. To understand this choice, it is instructive to look again into the difference between the monoenergetic and thermal transport coefficients pointed out in the last section. Equation (3.6) involves calculating at each optimisation loop the thermal radial transport coefficients at least for one ion species and electrons at one radial position, one value of collisionality and one value of the electric field. Using the 4-D approach, to meet such requirements, one needs to undergo only two drift kinetic solver iterations per loop (one per species). However, in the monoenergetic approach, we need more iterations. As seen in the previous section, the monoenergetic coefficients are species-independent through the dependence on $\nu /v$ instead of $\nu /v_{th}^a$. Thus, to obtain the thermal transport coefficients $L_{ij}^a$ from the monoenergetic coefficients, one needs to perform an energy convolution of the monoenergetic coefficients for which we need to obtain $D_{ij}(v)$ for different values of $\nu /v$. We could simplify such a method by performing an interpolation for each species. Because in a reactor electrons are usually close to the $1/\nu ^*$ transport regime, the behaviour of the $D_{11}$ coefficients for the electrons can be approximated by a linear interpolation and thus at least two values of $D_{11}$ are necessary to be computed at $E^*_{\psi }=0$. For the ion species, however, these can be in different collisionality regimes such as the $1/\nu ^*$, $\nu ^*$ or $\sqrt {\nu ^*}$. Therefore, at least around five iterations of the solver would be needed to use interpolation methods. If we would choose the monoenergetic approximation to calculate the cost function, at each optimisation loop, a minimum of seven drift kinetic solver iterations would be needed, together with the interpolation and convolution steps. It is therefore more cost-effective to use the 4-D approach in the optimisation loop. This leads to the choice of the SFINCS code to calculate $J_{E_{\psi }}$ during optimisation. The thermal transport coefficients $L^a_{{11}_S}$ outputted by SFINCS are defined to be non-dimensional and thus they are slightly different from the definitions already presented. In terms of these SFINCS thermal transport coefficients, we may write the cost function in (3.6) as

(3.7)\begin{equation} J_{E_{\psi}} = \left(\frac{L^i_{{11}_S}}{L^e_{{11}_S} T_e^{{-}2}v_{th}^e} \right)^2= \left(\frac{L^D_{{11}_S}T_D^{{-}2}v_{th}^D+\gamma L^T_{{11}_S}T_T^{{-}2}v_{th}^T }{(1+\gamma)L^e_{{11}_S}T_e^{{-}2}v_{th}^e} \right)^2. \end{equation}

Furthermore, even if we are aiming for a reactor scenario, we may consider only deuterium as the ion species in calculating the electron root cost function. Note that the tritium species has a slightly smaller thermal velocity and thus its normalised electric field is larger than that of deuterium, for the same electric field value. Therefore, the radial transport coefficient of the tritium species will always be equal to or slightly smaller than its deuterium counterpart. This makes the weighted average ion species radial transport coefficient, (2.37), equal or smaller when considering both species than when considering deuterium alone. Therefore, for minimising the ratio of ion to electron radial transport coefficients in a reactor, the worst-case scenario occurs when the tritium species is absent and we may substitute the ion-weighted average radial thermal transport coefficient in (2.37) by the radial thermal transport coefficient of deuterium in the cost function ($\gamma =0$), and skip one drift kinetic solver iteration. The total cost function for electron root optimisation can then be defined as

(3.8)\begin{equation} J_{2} =\sum_{j=1}^3 w_{{E_{\psi}}_j}\left[ \frac{L^i_{{11}_S} (s=s_j)}{L^e_{{11}_S} (s=s_j)}\frac{T_e^2 v_{th}^i}{T_i^2 v_{th}^e}\right]^2, \end{equation}

in which we consider deuterium $D$ as the only ion species $i$ and we choose to optimise for at least three flux surfaces, in order to better target the electron root over a wider extent of the plasma radius. More surfaces could of course be targeted, but three surfaces were found to be a reasonable choice as we will see in the next sections.

Collisionality and electric field values are necessary as input for SFINCS to solve the drift kinetic equation. The collisionality values are defined by the temperature and density profiles of the deuterium and electron species. In fusion reactors, one would envisage to have an equithermal plasma with $T_e=T_i$. However, a clamp of the ion temperature is often observed in experiments (Beurskens et al. Reference Beurskens, Bozhenkov, Ford, Xanthopoulos, Zocco, Turkin, Alonso, Beidler, Calvo and Carralero2021) so that $T_i< T_e$, but such behaviour is not expected in a reactor plasma. Regardless, during the optimisation, we consider the worst-case scenario to achieve an electron root, which occurs when $T_e=T_i$ (Beidler et al. Reference Beidler, Drevlak, Geiger, Helander, Smith and Turkin2024) and thus we consider equal temperature profiles for ions and electrons. We also consider quasineutrality to be satisfied and thus $n_e=n_i$. To follow previous work on neoclassical transport optimisation (Landreman, Buller & Drevlak Reference Landreman, Buller and Drevlak2022) we use the following profiles for temperature:

(3.9)\begin{equation} T=(T_0-T_b)(1-s)+T_b, \end{equation}

and density

(3.10)\begin{equation} n=(n_0-n_b)(1-s^5)+n_b. \end{equation}

Here, $T_0$ and $n_0$ are the temperature and density on-axis. Here $T_b$ and $n_b$ are the temperature and density at the boundary, which allow the collisionality parameter $\nu ^*$ to remain finite at this position. We use an on-axis temperature $T_0=17.8\ \textrm {keV}$ and density $n_0=4.21\times 10^{20}\ \textrm {m}^{-3}$ yielding a volume-averaged plasma $\beta$ of $4\,\%{-}5\,\%$ for a magnetic field configuration that does not deviate much from $B=8\ \textrm {T}$ and $a=1.2\ \textrm {m}$, which are the reactor-like targets used for scaling the initial configuration used for the optimisations in this work. At the boundary, we choose the temperature $T_b=0.7\ \textrm {keV}$ and the density $n_b=0.6\times 10^{20}\ \textrm {m}^{-3}$ which are consistent with the electron root solution in Beidler et al. (Reference Beidler, Drevlak, Geiger, Helander, Smith and Turkin2024). While different methods could have been chosen to obtain the boundary values, different choices of these values are not seen to have a significant impact on the electron root solution optimisation which mostly takes place in the core region. We only introduce these finite boundary values so that at the testing step, we may use the same profiles as during the optimisation loop, as will be shown in the next section. The electric field profile is unknown at the optimisation stage and thus we have to make a guess to this input. With regards to the choice of electric field values, we want to optimise for a finite electric field value that ideally would be as large as possible in order to obtain a large electron root. In this work, we use the value of electric field values for the three radial coordinates considered such that $eE_{r}/T=1$, which yielded favourable results for every optimisation considered.

Finally, we discuss the three remaining sets of inputs for SFINCS. The first is which trajectory set to choose when using the 4-D approach. This approach allows the usage of any of the three types of trajectories discussed in the previous section. However, as we show in the next section, the testing of the optimisation is simplified by making use of the monoenergetic radial transport coefficients. The only set of trajectories compatible with the monoenergetic approach are the ‘DKES trajectories’ as discussed in the previous section. Thus, to be consistent between the optimisation and testing steps, we choose the ‘DKES trajectories’ during optimisation even if the 4-D approach is used. The same argument is used to discuss the second set of inputs which concerns the type of collision operator to use. A Lorentz collision operator is used for consistency with the testing method. Nevertheless, the absence of the momentum correction term and energy terms in the collision operator are not expected to significantly change the results as these terms only become important at higher collisionalities and we expect our reactor plasma species to be at low collisionality values. The third set of inputs concerns the resolution used during the optimisation. The convergence of drift kinetic solvers such as SFINCS and DKES is extremely dependent on the magnetic configuration. During the optimisation, the magnetic configuration is always changing, making it difficult to set a fixed resolution for convergence. In this work, we opted to use a resolution of $N_{\theta }=25$, $N_{\phi }=31$, $N_{\xi }=34$ and $N_x=4$. These are not enough to ensure complete convergence of SFINCS for every configuration. Yet, complete convergence of the radial coefficients is not strictly necessary for the electron root optimisation as we are interested in the minimisation of the ratio of the coefficients and we are not looking for an exact match of the value of radial particle flux of each species. Thus, the values of resolution are chosen as the minimum required resolution for a usual W7-X-like configuration in such a way that doubling the resolution for each coordinate will lead to less than a $6\,\%$ variation in the transport coefficients. As this is configuration-dependent, we leave here some general guidelines about the resolution parameters. At low collisionality, we expect $N_x=4$ and $N_{\theta }=25$ to be generally good values for convergence in SFINCS. At the same time, the number of points in the $\phi$ direction should usually be equal or larger than in the $\theta$ direction and the resolution in the pitch angle variable $\xi$ should be equal or larger than the resolution in $\phi$.

The optimisations presented in this work are optimised using the total cost function $J=J_{1}+J_{2}$ subject to the ‘trust region reflective’ optimisation algorithm. We use boundary Fourier coefficients RBC and ZBS as independent variables with $M_\textrm {pol}=N_\textrm {tor}=2$. We use 24 flux surfaces for VMEC resolution during optimisation, unless stated otherwise. The QI cost function is targeted at $4$ flux surfaces, namely $s=1/24$, $3/24$, $9/24$ and $13/24$. For the calculation of this cost function, $141$ points are measured along each well, with $27$ wells and $51$ bounce points per well. For the Boozer coordinate transformation, $18$ poloidal and $18$ toroidal mode numbers are used. A maximum mirror ratio of $0.19$ is targeted and the target aspect ratio is $10$. To avoid low-order rational surfaces we target a mean iota of $\bar {\iota }=0.61$. The electron root optimisation is targeted at the radial positions $\rho =\sqrt {s}=r/a=0.2$, $0.29$ and $0.35$ with $a$ the minor radius as defined by VMEC. Unity weights are used for all cost functions except when targeting the mirror ratio in which a weight of $w_{M}=100$ is used. The optimisation method is carried out with SFINCS integrated into the SIMSOPT framework. All the optimisation scripts used in this work, and the configurations and datasets obtained can be found in Lascas Neto et al. (Reference Lascas Neto, Jorge, Beidler and Lion2024).

4 Validation method

In order to measure the efficiency of the electron root cost function, it is necessary to obtain the electric field profile for the configurations obtained from the optimisation and check for the presence of an electron root. The 4-D approach in SFINCS is able to reconstruct the ion and electron fluxes as a function of the electric field at different radial positions and assess the presence of the solutions of the ambipolarity constraint, i.e. the electric field values at which the ion and electron particle-flux-density curves cross. However, such a method is tedious to apply for several radial positions, is local and does not provide any information about the transition zones between negative electric field values (ion root) and positive electric field values (electron root). Instead, we solve the electric field transport equation, (2.8), which holds information about ambipolarity and the ‘diffusion’ of the electric field. This becomes important near root transitions. A code that solves this equation is NTSS (Turkin et al. Reference Turkin, Beidler, Maaberg, Murakami, Tribaldos and Wakasa2011).

The NTSS code is capable of solving the one-dimensional radial transport equations, (2.9)–(2.10), for deuterium, tritium and electrons, and the electric field transport equation, (2.8). Helium ash can also be considered by solving the continuity equation, (2.10), for helium. The NTSS code requires as inputs the initial temperature and density of the species of interest, the magnetic configuration obtained from VMEC and converted in Boozer coordinates with BOOZ_XFORM, and a database of monoenergetic radial transport coefficients at different radial positions. As a remark, we note that if the magnetic configuration is fixed, a database of the monoenergetic coefficients, which are species-independent (and therefore independent of the temperature and density profiles), is more relevant for the calculation of the particle-flux-densities and energy-flux-densities as these can be calculated solely by interpolating the monoenergetic coefficients database for different species, collisionality values, electric field values and radial positions. This procedure is performed in NTSS in order to be able to solve the transport equations, in which the temperature and density can evolve over time. The time evolution of the temperature and density would render the use of the radial thermal transport coefficients much less efficient than in the case of the optimisation loop. Historically, NTSS obtains such a database from DKES. Therefore, all the input files make use of DKES conventions. Since our optimisation is performed with the SFINCS code, to be consistent, we would preferably use SFINCS as a generator of the monoenergetic transport coefficients database. To replicate the DKES-like database with SFINCS, we use the typical range of values of the normalised $E\times B$ velocity $E_r/(vB)$ and inverse mean free path $\nu /v$ inputs used for generating a DKES-like NTSS database. These values are chosen to be the ones used to obtain the DKES database for the magnetic configuration Hydra-Np04-20190108 in Beidler et al. (Reference Beidler, Drevlak, Geiger, Helander, Smith and Turkin2024), which contains an electron root. We convert these values to the corresponding values of the SFINCS inputs $\nu '$ and $E_{\psi }^*$ for the magnetic configuration of interest. We then perform a convergence study in order to find a reasonable resolution to describe the radial transport monoenergetic coefficients for the magnetic configuration of interest at the different electric field, collisionality and radial points. This is achieved by varying each resolution parameter by up to $150\,\%$ of the initial values and ensuring that the difference between the initial and final radial transport coefficient values varies $3\,\%$ or less. This convergence study is done for $E_{\psi }=0$ at the radial position which is nearest to the magnetic axis ($r/a=0.12247$), and for the lowest collisionality value of the database. This point in the database space is chosen because it is the point which is expected to need a larger resolution, and thus, it is the best point to ensure we are close to convergence at all points in the database. In practice, it is possible that this level of convergence will not be met at all points. However, doing such a convergence study for every database point is expensive and impractical. Therefore, after setting the resolution through this process and obtaining the database scan with that resolution, we check for any point that does not follow the expected trend in the scan and increase the resolution accordingly. Note that we are not interested in the convergence of $D_{31}$, $D_{13}$ and $D_{33}$ as they are not relevant for calculating the radial particle and energy-flux-densities of interest for this work if we neglect the collisional momentum correction which is expected to be small at the collisionalities of interest. As such, the resolution is based on the convergence of only the radial transport coefficient $D_{11}$ to reduce unnecessary computation time. After obtaining the monoenergetic radial coefficient database with SFINCS, we convert this database into a DKES-like database file which can be read as an input by NTSS. The typical DKES database file used in Beidler et al. (Reference Beidler, Drevlak, Geiger, Helander, Smith and Turkin2024) uses seven radial positions at $r/a=[0.12247,0.25,0.375,0.5,0.625,0.75,0.875]$, and we follow the same pattern in this work. Moreover, such a database calculation is applied for the optimised configurations with a radial VMEC resolution of $201$ which is larger than the resolution of $24$ points used during optimisation.

The use of SFINCS to obtain a DKES-like database requires two conversions of quantities due to the different normalisations used in the two codes. Such normalisations are detailed in the Appendix (A). To confirm that these conversions are being performed correctly, a benchmark against the DKES results for the electron root configuration Hydra-Np04-20190108 shown in Beidler et al. (Reference Beidler, Drevlak, Geiger, Helander, Smith and Turkin2024) is provided here. Figure 1 shows the comparison of the radial transport monoenergetic coefficients at two different radial positions between SFINCS and DKES. These results were obtained starting from the same DKES-like database set of inputs $\nu /v$ and $E_r/(vB)$. Note that due to historical reasons the $E_{r}/(vB)$ parameter in the DKES-like database files is written using the electric field parameter $E_{\tilde {r}}$ in which $\tilde {r}$ is defined in the Appendix (A). The quantity $\hat {\varGamma }_{11}$ is the monoenergetic radial transport coefficient normalised as per the conventions of an NTSS input (see the Appendix (A)). The resolution used in these scans was obtained from the convergence test as described in the previous section and is $N_{\theta }=39$, $N_{\phi }=59$ and $N_{\xi }=116$. We can see that the overall values of $\hat {\varGamma }_{11}$ match quite well between the two codes at the two radial positions. Especially, if we account that these results are plotted on a logarithmic scale and thus span several orders of magnitude. At low collisionality, for the larger values of the normalised electric field at $r/a=0.12247$ in figure 1, there is a mismatch between DKES and SFINCS. However, SFINCS produces data points closer to the expected $\sqrt {\nu ^*}-\nu ^*$ behaviours than DKES. In figure 2 a similar comparison is performed for the larger electric field values. These values of the electric field are the values in which a resonant behaviour is observed. This can be observed, for example, by a sudden change of the radial transport coefficient when compared with the curves for the other values of the electric field (especially at middle values of collisionalities in which the electric field has, for the smaller non-resonant values, no effect on the radial coefficient). This resonance is due to the $E\times B$ poloidal drift being large enough to compete with the term which describes motion along the field line. We can see that in the case of these electric fields (see figure 2), there is a large mismatch. We expect the two codes not to agree in this region. The most important physics near this resonance is described by the vanishing of the terms in the drift kinetic equation, and thus the numerical result depends on the exact calculation of numerical values very close to zero. Since the two codes use distinct numerical methods, it is expected that they will present a different behaviour in this region.

Figure 1. Comparison between SFINCS (crosses) and DKES (circles) monoenergetic radial transport coefficients (normalised to DKES conventions) at $r/a=0.12247$ and $r/a=0.5$ for the Hydra-Np04-20190108 configuration presented in Beidler et al. (Reference Beidler, Drevlak, Geiger, Helander, Smith and Turkin2024). Values are shown for the electric fields away from resonant values. A SFINCS resolution of $N_{\theta }=39$, $N_{\phi }=59$ and $N_{\xi }=116$ is used.

Figure 2. Comparison between SFINCS (crosses) and DKES (circles) monoenergetic radial transport coefficients at $r/a=0.12247$ and $r/a=0.5$ for the Hydra-Np04-20190108 configuration presented in Beidler et al. (Reference Beidler, Drevlak, Geiger, Helander, Smith and Turkin2024). Comparison is shown for large electric fields near resonant values. A SFINCS resolution of $N_{\theta }=39$, $N_{\phi }=59$ and $N_{\xi }=116$ is used.

Still, to better understand if this mismatch is inherent to the differences between the codes and not because we are using a resolution based on the lowest collisionality point with a zero electric field, we redo the convergence test using a point at the lowest collisionality value, but at a large normalised electric field value near the resonance for $r/a=0.12247$. Indeed, the convergence study at this point indicates a need for a larger resolution of $N_{\theta }=39$, $N_{\phi }=431$, $N_{\xi }=350$. Using this larger resolution, we redo the transport coefficients scans for these large electric fields presented in figure 2. The results are again compared with the DKES results in figure 3. We can see that for $r/a=0.12247$, the increase in resolution provides a small improvement but the two codes still present a mismatch, which agrees with the expectation that the behaviour of the resonance will be captured differently by the two codes. At $r/a=0.5$, we can see from comparing figures 2 and 3 that the increase of resolution has almost no effect and the mismatch between the codes remains. Nevertheless, these large resonant electric fields are only expected to be important for large mass impurity species and thus should not affect the analysis of our optimised configurations in this work. In fact, for a deuterium ion at $T=17.8\ \textrm {keV}$ the electric field would need to be of the order of $E_r\sim 100\ \textrm {kV}\ \textrm {m}^{-1}$ to experience the resonant behaviour seen in figures 2 and 3.

Figure 3. Comparison between SFINCS (crosses) and DKES (circles) monoenergetic radial transport coefficients at $r/a=0.12247$ and $r/a=0.5$ for the Hydra-Np04-20190108 configuration presented in Beidler et al. (Reference Beidler, Drevlak, Geiger, Helander, Smith and Turkin2024). Comparison is shown for large electric field values near resonance. A SFINCS resolution of $N_{\theta }=39$, $N_{\phi }=431$ and $N_{\xi }=351$ is used.

As a final check, we run NTSS with the DKES- and SFINCS-generated databases. We compare the two results of the electric field profile in figure 4, in which the SFINCS database for the resolution $N_{\theta }=39$, $N_{\phi }=59$ and $N_{\xi }=116$ is used. We see that the electron root produced by the SFINCS-generated database is slightly different. This difference is expected since there are small differences between the values of the monoenergetic coefficients as was seen in figure 1. However, such a difference in both the radial scans and the electron root is small. The DKES electron root has a maximum value of $E_{{r}_\textrm {DKES}}^{max}=22.521 \ \textrm {kV}\ \textrm {m}^{-1}$ and in the case of SFINCS, the maximum is $E_{{r}_\textrm {SFINCS}}^\textrm {max}=22.883 \ \textrm {kV}\ \textrm {m}^{-1}$. The electron root to ion root transition is calculated by taking the mean value of the two radial grid points in which the electric field observes a sign change. The root transition location is with DKES $r_\textrm {DKES}=0.867\ \textrm {m}$ and with SFINCS $r_\textrm {SFINCS}=0.830\ \textrm {m}$. Thus, the relative error of the maximum electric field calculation is of $1.607\,\%$ while the error of the root transition position is $-4.25\,\%$. The largest error is found to be in the transition zone due to the sensitivity of this region to small variations of the inputs. Nevertheless, the errors are relatively small given that the calculation of such a profile involves using the results of two drift kinetic numerical solvers with different numerical methods and interpolating such results in a three-dimensional parameter space before solving the transport equations. Therefore, despite these small differences between the two codes, we expect the validation method proposed to be appropriate to test the optimisations for the presence of an electron root.

Figure 4. Comparison of the radial electric field solution obtained with NTSS from using as an input the DKES and the SFINCS generated databases of the monoenergetic radial transport coefficients for the resolution $N_{\theta }=39$, $N_{\phi }=59$ and $N_{\xi }=116$.

5 Results

In this section, optimisation studies in order to assess the use of the electron root cost function proposed in (3.8) are performed. All the optimisations shown in this section will be done assuming a vacuum solution of the MHD force balance equation except when stated otherwise.

5.1 The initial configuration

The efficiency of the QI cost function (3.3), was shown in Goodman et al. (Reference Goodman, Camacho Mata, Henneberg, Jorge, Landreman, Plunk, Smith, Mackenbach, Beidler and Helander2023) to be highly dependent on the initial condition used. Therefore, we use one of the QI-optimised equilibria obtained in Goodman et al. (Reference Goodman, Camacho Mata, Henneberg, Jorge, Landreman, Plunk, Smith, Mackenbach, Beidler and Helander2023) as a starting point. We thus select the QI-optimised configuration obtained in that work for $n_\textrm {fp}=2$. From that configuration, we remove the boundary coefficients with $M>2$ and $N>2$ and optimise with maximum mode numbers $M_\textrm {pol}=2$ and $N_\textrm {tor}=2$. We also scale the minor radius to $a=1.2\ \textrm {m}$ and $B=8\ \textrm {T}$ while keeping the aspect ratio $A_\textrm {initial}=10$. The magnetic field strength of the resulting configuration is depicted in figure 5 at different flux surfaces. The magnetic strength contours in Boozer coordinates at $s=0.1025$ indicate some deviations from precise QI. This is due to the removal of the boundary harmonics with $M>2$ and $N>2$ from the original configuration. We then optimise this initial configuration using the QI cost function alone, the electron root cost function alone and both cost functions. This allows us to compare the different configurations that result from such optimisations. We note, however, that the mirror ratio, aspect ratio and mean iota targets will be maintained fixed in all of the different optimisations.

Figure 5. Contours of the magnetic field strength of the initial configuration at $s=0.1025$ in Boozer coordinates and the same quantity at the plasma boundary.

Regarding this initial condition, the monoenergetic radial transport coefficients scans at $r/a=0.12247$ and $r/a=0.5$ can be seen in figure 6. With the monoenergetic radial transport coefficient database, we can solve the electric field diffusion equation, (2.8), with the NTSS code. We note here that the electric field diffusivity $D_{E_{\psi }}^a$ is an input of the NTSS code. In this work, we follow Beidler et al. (Reference Beidler, Drevlak, Geiger, Helander, Smith and Turkin2024) and use $D_{E_{\psi }}^a=2\ \textrm {m}^2\ \textrm {s}^{-1}$ in every NTSS simulation. This parameter will act in the root transition to select a unique ambipolarity root. The transition is sharper or smoother if the $D_{E_{\psi }}^a$ is smaller or larger. The value of $D_{E_{\psi }}^a=2\ \textrm {m}^2\ \textrm {s}^{-1}$ was seen to be enough to deliver a smooth transition in the simulations performed. A more sophisticated way to obtain such profiles would be to solve the drift kinetic equation at higher orders in ways similar to Kuczynski et al. (Reference Kuczynski, Kleiber, Smith, Beidler, Borchardt, Geiger and Helander2024), but that method is computationally expensive and out of the scope of this paper. Maintaining the constant density and temperature profiles equal to the ones used during optimisation (see (3.9) and (3.10), with $T_0=17.8\ \textrm {keV}$ and $n_0=4.21\times 10^{20}\ \textrm {m}^{-3}$), and assuming equal concentrations of deuterium and tritium as pertain to reactor conditions, we obtain the electric field profile which can be seen in figure 7 where $a=1.191\ \textrm {m}$. The electric field solution is an ion root for the initial configuration at the temperature and density profiles used during optimisation.

Figure 6. Scans of the monoenergetic radial transport coefficients at $r/a=0.12247$ and $r/a=0.5$ for the initial magnetic configuration.

Figure 7. Radial electric field solution at the temperature and density profiles in (3.9) and (3.10), with $T_0=17.8\ \textrm {keV}$ and $n_0=4.21\times 10^{20}\ \textrm {m}^{-3}$, for the initial magnetic configuration. Solution calculated using the NTSS code by solving (2.8).

5.2 The QI optimisation

We now optimise the initial configuration regarding a QI target only. While performing this optimisation, we still maintain the previously discussed targets of the mirror ratio, mean iota and aspect ratio in the optimisation. The optimisation resulted in a final cost function of $3.003\times 10^{-5}$. The final aspect ratio is $10.000$ and the volume-averaged magnetic field is $7.962 \ \textrm {T}$ and a final mirror cost function of $2.90\times 10^{-6}$. The magnetic field strength can be seen in figure 8. We can see that the magnetic field strength in Boozer coordinates at $s=0.1025$ has no local minima, contrary to what was observed in figure 5 for the initial configuration. This indicates that the QI optimisation was successful in making the magnetic field closer to a QI stellarator. At other flux surfaces, the field still presents some deviations from QI. However, this is not surprising as the QI target is applied at only four flux surfaces $s=[1/24,5/24,9/24,13/24]$. It is thus expected that the field outside the exact location of these flux surfaces may be less close to a QI field. The monoenergetic radial transport coefficients scans for this optimised configuration at $r/a=0.12247$ and $r/a=0.5$ are shown in figure 9. At both radial positions we see that the plateau regime at intermediate collisionality values disappears after the QI optimisation when compared with the initial configuration (see figure 6). While the overall magnitude $1/\nu ^*$ regime is roughly the same at $r/a=0.12247$ between the two configurations, the reduction of transport for the same absolute value of the electric field $E_r$ is larger in the QI-optimised configuration due to the disappearance of the $1/\nu ^*$ regimes at finite $E_r$ and low collisionalities for the optimised configuration. At $r/a=0.5$ the value of the transport coefficients in the $1/\nu ^*$ regime is reduced (as expected for a magnetic field closer to perfect QI) even if at finite electric field we still see the presence of a $1/\nu ^*$-like regime. By solving the electric field diffusion equation with the temperature and density profiles used during optimisation (see (3.9) and (3.10) with $T_0=17.8\ \textrm {keV}$ and $n_0=4.21\times 10^{20}\ \textrm {m}^{-3}$) we obtain the electric field solution shown in figure 10 where $a=1.190\ \textrm {m}$. We see that an ion root is again observed for the density and temperature values chosen for the optimisations in this work. Thus, we conclude that the QI cost function alone is not enough to optimise the magnetic configuration to have an electron root solution at such plasma parameters. However, we note that the magnitude of the radial electric field solution obtained for this optimised configuration is smaller overall than the solution in figure 7 obtained for the initial configuration. This indicates that, for this configuration, ambipolarity is achieved for smaller values of the radial electric field strength. Such a result is a consequence of the reduction of the ion-relevant radial transport observed in this configuration.

Figure 8. Contours of the magnetic field strength at $s=0.1025$ in Boozer coordinates, for the magnetic configuration optimised only for the QI cost function and the same quantity at the plasma boundary.

Figure 9. Scans of the monoenergetic radial transport coefficients at $r/a=0.12247$ and $r/a=0.5$ for the magnetic configuration optimised only for the QI cost function.

Figure 10. Radial electric field solution at the temperature and density profiles in (3.9) and (3.10), with $T_0=17.8\ \textrm {keV}$ and $n_0=4.21\times 10^{20}\ \textrm {m}^{-3}$, for the magnetic configuration optimised only for the QI cost function. Solution obtained using the NTSS code by solving (2.8).

5.3 Electron root optimisation

We now analyse the result of an optimisation of the initial condition when we target only the electron root cost function while keeping the mirror, aspect ratio and mean iota targets the same. The final aspect ratio of this configuration is $9.985$, the final volume average magnetic field is $7.906 \ \textrm {T}$, the final total cost function is $0.079$ and the final mirror cost function is $2.67\times 10^{-3}$. In particular, the final electron root cost function at $r/a=0.2$ is $0.154$, at $r/a=0.29$ is $0.117$ and at $r/a=0.35$ is $0.142$. The magnetic field strength at different flux surfaces is shown in figure 11. We can see from the magnetic field strength contours in Boozer coordinates at $s=0.1025$ that the resulting magnetic field has a preference for poloidal closed contours. However, the resulting magnetic field is quite non-optimised in regards to omnigenity. Indeed, looking more closely at the magnetic field strength at $s=0.1025$, we can observe that the magnetic spectrum has three local wells per field period and that these wells have different shapes for different field lines.

Figure 11. Contours of magnetic field strength at $s=0.1025$ in Boozer coordinates for the magnetic configuration optimised only for the electron root cost function, and the same quantity at the plasma boundary.

The monoenergetic radial transport coefficients scans at different radial positions can be seen in figure 12. We can see that there is a large increase of the overall magnitude of the $1/\nu ^*$ regime when compared with the QI-optimised configuration (see figure 9). However, at finite electric field and low collisionality, there are well-defined transitions from a $1/\nu ^*$ to a linear $\nu ^*$ regime, with the $\sqrt {\nu ^*}$ regime disappearing almost completely. The electric field profile obtained from solving the electric field transport equation is shown in figure 13. The temperature and density profiles used to obtain such a solution were the same as the ones used during optimisation, which we recall here, are given by (3.9) and (3.10), with $T_0=17.8\ \textrm {keV}$ and $n_0=4.21\times 10^{20}\ \textrm {m}^{-3}$. It is observed that this configuration has an electron root at these temperature and density profiles which is the original purpose of the electron root cost function. Nevertheless, as we do not target any symmetry constraints, the magnetic field can freely adapt to satisfy the electron root function and the final result is not a QI field as shown in figure 11. The electron root for this configuration has a maximum of $E_{r}^\textrm {max}=28.557\ \textrm {kV}\ \textrm {m}^{-1}$ and a root transition at $r/a=0.570$ with $a=1.195\ \textrm {m}$.

Figure 12. Scans of the monoenergetic radial transport coefficients at $r/a=0.12247$ and $r/a=0.5$ for the magnetic configuration optimised only for the electron root cost function.

Figure 13. Radial electric field solution at the temperature and density profiles in (3.9) and (3.10), with $T_0=17.8\ \textrm {keV}$ and $n_0=4.21\times 10^{20}\ \textrm {m}^{-3}$, for the magnetic configuration optimised only for the electron cost function. Solution obtained using the NTSS code by solving (2.8).

5.4 The QI and electron root optimisation

We now optimise for both QI and electron root targets keeping the same aspect ratio, mean iota and mirror ratio targets as before. The final aspect ratio of the optimised configuration is $10.015$, the volume-averaged magnetic field strength is $8.373\ \textrm {T}$ and the final total cost function is $0.043$. The final electron root cost function is $0.130$ at $r/a=0.2$, $0.0917$ at $r/a=0.29$ and $0.120$ at $s=0.35$ and the final mirror cost function is $2.79\times 10^{-3}$. The magnetic field strength contours in Boozer coordinates at $s=0.1025$ and the magnetic field strength at the plasma boundary are depicted in figure 14. At $s=0.1025$ we can see that the magnetic field strength presents the usual contours of a QI configuration with a small local minimum. When comparing such contours with the previous ones shown for the QI-only optimised configuration in figure 8, it is clear that adding the optimisation for the electron root adds this local minimum, making the magnetic field less close to perfect QI. Nevertheless, such defects at the finer structure of the magnetic field can be tuned by changing the weights of both QI and electron root cost functions, which were maintained at unity throughout this work. The monoenergetic radial transport coefficients scans at $r/a=0.12247$ and $r/a=0.5$ for this configuration are shown in figure 15. From the radial transport coefficient scans we see again the increase of the $1/\nu ^*$ transport when compared with the configuration optimised only for QI (see figure 9). However, such an increase is not as large as in the case of the optimisation performed only for an electron root (see figure 12). At the same time, the behaviour at finite electric field is again characterised by a well-defined and quick transition from the $1/\nu ^*$ regime to a $\nu ^*$ regime at both radial positions. As in the previous optimisation, which targeted only the electron root cost function, such behaviour is a consequence of the disappearance of the $\sqrt {\nu ^*}$ transport regime at low collisionality and finite electric field. In this case, though, the $\nu ^*$ regime at low collisionality for the scan at $r/a=0.5$ is better defined than in the previous optimisation done only for an electron root (see figure 12). We solve again the electric field transport equation while keeping the same temperature and density profiles used during optimisation, which we again recall are given by (3.9) and (3.10), with $T_0=17.8\ \textrm {keV}$ and $n_0=4.21\times 10^{20}\ \textrm {m}^{-3}$. The electric field solution obtained is depicted in figure 16. We can see that even with a configuration with relatively good QI properties we can obtain an electron root at reactor parameters. This electron root has a maximum of $E_r^\textrm {max}=33.791\ \textrm {kV}\ \textrm {m}^{-1}$ and a root transition at $r/a=0.6300$ with $a=1.178\,\textrm {m}$. Notice also that the maximum of the electron root obtained in this case is larger than the one obtained when optimising only for the electron root (see figure 13). In fact, the final electron root cost function obtained was overall larger when optimising only for an electron root than in this case. These results indicate that the electron root optimisation can get better results when used together with the QI optimisation. Similar observations were made in Kim et al. (Reference Kim, Buller, Conlin, Dorland, Dudt, Gaur, Jorge, Kolemen, Landreman and Mandell2024) for the case of turbulence optimisation and in Bader et al. (Reference Bader, Drevlak, Anderson, Faber, Hegna, Likin, Schmitt and Talmadge2019) for the optimisation of fast particles. In the case of the electron root optimisation, such behaviour may occur because, in the optimisation using only the electron root cost function, a path that increases the $1/\nu ^*$ very quickly is first chosen by the optimiser, leading to a magnetic field strength with a lack of any symmetry (see figure 11). In the case in which both QI and electron root cost functions are used, the optimiser first selects a magnetic field with good QI properties limiting the increase of the $1/\nu ^*$ regime. Then, in both optimisations small adjustments are done to the magnetic field, in order to generate a well-defined $\nu ^*$ regime at low collisionality. From comparing figures 12 and 15, these small adjustments look easier to achieve by the optimisation loop when starting from a configuration closer to the QI magnetic geometry.

Figure 14. Contours of the magnetic field strength in Boozer coordinates at $s=0.1025$ for the magnetic configuration optimised for both QI and the electron root. The same quantity at the plasma boundary is also shown.

Figure 15. Scans of the monoenergetic radial transport coefficients at $r/a=0.12247$ and $r/a=0.5$ for the magnetic configuration optimised for both the electron root and QI.

Figure 16. Radial electric field solution at the temperature and density profiles in (3.9) and (3.10), with $T_0=17.8\ \textrm {keV}$ and $n_0=4.21\times 10^{20}\ \textrm {m}^{-3}$, for the magnetic configuration optimised for both QI and the electron root. Solution obtained using the NTSS code by solving (2.8).

In order to check the robustness of this optimised magnetic configuration with respect to the presence of an electron root for different plasma parameters we solve the electric field diffusion equation for different temperatures and density values. The profiles used are obtained by scaling the profiles in (3.9) and (3.10) used during the optimisation. We plot the results of this scan for the maximum of the electric field and the location of the root transition (normalised to the minor radius). The results can be seen in figure 17. We note that when the electric field has a maximum of zero, no electron root is present. At the same time, a zero value is attributed to the relative position of the transition to indicate the cases in which no electron root is present. As expected, the region of plasma parameters in which the electron root ceases to exist lies at low temperatures and high densities (delimited by the darkest blue contour), which represents the higher collisionality region. Nevertheless, the presence of the electron root for this optimised configuration is seen to be robust as it is present for a large number of temperatures and densities relevant to reactor scenarios. In particular, for a density on-axis of $n_0=4.21 \times 10^{20} \ \textrm {m}^{-3}$ we can see that for the entire range of temperatures scanned ($14.24{-} 21.36 \ \textrm {keV}$) the electron root is present. The trend in the scan also indicates that larger temperatures than $21.36\ \textrm {keV}$ (which were not simulated here but can also be relevant for reactor scenarios) should have an electron root at this on-axis density. The maximum of the radial electric field strength increases with increasing temperature and decreasing density, with the maximum values of $E_{r}^\textrm {max}=50.789 \ \textrm {kV}\ \textrm {m}^{-1}$ being achieved at the largest temperature of $21.36 \ \textrm {keV}$ and lower density of $1.05\times 10^{20} \ \textrm {m}^{-3}$ in the scan. More surprising perhaps is the position of the root transition, which from observing the scan in figure 17, increases more with decreasing density than with increasing temperature. At the lowest densities in the scan, the electron–ion root transition occurs quite near the plasma boundary, with the electron root spanning through the majority of the minor radius extension.

Figure 17. Scans of the maximum of the radial electric field and the location of the root transition for different temperatures and densities on-axis. The magnetic field and volume of the plasma are maintained fixed in this scan. The plotted line is a contour of constant $\beta =4.24\,\%$.

A similar scan was done by scaling the major radius and the magnetic field. In this scan, the aspect ratio is maintained constant by scaling the minor radius with the major radius, and $\beta =4.24\,\%$ is also kept constant by scaling the density and temperature profiles by the same scaling factor applied to the magnetic field. The resulting scans are shown in figure 18. We see that the area in which the electron root ceases to exist is located in the region of large major radius and low magnetic field (delimited by the darkest blue contour). It is observed that in the range of scanned major radius and magnetic field values, the electron root exists at almost all parameters. It is interesting to see that the maximum of the electron root has its larger value at a high magnetic field and small major radius values, with a maximum electric field of $E_{r}^\textrm {max}=69.307 \ \textrm {kV}\ \textrm {m}^{-1}$ being obtained. The position of the root transition increases more with decreasing major radius than with the increasing magnetic field, and a value of up to $r/a=0.810$ can be obtained showing that a large region of the plasma can achieve an electron root for a small major radius plasma at a wide range of magnetic field values. Both the density and temperature scans in figure 17, and the radius and magnetic field scans in figure 18 show the possibility of an electron root solution at a wide range of both plasma and configuration parameters, which indicates the robustness of the optimised magnetic configuration to the presence of an electron root. It is also observed that the electron root maximum and the position of the root transition increase with the magnetic field. However, one would expect the radial transport coefficient in the $1/\nu ^*$ regime to decrease with the square of the magnetic field strength and the $\nu ^*$ regime radial transport coefficient to be independent of $B$ (Beidler et al. Reference Beidler, Allmaier, Isaev, Kasilov, Kernbichler, Leitold, Maaßberg, Mikkelsen, Murakami and Schmidt2011). These scalings are expected to make the electron root disappear as the magnetic field increases. To better investigate the presence of such behaviour, a scan in the values of the magnetic field and major radius allowing $\beta$ to vary is performed and can be seen in figure 19. This is achieved by maintaining the on-axis density and temperature fixed at the same values used in the optimisations, i.e. $T_0=17.8\ \textrm {keV}$ and $n_0=4.21\times 10^{20}\ \textrm {m}^{-3}$. We point out here that in these two last scans in figure 19, $\beta$ varies between small $\beta =1.88\,\%$ and quite large $\beta =16.96\,\%$ values. The larger beta values obtained for the lower magnetic fields are outside of the expected values for a stellarator fusion reactor and thus we do not expect these to be achievable in reality. Nevertheless, we add this scan for completeness and to see the behaviour of the electron root when the magnetic field is increased without adjusting the density and temperature profiles.

Figure 18. Scans of the maximum of the radial electric field and the location of the root transition for different scaling of the major radius and magnetic field strength on-axis. The plasma $\beta =4.24\,\%$ is constant for the entire scan. The plotted line is a contour of constant toroidal flux at the boundary.

Figure 19. Scans of the maximum of the radial electric field and the location of the root transition for different scaling of the major radius and magnetic field strength on-axis. Contrary to the scans in figure 18, the plasma $\beta$ varies in these scans. The plotted line is a contour of fixed toroidal flux at the boundary.

In the scan shown in figure 19, we see that the electron root ceases to exist in the region of large magnetic field and large major radius. Accordingly, it is observed that the increase of the magnetic field at a large major radius decreases the values of the maximum of the electron root. However, at small major radius values an increase in the magnetic field still increases the electron root maximum for the values scanned here. We note that for the smaller major radius, the normalised position of the root transition decreases with the magnetic field. Such behaviour indicates that at smaller sizes of the machine, the region of the plasma in which the electron root exists will get smaller with increasing magnetic field even if the electron root maximum may get large. Thus, at some high magnetic field, larger than the ones scanned here, the electron root transition will eventually disappear. The same will happen at a sufficiently large major radius. The electron roots with a larger position of the transition are located at a low major radius and low magnetic field. These scans retrieve the expected behaviour of the electron root being more difficult to achieve with an increasing magnetic field. Nevertheless, for this electron root optimised configuration, an electron root is still achievable at a value of magnetic field of $B=11\ \textrm {T}$ for a wide range of major radii. We point out that such magnetic field value is already quite large for state-of-the-art technologies. Outside of these scans, for the lowest value of the major radius shown here, it was found that the electron root ceases to exist between $B=37\ \textrm {T}$ and $B=38\ \textrm {T}$ which are unrealistic values for the currently available technologies.

5.5 Finite-$\beta$ optimisation

Up to now, the optimisations shown have considered vacuum VMEC solutions only. This does not account for the finite-$\beta$ corrections from the plasma pressure in the MHD equation. Such considerations are worth investigating because the electron root optimisation depends on the solution of the drift kinetic equation, which is obtained by considering the existence of finite prescribed density and temperature profiles. A pressure gradient is added to the MHD force balance, and we optimise for both QI and electron root targets with the pressure profile given by $P=2\ \textrm {nT}$, where the density and temperature profiles are again the ones in (3.9) and (3.10), with $T_0=17.8\ \textrm {keV}$ and $n_0=4.21\times 10^{20}\ \textrm {m}^{-3}$. We note that the optimisation for this case is done using $51$ flux surfaces in VMEC instead of the standard $24$ flux surfaces because it was observed that a smaller resolution would lead to an optimised configuration which would present significant variations in the $\iota$ profile when the final optimised configuration radial resolution was increased from $24$ to $201$ flux surfaces. Here we chose to maintain $s=[1/24,5/24,9/24,13/24]$ as the flux surfaces used for the QI optimisation even though we increased the number of flux surfaces to $51$. This choice means that we enforce the QI constraint to a higher degree in the core which may result in more competition between the QI cost function and the electron root cost function. The resulting configuration has an aspect ratio of $10.004$ and a volume-averaged magnetic field of $7.760\ \textrm {T}$. The final total cost function is $0.053$ and the final electron root cost function is $0.113$ at $r/a=0.2$, $0.116$ at $r/a=0.29$ and $0.1452$ at $r/a=0.35$. The final mirror cost function is $3.81\times 10^{-3}$. The magnetic field strength in Boozer coordinates at $s=0.1025$ and magnetic field strength at the plasma boundary for this configuration can be seen in figure 20. By comparing the contours of the magnetic field strength at $s=0.1025$ with the one obtained with the optimisation at a vacuum (see figure 14), we can see that both have similar behaviours, with the finite beta contours also showing a local minimum. The monoenergetic radial transport coefficient scans at $r/a=0.12247$ and $r/a=0.5$ can be seen in figure 21. From solving the electric field diffusion equation with the same profiles as in (3.9) and (3.10) used at optimisation, with $T_0=17.8\ \textrm {keV}$ and $n_0=4.21\times 10^{20}\, \textrm {m}^{-3}$, we obtain the electric field profile which is shown in figure 22. We can see that an electron root is achievable when considering a consistent finite-$\beta$ correction in the MHD force balance during optimisation. This electron root is smaller than the one obtained at the same temperature and densities with the vacuum optimisation (see figure 16). This result can be a consequence of the change of the number of flux surfaces while maintaining the same values of $s$ at which the QI cost function is targeted, which can result in more direct competition of the QI and electron root cost functions. The maximum of the electron root is in this case $E_{r}^\textrm {max}=22.569\ \textrm {kV}\ \textrm {m}^{-1}$ and the root transition occurs at $r/a=0.530$ with $a=1.200\ \textrm {m}$. This result shows that an electron root optimisation is possible with finite-$\beta$ corrections, especially when we consider that this optimised configuration has a volume averaged $\beta \approx 5\,\%$.

Figure 20. Contours of the magnetic field strength in Boozer coordinates at $s=0.1025$ for the magnetic configuration optimised for both QI and electron root with finite-$\beta$ corrections. The same quantity at the plasma boundary is also shown.

Figure 21. Scans of the monoenergetic radial transport coefficients at $r/a=0.12247$ and $r/a=0.5$ for the magnetic configuration optimised for both electron root and QI with finite-$\beta$ effects.

Figure 22. Radial electric field solution at the temperature and density profiles in (3.9) and (3.10) with $T_0=17.8\ \textrm {keV}$ and $n_0=4.21\times 10^{20}\ \textrm {m}^{-3}$ for the magnetic configuration optimised for both QI and electron root with finite-$\beta$ effects. Solution obtained with the NTSS code by solving (2.8).

5.6 Confinement properties

We now take a look at the two common figures of merit used to evaluate the confinement of fast particles, the effective ripple (epsilon effective), $\varepsilon _\textrm {eff}$, and the fraction of lost alpha particles from fusion reactions. The epsilon effective can be inferred from the slope of the $1/\nu ^*$ previous monoenergetic radial scans in figures 6, 9, 12, 15 and 21, which was discussed in the previous sections. Nevertheless, here we show the value of $\varepsilon _\textrm {eff}$ obtained with the code NEO (Nemov et al. Reference Nemov, Kasilov, Kernbichler and Heyn1999). The loss fraction of alpha particles is obtained with the code SIMPLE (Albert, Kasilov & Kernbichler Reference Albert, Kasilov and Kernbichler2020). This code does not include the effects of a strong electric field, but such an effect is expected to be small for the fast alpha particles due to their large velocity. We note that we did not optimise directly for either of these metrics, only indirectly using the QI cost function, as the main objective was to isolate and understand the effects and effectiveness of the electron root cost function. Indeed, as we observe in figure 23, we can see that for both vacuum and finite-$\beta$ configurations optimised for both QI and the electron root, the epsilon effective is less than $0.05$ for all the flux surfaces, but is around three times larger than the same quantity for W7-X. However, the loss fraction of fast particles saturates at a lower level for both vacuum and finite-$\beta$ configurations when compared with the vacuum standard W7-X configuration. The loss fraction saturates around $12\,\%$ for the vacuum configuration and, while not saturating, it is around $3\,\%$ for the finite-$\beta$ configuration at $t=0.01\ \textrm {s}$. The decrease in the loss of fast particles of the finite-$\beta$ configuration when compared with the vacuum configuration is probably related to the effect of the tangent magnetic fields which originate from the current due to the finite pressure gradient. We conclude that even if we do not achieve zero losses for precise QI-only optimised configurations, we still obtain a relatively small loss of fast particles compared with W7-X. We also note that the configuration optimised only for the QI cost function shows values of epsilon effective smaller than the ones observed for W7-X and the loss of fast particles is also the smallest of all the configurations.

Figure 23. Epsilon effective and loss fraction of fusion alpha particles for the vacuum standard W7-X configuration, for the configuration optimised only for QI and for the vacuum and finite-$\beta$ configurations optimised for both QI and electron root. For the loss fraction of alpha particles, the volume-averaged magnetic field and minor radius of the configurations were scaled to $5.7\ \textrm {T}$ and $1.7\ \textrm {m}$, respectively, to mimic the ARIES-CS scaling (Najmabadi et al. Reference Najmabadi, Raffray, Abdel-Khalik, Bromberg, Crosatti, El-Guebaly, Garabedian, Grossman, Henderson and Ibrahim2008). The alpha particles are initialised at $s=0.25$.

5.7 Power balance considerations

So far, we have calculated the electric field solution by solving only the electric field diffusion equation, (2.8), which means that we are solving NTSS transport equations with density and temperature profiles kept constant in time so that $\partial n_aT_a/ \partial t = 0$ in (2.9). It is also important to understand if an electric field can be realised while achieving the required target of fusion power for a fusion reactor. While an extensive study of such conditions is not the main topic of this work, we would like to understand if such solutions of the transport equations can be achieved with the optimised configurations for the QI and electron root optimised configurations shown in the previous sections. This is of particular importance if we take into account that we did not optimise for a small $\varepsilon _\textrm {eff}$. This means that the $1/\nu ^*$ regime, in which the heat flux scales with $T^{9/2}$, is not directly minimised in these optimisations. While we expect an electron root in the core region to reduce the ion transport due to the effect of the electric field in bounding the radial drift (Beidler et al. Reference Beidler, Allmaier, Isaev, Kasilov, Kernbichler, Leitold, Maaßberg, Mikkelsen, Murakami and Schmidt2011) as previously discussed in § 2, it is still true that the $1/\nu ^*$ regime is important in the electron to ion root transition zone and possibly in the ion root region if the electric field there is small. Therefore, it is important to verify if, at the reactor temperatures we are aiming for, we can still contend with the heat flux generated by the possible $1/\nu ^*$ regime in the root transition region.

With this goal in mind, we perform simulations with the NTSS code for the configuration optimised only for QI and for the vacuum and finite-$\beta$ configurations optimised for both QI and electron root, to target an alpha power close to $300\ \textrm {MW}$. We perform these simulations by fixing only the density profiles of deuterium and tritium while allowing the temperatures of all species and the electric field to evolve. We also allow the helium density to evolve as in the simulations performed in Beidler et al. (Reference Beidler, Drevlak, Geiger, Helander, Smith and Turkin2024), with a relative concentration of deuterium and tritium of $0.4$ and a helium relative concentration of $0.1$. Following Beidler et al. (Reference Beidler, Drevlak, Geiger, Helander, Smith and Turkin2024) we also add anomalous transport via a model in which $\chi _a=\mathcal {C}(P \textrm {[MW]})^{3/4}/(n [10^{20}\ \textrm {m}^{-3}])$, with the choice of $\mathcal {C}=6.5\times 10^{-3}\ \textrm {m}^2\ \textrm {s}^{-1}$ and $\mathcal {D}_a=0.01\chi _a$. We note that this is a very simple model but it is added here to capture the same conditions used for the simulation shown in Beidler et al. (Reference Beidler, Drevlak, Geiger, Helander, Smith and Turkin2024) in which an electron root was observed, thus allowing for a better comparison. For the configuration optimised only for QI we use (3.9) and (3.10) scaled to $T_0=7.12\ \textrm {keV}$ and $n_0=1.47\times 10^{20}\ \textrm {m}^{-3}$ as the initial profiles. These values are obtained by performing several simulations in which the initial temperature is varied (by scaling $T_0$) at a fixed density $n_0$, until the final solution provides the required target power. We obtain a plasma with $302\ \textrm {MW}$ of volume-averaged alpha power after subtracting the Bremsstrahlung losses and an energy confinement time of $\tau _E=1.68\ \textrm {s}$ is obtained which is a factor of $2.32$ above the ISS04 ( International Stellarator Scaling ) scaling (Yamada et al. Reference Yamada, Harris, Dinklage, Ascasibar, Sano, Okamura, Talmadge, Stroth, Kus and Murakami2005). The temperature and density profiles for the different species can be observed in figure 24. The radial electric field solution obtained is shown in figure 25 and is an ion root despite the low density and high temperatures of the ion species. Such a result is expected as this configuration was not optimised for electron root. The neoclassical and total (neoclassical and anomalous) energy fluxes, normalised to the volume-averaged alpha power $P_{\alpha }=302\ \textrm {MW}$, are depicted in figure 26. The neoclassical energy fluxes are seen to be small for both ion and electron species as expected from a QI-optimised configuration, with the sum of ion and electron neoclassical energy fluxes staying below $0.15P_{\alpha }$. The sum of ion and electron total energy fluxes (neoclassical and anomalous) stays below the value of the volume-averaged alpha power, except at the boundary, at which, the large energy flux variation is likely due to the artificial choice of the temperature and density values at the boundary.

Figure 24. Density (a) and temperature (b) profiles obtained with the NTSS code by performing a power balance simulation for the configuration optimised only for QI.

Figure 25. Radial electric field profile obtained from the NTSS code by performing a power balance simulation for the configuration optimised only for QI.

Figure 26. Neoclassical (a) and total (b) energy fluxes obtained for the power balance simulation performed for the magnetic configuration optimised only for QI. The simulation was performed by the NTSS code and the energy fluxes are normalised to the volume-averaged alpha power of $P_{\alpha }=302\ \textrm {MW}$.

We now perform a similar simulation for the vacuum configuration optimised for both QI and electron root. We use again the profiles shape given in (3.9) and (3.10) for the initial profiles. In this case, such profiles are scaled so that the initial density on-axis is $n_0=2.74\times 10^{20} \ \textrm {m}^{-3}$ and the initial temperature on-axis for all species is $T_0=12.46\ \textrm {keV}$. We obtain an equilibrium plasma with $305\ \textrm {MW}$ of volume-averaged alpha power after subtracting the Bremsstrahlung losses. An energy confinement time of $\tau _E=1.38\ \textrm {s}$ is obtained which is a factor of $1.50$ above the ISS04 scaling. The final density and temperature profiles can be seen in figure 27. The electric field profile obtained is depicted in figure 28, where we observe that an electron root is obtained while achieving a target alpha power close to $300\ \textrm {MW}$. This electron root has a maximum of $E_r^\textrm {max}=40.297\ \textrm {kV}\ \textrm {m}^{-1}$ and the root transition is located at $r/a=0.410$ with $a=1.178\ \textrm {m}$. We note that such an electron root is obtained for lower temperature and higher density values (and thus higher collisionalities) than the ones obtained in figure 24 for the configuration optimised only for QI. This shows the effect of the electron root cost function in helping to achieve an electron root, even when the power balance is considered. The neoclassical and total (neoclassical and anomalous) energy fluxes normalised to the volume-averaged alpha power $P_{\alpha }=305\ \textrm {MW}$ are shown in figure 29. We observe that inside the electron root region, there is a large reduction of the neoclassical energy fluxes, as expected. This reduction leads to the sum of ion and electron neoclassical energy fluxes being below $0.2P_{\alpha }$ in the electron root region. We point out that such values of neoclassical energy fluxes are very close to the maximum values of energy fluxes shown in figure 26 for the configuration optimised only for QI. Such observation is quite interesting as it shows that such reduced energy fluxes in the electron root zone are obtainable for this configuration, despite the fact that it has an overall much larger $\varepsilon _\textrm {eff}$ than the configuration optimised only for QI (see figure 23). Additionally, such values of energy flux obtained in the electron root are also below the maximum values observed for W7-X standard and high mirror configurations in Beidler et al. (Reference Beidler, Smith, Alonso, Andreeva, Baldzuhn, Beurskens, Borchardt, Bozhenkov, Brunner and Damm2021). Such a result indicates that if we are able to obtain a configuration with an electric field solution such that an electron root exists and spans over a large region of the plasma minor radius, we may be able to obtain reduced levels of neoclassical energy transport without the need of targeting the $\epsilon _\textrm {eff}$. The total energy fluxes for this configuration are again below $P_{\alpha }$, except at the boundary. This happens despite the neoclassical energy flux being larger in the ion root region than for the configuration optimised only for QI. This is just a consequence of the $300\ \textrm {MW}$ solution for the QI-only optimised configuration (see figure 24) being obtained at a smaller density, which makes the anomalous transport larger for the simple anomalous transport model considered. In the root transition, the energy fluxes have a small discontinuity due to the abrupt variation of the electric field. The maximum value of energy flux at the discontinuity is around $0.6P_{\alpha }$ which is still close to $0.5P_{\alpha }$ despite the large value of the epsilon effective seen in figure 23.

Figure 27. Density (a) and temperature (b) profiles obtained from a power balance simulation, performed by the NTSS code for the vacuum configuration optimised for both QI and electron root.

Figure 28. Radial electric field profile obtained from a power balance simulation, performed with the NTSS code for the vacuum configuration optimised for both QI and electron root.

Figure 29. Neoclassical (a) and total (b) energy fluxes obtained for the power balance simulation for the vacuum magnetic configuration optimised for both QI and electron root. The simulation was performed by the NTSS code and the energy fluxes are normalised to the volume-averaged alpha power of $P_{\alpha }=305\ \textrm {MW}$.

A similar simulation is now performed using the finite-$\beta$ configuration optimised for both QI and electron root. In this case, because we optimised for a specific pressure profile, using density and temperature profiles that do not match such pressure is not fully consistent. However, due to the aspect ratio and volume of the initial configuration, using the density and temperature profiles consistent with optimisation will lead to a very large alpha power that will be unrealistic from the engineering point of view. Moreover, even if NTSS is initialised with consistent temperature and density profiles, it will not lead to final profiles that are consistent with the pressure profile used in the MHD force balance. The most correct approach would be to obtain several optimised configurations for increasing $\beta$ and use them in an equal number of NTSS simulations to mimic a startup phase of the plasma, but such a task is out of the scope of this paper. We thus resort to the same strategy as in the previous vacuum configurations. We use again as initial profiles the density and temperature profiles of (3.9) and (3.10) scaled to an on-axis temperature of $T_0=12.55 \ \textrm {keV}$ and density of $n_0=2.31\times 10^{20} \ \textrm {m}^{-3}$. We obtain a volume-averaged alpha power of $P_{\alpha }=314\ \textrm {MW}$ after subtracting the Bremsstrahlung losses and an energy confinement time of $\tau _E=1.45\ \textrm {s}$ which is a factor of $1.79$ above the ISS04 scaling. The final temperature and density profiles obtained in this simulation are presented in figure 30 and the electric field solution can be seen in figure 31. In this case, we see that an electron root region is again present in the core region, but two ion root regions are present. One ion root region is the usual region spanning from the boundary of the plasma until the root transition occurs at $r=0.564 \ \textrm {m}$. Another one is an ion root living between the magnetic axis and the root transition at $r=0.228 \ \textrm {m}$. The electron root has a maximum of $E_r^\textrm {max}=28.577\ \textrm {kV}\ \textrm {m}^{-1}$. The electron root and ion root regions situated between $r=0.228 \ \textrm {m}$ and the plasma boundary are similar to the ones observed in figure 28. The ion root zone situated between the magnetic axis and $r=0.228 \ \textrm {m}$ has very small (almost residual) electric field values and is related to a flattening of the temperature and electron density profiles in this region (see figure 30). In particular, due to the flattening of such profiles in the first ion root region, there is a reduction in the magnitude of the peak of the electron density profile (see figure 30) with a positive electron density gradient occurring in the region of the electron root situated between $r=0.228 \ \textrm {m}$ and $r=0.564 \ \textrm {m}$. Such behaviour shares similarities with LHD experimental scenarios analysed in Fujiwara et al. (Reference Fujiwara, Kawahata, Ohyabu, Kaneko, Komori, Yamada, Ashikawa, Baylor, Combs and deVries2001), Fujita et al. (Reference Fujita, Satake, Nunami, García-Regaña, Velasco and Calvo2021) and Velasco et al. (Reference Velasco, Calvo, Satake, Alonso, Nunami, Yokoyama, Sato, Estrada, Fontdecaba and Liniers2016). The final neoclassical and total (neoclassical and anomalous) energy fluxes normalised to the volume-averaged alpha power $P_{\alpha }=314\ \textrm {MW}$ can be seen in figure 32. Two discontinuity regions can be seen in the energy flux solutions at the location of the two root transitions situated at $r=0.228 \ \textrm {m}$ and $r=0.564 \ \textrm {m}$ as expected. We observe again that in the electron root regions, the neoclassical energy fluxes are reduced with this case showing overall values of the sum of electron and ion neoclassical energy fluxes below $0.25P_{\alpha }$ in the electron root regions. In the root transition situated at $r=0.564 \ \textrm {m}$, the sum of ion and electron neoclassical energy fluxes reaches a value of around $0.6P_{\alpha }$ as seen for the previous case of the vacuum configuration optimised for QI and electron root (see figure 29). The total energy fluxes have a similar behaviour as in the previous simulation staying below $P_{\alpha }$ except at the boundary. We conclude this section by emphasising that despite the large values the $\varepsilon _\textrm {eff}$ observed for the electron root optimised configurations, the maximum heat flux in the electron root region of these configurations is around the maximum level of heat flux observed for the QI-only optimised configuration. This is an important observation as it indicates that using $\varepsilon _\textrm {eff}$ as an optimisation metric is not strictly necessary if we are looking for an electron root plasma operation.

Figure 30. Density (a) and temperature (b) profiles obtained with the NTSS code by performing a power balance simulation for the finite-$\beta$ configuration optimised for both QI and electron root.

Figure 31. Radial electric field profile obtained with the NTSS code by performing a power balance simulation for the finite-$\beta$ configuration optimised for both QI and electron root.

Figure 32. Neoclassical (a) and total (b) energy fluxes obtained for the power balance simulation performed for the finite-$\beta$ magnetic configuration optimised for electron root and QI. The simulation was performed with the NTSS code and the energy fluxes are normalised to the volume-averaged alpha power of $P_{\alpha }\approx 314\ \textrm {MW}$.

5.8 Impurity transport properties

In this section, we discuss briefly the properties of impurities. Two types of impurities are of concern for future stellarator reactors, helium ash and heavy tungsten impurities. We compare the properties of these two species for the configuration optimised only for QI and the vacuum and finite-$\beta$ configurations optimised for both QI and the electron root. To analyse the helium ash, we follow Beidler et al. (Reference Beidler, Drevlak, Geiger, Helander, Smith and Turkin2024) and use the ratio $S_{\alpha }\tau _\textrm {He}/n_\textrm {He}$ to quantify if the helium ash is transported effectively outwards. Here the source of helium ash $S_{\alpha }$ is taken to be the source of fusion alpha particles and the helium ash confinement time is calculated as $\tau _\textrm {He}=\int _0^1 \textrm {d}\rho \,n_\textrm {He}\rho / \int _0^1 \textrm {d}\rho \,S_{\alpha } \rho$. Such quantities can be obtained from the power balance simulations performed in the previous section. As mentioned in Beidler et al. (Reference Beidler, Drevlak, Geiger, Helander, Smith and Turkin2024), we expect that, if the quantity $S_{\alpha }\tau _\textrm {He}/n_\textrm {He}$ is well above one, the helium ash will be transported more effectively outwards. The quantities $S_{\alpha }$ and $S_{\alpha }\tau _\textrm {He}/n_\textrm {He}$ are shown in figure 33 for the three configurations of interest. We can see that in the electron root region, $S_{\alpha }\tau _\textrm {He}/n_\textrm {He}$ is considerably larger for both the vacuum and finite-$\beta$ configurations optimised for both QI and the electron root when compared with the configuration optimised only for QI. We can see in figure 33, that in the configurations optimised for both QI and the electron root, the quantity $S_{\alpha }\tau _\textrm {He}/n_\textrm {He}$ can reach values greater than two. Such an increase of $S_{\alpha }\tau _\textrm {He}/n_\textrm {He}$ in those configurations, happens both in zones, which show a value $S_{\alpha }$ greater and lower than in the configuration optimised only for QI. Such observations indicate that the electron root region obtained for the QI and the electron root optimised configurations is beneficial with respect to the exhaust of helium ash.

Figure 33. Source of fusion $\alpha$ particles (a) and the ratio $S_{\alpha }\tau _\textrm {He}/n_\textrm {He}$ (b) for the configuration optimised only for QI, and the vacuum and finite-$\beta$ configurations optimised for both QI and electron root. Such quantities are obtained from the power balance simulations performed with NTSS.

We now look into the properties of heavy tungsten particles, which we take here to be in a single ionisation state with charge $Z_W=40$ and a mass number of $184.840$. We also consider $T_W=T_i=(T_D+T_T)/2$ so that the tungsten temperature is such that $T_W\sim T_i$. To keep this treatment simple we assume that the tungsten is a trace species and its density is such that $n_W=10^{-5}n_i$. This assumption allows $Z_W^2n_W/n_i\sim 10^{-2}\ll 1$ so that tungsten will not have an influence in the background species either through modification of the ambipolarity constraint or by collisions. This consideration is in line with the ‘safe’ concentration limit of tungsten seen in tokamak experiments (Pütterich et al. Reference Pütterich, Dux, Neu, Bernert, Beurskens, Bobkov, Brezinsek, Challis, Coenen and Coffey2013) and it should serve here to elucidate about tungsten physics before any accumulation or peaked density profiles are observed. Such an assumption also allows the usage of the profiles obtained from the previous power balance simulations. We use such profiles in SFINCS to solve the drift kinetic equation in the form (2.11), considering four species, i.e. tungsten, deuterium, tritium and electrons, with a full linearised Fokker–Planck collision operator. Such an approach to solve the drift kinetic equations is necessary here because tungsten is expected to be highly collisional, due to its high charge, and its particle-flux-density is expected to be affected by collisions with the main thermal ions. The tungsten particle-flux-density for the QI-only optimised configuration and the vacuum and finite-$\beta$ configurations optimised for both QI and electron root can be seen in figure 34. We observe that the tungsten particle-flux-density for the configuration optimised only for QI has a region of positive values despite having an ion root (see figure 25). This is due to a combination of the screening provided by the temperature gradient and the fact that the values of the electric field, in such region of the ion root solution of this configuration, are small. However, the tungsten particle-flux-density for this configuration is overall around one order of magnitude smaller than the particle-flux-densities for the other configurations. We see that for both the vacuum and finite-$\beta$ configurations optimised for QI and the electron root, there is a large value of the tungsten particle-flux-density in the electron root region (see figures 28 and 31). Therefore, as expected such configurations show stronger beneficial effects on the removal of tungsten particles from the plasma core. We point out that such results are similar to the ones observed in Fujita et al. (Reference Fujita, Satake, Nunami, García-Regaña, Velasco and Calvo2021) and Velasco et al. (Reference Velasco, Calvo, Satake, Alonso, Nunami, Yokoyama, Sato, Estrada, Fontdecaba and Liniers2016), which model impurity hole plasmas in LHD scenarios, but here, such behaviour is found at reactor-relevant plasma and configuration parameters. We end this section by mentioning that, though small, the positive values of the QI configuration could provide for this plasma equilibrium an outward flux of tungsten. Nevertheless, such an effect is weaker than for the other configurations in the electron root region and is dependent on the temperature profiles achieved. The large tungsten particle-flux-densities observed in the electron root region for the configurations optimised for electron root indicate that in such conditions the beneficial outward impurity flux should be present for most temperature profiles.

Figure 34. Neoclassical particle-flux-density of tungsten obtained from a SFINCS simulation with tungsten, deuterium, tritium and electrons. The particle-flux-density is shown for the configuration optimised only for QI and both the vacuum and finite-$\beta$ configurations are optimised for QI and electron root. A full Fokker–Planck collision operator is used for the four species and the profiles obtained in the power balance simulations are used.

5.9 Coil configurations

For completeness, we show two sets of coils obtained for the QI and electron root optimised configurations. One set of coils is obtained for the vacuum configuration and another for the finite-$\beta$ configuration. Both of these were obtained with the SIMSOPT code using a second-stage optimisation approach. In the finite-$\beta$ configuration, a virtual casing principle is used to separate the magnetic field due to the coils from the plasma contribution. The coils obtained for the vacuum configuration case together with the plasma boundary are presented in figure 35 for half of a field period. The surface contour plot of the boundary represents the residual normal magnetic field at the boundary of the plasma (normalised to the magnetic field) which should vanish when the coil field matches the target field from the optimisation. The set of coils is composed of four independent sets of 10 coils, each one representing half of a field period, due to the constraints given by assuming stellarator symmetry and $n_\textrm {fp}=2$. They have a minimum separation between them of 0.1 m, and each has a total length of 38.5, 36.9, 28.2, 26.7, 28.2, 31.4, 34.3, 38.2, 37.2 and 40.2 m, respectively. They have a maximum curvature of $1.9\ \textrm {m}^{-1}$ and a maximum mean squared curvature of $0.6\ \textrm {m}^{-1}$. One of the coils has a fixed current of $10^6$ A, and the remaining coils have a maximum current of $1.31 \times 10^6$ A and a minimum current of $5.6 \times 10^5$ A. The poloidal cross-section, at toroidal angle $\phi =0$, of the respective Poincaré plots is shown in figure 36 that demonstrates that such coils can accurately reproduce the optimised fixed boundary equilibrium.

Figure 35. Optimised coils obtained for half of a field period for the vacuum configuration optimised for both QI and electron root. A surface plot of the plasma boundary representing the normal magnetic field at the boundary normalised to the magnetic field is also shown.

Figure 36. Poloidal cross-section at $\phi =0$ of the three-dimensional Poincaré plots for the optimised set of coils obtained for the vacuum configuration optimised for both QI and electron root.

The set of optimised coils for half of a field period obtained for the finite-$\beta$ configuration optimised for QI and electron root can be seen in figure 37. This set of coils is again one of the four equal sets of 10 different coils that compose the total array of coils for the entire plasma. They have a minimum separation between them of 0.26 m, and each has a total length of 25.6, 26.7, 27.2, 28.1, 32.2, 37.3, 39.1, 41.5, 39.5 and 38.2 m, respectively. They have a maximum curvature of $1.7\ \textrm {m}^{-1}$ and a maximum mean squared curvature of $0.5\ \textrm {m}^{-1}$. The coils have a fixed total sum of currents of $1.42 \times 10^8$ A with each current of approximately 10 % of the total current. The surface contour plot of the magnetic field normal to the plasma boundary normalised to the magnetic field is also shown in figure 37 for half of a field period. The normal magnetic field shown is comprised of the coil magnetic field computed using the Biot–Savart law minus the target magnetic field obtained using the virtual casing principle. No Poincaré plots are presented for the finite-$\beta$ configuration because they would only represent the vacuum contribution of the magnetic field. To consider both the vacuum and the plasma contributions a finite-$\beta$ optimisation considering the free boundary solution of the MHD force balance is required. However, such a simulation is out of the scope of this paper.

Figure 37. Optimised coils obtained for the finite-$\beta$ for half of the field period of the configuration optimised for both QI and electron root. A surface plot of the plasma boundary representing the normal magnetic field at the boundary normalised to the magnetic field is also shown.

5.10 Intuitive picture of electron root optimisation

To obtain an intuitive picture of what is being achieved by the electron root cost function during optimisation, we may look into the behaviour of the monoenergetic radial transport scans of the configurations which are optimised with the electron root cost function and into the cost function itself. There are two ways the ratio $L_{11}^i/L_{11}^e$ can be minimised. One is by increasing the radial transport coefficient of the electrons and another one is by decreasing the ion radial transport coefficient at the values of collisionality of interest. Ideally, we expect the optimisation to achieve an electron root by targeting both ways of minimising the cost function. From the monoenergetic radial transport coefficient scans in figures 12, 15 and 21, we find that two trends appear. The first is the increasing of the overall magnitude of the $1/\nu ^*$ regime and the second is the existence of two well-defined regimes at low collisionality and finite electric field. These are the $1/\nu ^*$ and $\nu ^*$ regimes, as a consequence of the disappearance of the $\sqrt {\nu ^*}$ regime. The electron transport increase is related to the increase in the $1/\nu ^*$ regime, in which this species usually lives due to its large thermal velocity. Therefore, the disappearance of the $\sqrt {\nu ^*}$ must be related to the ion species’ role in minimising the electron root cost function. To better understand these trends we can look at how the distribution function solution is obtained for these regimes. Following Shaing et al. (Reference Shaing, Ida and Sabbagh2015), the contribution to the non-ambipolar radial transport is completely described by the solution of the bounce-averaged drift kinetic equation. By using the variables $(\psi,\alpha,\phi )$, and neglecting the mirror term (which is small at large aspect ratios), the local bounce-averaged drift kinetic equation can be written as

(5.1)\begin{equation} \langle \boldsymbol{v_d}\boldsymbol{\cdot}\boldsymbol{\nabla} \alpha\rangle_B\frac{\partial f_{1}}{\partial \alpha}+\langle \boldsymbol{v_d}\boldsymbol{\cdot}\boldsymbol{\nabla} \psi\rangle_B\frac{\partial f_{M}}{\partial \psi}=\langle C(f_{1})\rangle_B, \end{equation}

where the bounce average of a quantity $X$ is defined as

(5.2)\begin{equation} \langle X\rangle_B =\int^{\phi_B}_{-\phi_B}\left.\frac{B X \,{\rm d} \phi}{v_{{\parallel}}}\right/\int^{\phi_B}_{-\phi_B}B\,{\rm d} \phi/v_{{\parallel}}. \end{equation}

For large aspect ratio stellarators, the most important drift across field lines is the poloidal $E\times B$ drift. In the $1/\nu ^*$ regime, the poloidal $E\times B$ drift frequency $\omega _{E\times B}$ is smaller than the collision frequency, and the bounce-averaged drift kinetic equation to be solved (Shaing et al. Reference Shaing, Ida and Sabbagh2015; Calvo et al. Reference Calvo, Parra, Velasco and Alonso2017) is

(5.3)\begin{equation} \langle \boldsymbol{v_d}\boldsymbol{\cdot}\boldsymbol{\nabla} \psi\rangle_B\frac{\partial f_{M}}{\partial \psi}=\langle C(f_{1})\rangle_B, \end{equation}

which represents a balance between the bounce-averaged radial drift and the collisions. The collisions limit the radial excursions of the helical trapped particles. Thus, if the collision frequency increases, the radial excursion decreases and we have a $1/\nu ^*$ scaling for the radial transport coefficient. The deeply trapped electrons with their large thermal velocity are usually well represented by this regime. At lower collisionalities, the collisions get weaker and the $E\times B$ drift becomes important in reducing the radial excursion resulting from the radial drift. An expansion in the small parameter $\nu _{aa}/\omega _{E\times B}$ is usually applied (Shaing et al. Reference Shaing, Ida and Sabbagh2015; Calvo et al. Reference Calvo, Parra, Velasco and Alonso2017). The equations to solve are then

(5.4)\begin{gather} \langle \boldsymbol{v_d}\boldsymbol{\cdot}\boldsymbol{\nabla} \alpha\rangle_B\frac{\partial f_{1}^0}{\partial \alpha}+\langle \boldsymbol{v_d}\boldsymbol{\cdot}\boldsymbol{\nabla} \psi\rangle_B\frac{\partial f_{M}}{\partial \psi}=0, \end{gather}

(5.5)\begin{gather}\langle \boldsymbol{v_d}\boldsymbol{\cdot}\boldsymbol{\nabla}\alpha\rangle_B\frac{\partial f_{1}^1}{\partial \alpha}=\langle C(f_{1}^{0})\rangle_B, \end{gather}

in which the superscript indicates the ordering in $\nu _{aa}/\omega _{E\times B}$. In (5.4) and (5.5), the balance between the poloidal and radial drifts is the main mechanism in defining the distribution function. Ions feel this effect before the electrons because their thermal velocity is smaller. Therefore, their radial excursion is smaller and easier to limit by the poloidal drift. The solution to (5.4) and (5.5) is discontinuous at the trapping–passing boundary which leads to unphysical diverging fluxes (Shaing et al. Reference Shaing, Ida and Sabbagh2015; Calvo et al. Reference Calvo, Parra, Velasco and Alonso2017) if the region closed to that boundary is considered. Thus, the solution to (5.4) and (5.5) is only well-defined for particles which are not close to the trapping–passing boundary. For these particles, the radial transport is related to the different total drift they feel when moving in the magnetic wells. These different values of the drifts allow trapped particles to be detrapped and then trapped again without collisions. Such a regime leads to a radial transport that scales linearly with $\nu ^*$ and decreases with the electric field (Shaing et al. Reference Shaing, Ida and Sabbagh2015; Calvo et al. Reference Calvo, Parra, Velasco and Alonso2017). This regime is not to be confused with the banana-like regime in which the radial transport also scales with $\nu ^*$, but is independent of the electric field. The banana regime also exists in stellarators but does not contribute to the solution of the bounce-averaged drift kinetic equation, because it generates an intrinsically ambipolar radial transport. The discontinuity of the solution to (5.4) and (5.5) at the passing–trapped boundary is resolved by adding the collision term effect at the trapped/passing boundary (Shaing et al. Reference Shaing, Ida and Sabbagh2015; Calvo et al. Reference Calvo, Parra, Velasco and Alonso2017). The behaviour of the barely trapped particles is only completely described when (5.1) is considered. This means that these particles are so barely trapped that even a very weak collision frequency is enough to detrap them. This resonance layer physics is well known and leads to the $\sqrt {\nu ^*}$ regime (Shaing et al. Reference Shaing, Ida and Sabbagh2015; Calvo et al. Reference Calvo, Parra, Velasco and Alonso2017). Because it has a resonant behaviour, the part of the distribution function that represents the solution of the barely trapped particles is usually the most important contribution to the total solution of the distribution, unless the magnetic configuration is such that the bounce-averaged radial drift is kept finite and small near $B_\textrm {max}$, leading to a small contribution of these particles to the bounce-averaged component of the distribution function. Thus, the vanishing of the $\sqrt {\nu ^*}$ regime observed in figures 12, 15 and 21 indicates that the electron root cost function is minimised by keeping the bounce-averaged radial drift contribution for the barely trapped ions small.

Therefore, the cost function for the electron root is a balance of minimising the radial transport of barely trapped ions and maximising the radial transport of the deeply trapped electrons. As seen in Rodríguez & Plunk (Reference Rodríguez and Plunk2023), the bounce-averaged drift can be measured by the magnetic quantity

(5.6)\begin{equation} \varDelta_{Y}=Y(\phi_B)-Y(-\phi_B), \quad Y=\frac{\boldsymbol{\nabla} \psi \times B \boldsymbol{\cdot}\boldsymbol{\nabla} B}{B\boldsymbol{\cdot}\boldsymbol{\nabla}B}. \end{equation}

We can obtain, for the different magnetic configurations presented in this work, approximated measures of such quantity for the barely trapped ions and the deeply trapped electrons, and see if a relation is seen between them and the electron root optimisation. To do so, we measure the bounce-averaged radial drift using the magnetic quantity in (5.6). We calculate such quantity at $r/a=0.5$ for the consecutive maxima, $\varDelta _{Y_\textrm {max}}$ or minima, $\varDelta _{Y_\textrm {min}}$, present in the magnetic field spectrum of the configuration between $0$ and $10{\rm \pi}$. We then perform an average over the absolute value of such quantity for all the pairs of maxima or minima, respectively. Because different field lines can have distinct magnetic wells, we also perform an average over different field lines, in which we consider $51$ values of the field line label $\alpha$ uniformly distributed between $0$ and $2{\rm \pi}$. The final result is shown for the different configurations in table 1. The ratio $\varDelta _{Y_\textrm {max}}/\varDelta _{Y_\textrm {min}}$, which represents the relative strength of the bounce-averaged drift of barely trapped particles to the bounce-averaged radial drift of the deeply trapped particles, is also shown in table 1. We see that configurations that do not have an electron root for the density and temperature profiles used during optimisation, i.e. the initial configuration and the configuration only optimised for QI (see figures 7 and 10), have a ratio $\varDelta _{Y_\textrm {max}}/\varDelta _{Y_\textrm {min}}$ close or larger than $1$. On the other hand, configurations that are optimised for an electron root, have a value of the ratio $\varDelta _{Y_\textrm {max}}/\varDelta _{Y_\textrm {min}}$ smaller than $0.5$. In fact, the vacuum configuration optimised for both QI and electron root has the lower value of $\varDelta _{Y_\textrm {max}}/\varDelta _{Y_\textrm {min}}$ and also shows the largest electron root for the density and temperature profiles used during optimisation. Therefore, it is seen that the ratio $\varDelta _{Y_\textrm {max}}/\varDelta _{Y_\textrm {min}}$ has similar behaviour to the electron root cost function, which then justifies the intuitive picture of the electron root optimisation presented. We also note that these observations agree with the work in Helander et al. (Reference Helander, Goodman, Beidler, Kuczynski and Smith2024), which tailors the magnetic field during optimisation to favour the presence of an electron root.

Table 1. Average of the quantity in (5.6) between consecutive maxima ($\varDelta _{Y_\textrm {max}}$) and minima ($\varDelta _{Y_\textrm {min}}$) over $51$ field lines at $r/a=0.5$. The ratio $\varDelta _{Y_\textrm {max}}/\varDelta _{Y_\textrm {min}}$ is also shown. Such quantities are presented for the different optimised magnetic configurations obtained in this work.

6 Conclusion

In this work, we propose a new optimisation method to obtain magnetic configurations that favour an electron root solution for the radial electric field. We motivate the cost function and the optimisation method using neoclassical transport theory and describe a method for testing the efficacy of the cost function. Using both methods, we show that the electron root optimisation leads to magnetic configurations which favour an electron root in the plasma. We obtain a vacuum magnetic configuration optimised for both QI and an electron root. This magnetic configuration has an electron root solution at a wide range of reactor-relevant density and temperature values. We show that it is possible to achieve an electron root in such configurations for strong values of magnetic field strength on-axis, such as $B_0=8 \ \textrm {T}$ and above. In fact, it is seen that such a configuration is robust to the presence of an electron root solution of the electric field for several values of $n_0$, $T_0$, $R_0$ and $B_0$, without the need to rely on tweaking the shape of the profiles. Such electric field solutions are obtained for the reactor-relevant case of $T_e=T_i$ and high densities, contrary to the case of current experimental devices, like W7-X and LHD, where such electron roots are obtained at low densities or by decoupling electron from ion temperatures and reaching very peaked electron temperatures with the use of ECRH. A finite-$\beta$ configuration optimised for both QI and electron root was also obtained. It was shown that such a configuration has an electron root at the temperature and density profiles consistent with the ones targeted during optimisation, showing that such solutions can also be pursued when considering finite-$\beta$ effects.

By performing power balance simulations we conclude that the vacuum and finite-$\beta$ configurations optimised for both QI and an electron root, show electron root solutions for a volume-averaged alpha power target of 300 MW. Such solutions indicate that the maximum value of neoclassical energy fluxes in the electron root region can be similar to the maximum values of neoclassical energy fluxes observed in the configuration optimised only for QI and even lower than the maximum of the neoclassical energy fluxes observed in Beidler et al. (Reference Beidler, Smith, Alonso, Andreeva, Baldzuhn, Beurskens, Borchardt, Bozhenkov, Brunner and Damm2021) for W7-X configurations. Such a result is observed, despite the fact that both W7-X and the configuration optimised only for QI have much lower $\varepsilon _\textrm {eff}$ values than the ones optimised for both QI and an electron root. Such observation indicates that if the configurations are optimised for an electron root to be present in a large extension of the plasma radius, the $\varepsilon _\textrm {eff}$ may not be necessary as an optimisation metric, which may allow for a less restrictive parameter space during multiobjective optimisations.

Concerning impurities, evidence of beneficial effects for the exhaust of helium ash and the reduction of core tungsten accumulation is observed for the solutions of these power balance simulations performed for the vacuum and finite-$\beta$ configurations optimised for both QI and an electron root. Additionally, the large electric field gradients observed in the root transitions of these solutions are expected to generate large $E\times B$ shear flows that may play a role in reducing turbulence. As side note, we remark here that if such a turbulence reduction is verified in large root transitions, it would mean that an ion and electron temperature decoupling would cease to exist, and thus, an electron root obtained via this decoupling would also cease to exist. Such argument further motivates the optimisations done in this work, in which configurations with robust electron root solutions at $T_e=T_i$ are achieved. We point out that varying the parameters $n_0$, $T_0$, $a_0$ and $B_0$ was seen to have an impact on the electron root maximum and root transition position. Thus, further tuning of such parameter space can be used to enhance the beneficial effects of the electron root on impurities and (possibly) turbulence. Both the impurity behaviour and the possible turbulence reduction can greatly affect plasma performance, thus making electron root solutions worth studying as an operation scenario for future stellarator fusion reactors.

Future paths for this work include doing multiobjective optimisations, by adding other cost functions to (3.5) and (3.8). Targets such as MHD stability, fast particle confinement, bootstrap current or turbulence could be considered. Such work would ideally be performed together with the analysis of the different trade-offs between the electron root cost function and the other cost functions. Another important subject is the study of the impact of these new optimised electron root configurations on heavy impurities such as tungsten. Though here we have shown the beneficial effects of the electron root solutions obtained on the improvement of helium ash exhaust and the reduction of tungsten accumulation, further impurity simulations should be performed considering the effect of tungsten density asymmetries. Such simulations can be used to analyse the effect of different electron root parameters, such as the maximum and root transition location, on heavy impurity transport. Ideally, this would result in optimal points in plasma and configuration parameters space for helium ash exhaust and tungsten reduction. Finally, a study on how the electron root solutions obtained with these optimisations can impact turbulence should be performed. A study of the impact of different electron root parameters, such as root transition steepness, location and width should be studied to find optimal values of the plasma and configuration space for the possible reduction of turbulence. As a final remark, an important question concerns whether we can optimise directly for such parameters and if we can define the limits to obtain a reduction of turbulence. Regarding such a topic, we point out that the gradient length of the root transition should ideally stay within an ion gyroradius to have a strong effect in reducing turbulence.

Acknowledgments

We thank M. Landreman, H.M. Smith, A.G. Goodman, B.F. Lee and J.L. Velasco for helpful discussions and the SIMSOPT, DKES, SFINCS and NTSS teams for their valuable contributions.

Editor P. Catto thanks the referees for their advice in evaluating this article.

Funding

This work has been carried out within the framework of the EUROfusion Consortium, partially funded by the European Union via the Euratom Research and Training Programme (grant agreement no. 101052200 — EUROfusion). Contributions to this work also include funding from Proxima Fusion. Views and opinions expressed are, however, those of the author(s) only and do not necessarily reflect those of the European Union, the European Commission or Proxima Fusion. Neither the European Union nor the European Commission nor Proxima Fusion can be held responsible for them. R.J. is supported by the Portuguese FCT-Fundação para a Ciência e Tecnologia, under grant 2021.02213.CEECIND and DOI 10.54499/2021.02213.CEECIND/CP1651/CT0004. IPFN activities were supported by FCT – Fundação para a Ciência e Tecnologia, I.P. by project reference UIDB/50010/2020 and DOI 10.54499/UIDB/50010/2020, by project reference UIDP/50010/2020 and DOI 10.54499/UIDP/50010/2020 and by project reference LA/P/0061/202 and DOI 10.54499/LA/P/0061/2020. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under contract no. DE-AC02-05CH11231 using NERSC award NERSC DDR-ERCAP0030134. This work used Jetstream2 at Indiana University through allocation PHY240054 from the Advanced Cyberinfrastructure Coordination Ecosystem: Services & Support (ACCESS) program, which is supported by National Science Foundation grants #2138259, #2138286, #2138307, #2137603 and #2138296.

Declaration of interests

The authors report no conflict of interest.

Appendix A. The DKES–SFINCS normalisations

The DKES-like database used as an input in the NTSS code, assumes the DKES code to use as input parameters, the inverse mean free path $\nu /v$ and the normalised electric field $E_{\tilde {r}}/(vB_0)$ where typically one sets $B_0=1$. The coordinate $\tilde {r}$ in DKES is defined so that

(A1)\begin{equation} \frac{{\rm d} \psi}{{\rm d}\tilde{r}}=rB_0, \quad \frac{{\rm d} \tilde{r}}{{\rm d} r}=\frac{2 \psi_B}{a^2B_0}, \end{equation}

with the minor radius $a=\sqrt {V(\psi _B)/(2{\rm \pi} ^2R_0(\psi _B))}$. Using these definitions the SFINCS monoenergetic input variables $\nu ^{\prime }$ and $E_{\psi }^*$ can be written in terms of the DKES inputs as

(A2)\begin{equation} \nu^{\prime}=\frac{G+\iota I}{B_0}\left[\frac{3\sqrt{\rm \pi}}{4}\left(\left.\text{erf}(x=1)-\frac{\text{erf}(x=1)}{2x^2}+\frac{2}{2x}\frac{{\rm d}[\text{erf}(x)]}{{\rm d}x}\right|_{x=1}\right)\right]^{{-}1}\frac{\nu}{v} \end{equation}

and

(A3)\begin{equation} E_{\psi}^*={-}\frac{G}{\iota}\left(\frac{{\rm d}\psi}{{\rm d}\tilde{r}}\right)^{{-}1}\frac{E_{\tilde{r}}}{v}, \end{equation}

respectively. In (A2), $\textrm {erf}(x)=2/\sqrt {{\rm \pi} }\int _0^x \exp (-y^2) \,{\textrm {d}y}$ is the error function, and the terms containing such function come from the definition of the deflection frequency in Helander & Sigmar (Reference Helander and Sigmar2002). The radial transport coefficient as per the DKES-like database files is written in terms of the SFINCS monoenergetic radial coefficient $D_{11}^S$ as

(A4)\begin{equation} \hat{\varGamma}_{11}=\frac{\sqrt{\rm \pi}G^2 B_0}{8(G+\iota I)\left(\dfrac{{\rm d} \psi}{{\rm d}\tilde{r}}\right)^2}D_{11}^S. \end{equation}

In Beidler et al. (Reference Beidler, Allmaier, Isaev, Kasilov, Kernbichler, Leitold, Maaßberg, Mikkelsen, Murakami and Schmidt2011) the normalised monoenergetic radial coefficient is defined as $D_{11}^*=D_{11}/D_{{11}_p}$, with $D_{{11}_p}$ the value of $D_{11}$ in the plateau regime of an equivalent tokamak. This quantity can be defined in terms of $D_{11}^S$ as

(A5)\begin{equation} D_{11}^*={-}\frac{\iota R_0 G^2 B_0^3\left(\dfrac{{\rm d}\tilde{r}}{{\rm d}r}\right)^2}{\sqrt{\rm \pi} (G+\iota I)\left(\dfrac{{\rm d}\psi}{{\rm d}\tilde{r}}\right)^2}D_{11}^S. \end{equation}

References

Albert, C.G., Kasilov, S.V. & Kernbichler, W. 2020 Accelerated methods for direct computation of fusion alpha particle losses within, stellarator optimization. J. Plasma Phys. 86 (2), 815860201.CrossRef Google Scholar

Bader, A., Drevlak, M., Anderson, D.T., Faber, B.J., Hegna, C.C., Likin, K.M., Schmitt, J.C. & Talmadge, J.N. 2019 Stellarator equilibria with reactor relevant energetic particle losses. J. Plasma Phys. 85 (5), 905850508.CrossRef Google Scholar

Beidler, C.D., Drevlak, M., Geiger, J., Helander, P., Smith, H.M. & Turkin, Y. 2024 Reduction of neoclassical bulk-ion transport to avoid helium-ash retention in stellarator reactors. Nucl. Fusion 64 (12), 126030.CrossRef Google Scholar

Beidler, C.D., Allmaier, K., Isaev, M.Yu., Kasilov, S.V., Kernbichler, W., Leitold, G.O., Maaßberg, H., Mikkelsen, D.R., Murakami, S., Schmidt, M., et al. 2011 Benchmarking of the mono-energetic transport coefficients—results from the international collaboration on neoclassical transport in stellarators (icnts). Nucl. Fusion 51 (7), 076001.CrossRef Google Scholar

Beidler, C.D., Smith, H.M., Alonso, A., Andreeva, T., Baldzuhn, J., Beurskens, M.N.A., Borchardt, M., Bozhenkov, S.A., Brunner, K.J., Damm, H., et al. 2021 Demonstration of reduced neoclassical energy transport in Wendelstein 7-X. Nature 596, 221–226.CrossRef Google Scholar PubMed

Beurskens, M.N.A., Bozhenkov, S.A., Ford, O., Xanthopoulos, P., Zocco, A., Turkin, Y., Alonso, A., Beidler, C., Calvo, I., Carralero, D., et al. 2021 Ion temperature clamping in Wendelstein 7-X electron cyclotron heated plasmas. Nucl. Fusion 61 (11), 116072.CrossRef Google Scholar

Boozer, A.H. 1981 Plasma equilibrium with rational magnetic surfaces. Phys. Fluids 24 (11), 1999.CrossRef Google Scholar

Buller, S., Smith, H.M., Helander, P., Mollén, A., Newton, S.L. & Pusztai, I. 2018 Collisional transport of impurities with flux-surface varying density in stellarators. J. Plasma Phys. 84 (4), 905840409.CrossRef Google Scholar

Buller, S., Smith, H.M. & Mollén, A. 2021 Recent progress on neoclassical impurity transport in stellarators with implications for a stellarator reactor. Plasma Phys. Control. Fusion 63 (5), 054003.CrossRef Google Scholar

Calvo, I., Parra, F., Velasco, J., Alonso, J. & García-Regaña, J.M. 2018 Stellarator impurity flux driven by electric fields tangent to magnetic surfaces. Nucl. Fusion 58 (12), 124005.CrossRef Google Scholar

Calvo, I., Parra, F.I., Velasco, J.L. & Alonso, J.A. 2017 The effect of tangential drifts on neoclassical transport in stellarators close to omnigeneity. Plasma Phys. Control. Fusion 59 (5), 055014.CrossRef Google Scholar

Camacho Mata, K., Plunk, G.G. & Jorge, R. 2022 Direct construction of stellarator-symmetric quasi-isodynamic magnetic configurations. J. Plasma Phys. 88 (5), 905880503.CrossRef Google Scholar

Cary, J.R. & Shasharina, S.G. 1997 Omnigenity and quasihelicity in helical plasma confinement systems. Phys. Plasmas 4 (9), 3323–3333.CrossRef Google Scholar

Dewar, R.L. & Hudson, S.R. 1998 Stellarator symmetry. Physica D 112 (1–2), 275.CrossRef Google Scholar

Freidberg, J. 2014 Ideal MHD. Cambridge University Press.CrossRef Google Scholar

Fu, J.Y., Nicolau, J.H., Liu, P.F., Wei, X.S., Xiao, Y. & Lin, Z. 2021 Global gyrokinetic simulation of neoclassical ambipolar electric field and its effects on microturbulence in W7-X stellarator. Phys. Plasmas 28 (6), 062309.CrossRef Google Scholar

Fujita, K., Satake, S., Nunami, M., García-Regaña, J.M., Velasco, J.L. & Calvo, I. 2021 Study on impurity hole plasmas by global neoclassical simulation. Nucl. Fusion 61 (8), 086025.CrossRef Google Scholar

Fujiwara, M., Kawahata, K., Ohyabu, N., Kaneko, O., Komori, A., Yamada, H., Ashikawa, N., Baylor, L.R., Combs, S.K., deVries, P.C., et al. 2001 Overview of LHD experiments. Nucl. Fusion 41 (10), 1355.CrossRef Google Scholar

Goodman, A.G., Camacho Mata, K., Henneberg, S.A., Jorge, R., Landreman, M., Plunk, G.G., Smith, H.M., Mackenbach, R.J.J., Beidler, C.D. & Helander, P. 2023 Constructing precisely quasi-isodynamic magnetic fields. J. Plasma Phys. 89 (5), 905890504.CrossRef Google Scholar

Hastings, D.E. 1984 Development of a differential equation for the ambipolar electric field in a bumpy torus in the low collision frequency limit. Phys. Fluids 27 (9), 2372–2374.CrossRef Google Scholar

Hastings, D.E., Houlberg, W.A. & Shaing, K.C. 1985 The ambipolar electric field in stellarators. Nucl. Fusion 25 (4), 445.CrossRef Google Scholar

Helander, P. 2014 Theory of plasma confinement in non-axisymmetric magnetic fields. Rep. Prog. Phys. 77 (8), 087001.CrossRef Google Scholar PubMed

Helander, P., Goodman, A.G., Beidler, C.D., Kuczynski, M. & Smith, H.M. 2024 Optimised stellarators with a positive radial electric field. J. Plasma Phys. 90 (6), 175900602.CrossRef Google Scholar

Helander, P., Newton, S.L., Mollén, A. & Smith, H.M. 2017 Impurity transport in a mixed-collisionality stellarator plasma. Phys. Rev. Lett. 118, 155002.CrossRef Google Scholar

Helander, P. & Sigmar, D.J. 2002 Collisional Transport in Magnetized Plasmas. Cambridge University Press.Google Scholar

Helander, P. & Simakov, A.N. 2008 Intrinsic ambipolarity and rotation in stellarators. Phys. Rev. Lett. 101, 145003.CrossRef Google Scholar PubMed

Hirshman, S.P. & Whitson, J.C. 1983 Steepest-descent moment method for three-dimensional magnetohydrodynamic equilibria. Phys. Fluids 26 (12), 3553.CrossRef Google Scholar

Jorge, R., Dorland, W., Kim, P., Landreman, M., Mandell, N.R., Merlo, G. & Qian, T. 2024 Direct microstability optimization of stellarator devices. Phys. Rev. E 110, 035201.CrossRef Google Scholar PubMed

Jorge, R., Goodman, A., Landreman, M., Rodrigues, J. & Wechsung, F. 2023 Single-stage stellarator optimization: combining coils with fixed boundary equilibria. Plasma Phys. Control. Fusion 65, 074003.CrossRef Google Scholar

Jorge, R., Plunk, G.G., Drevlak, M., Landreman, M., Lobsien, J.F., Camacho Mata, K. & Helander, P. 2022 A single-field-period quasi-isodynamic stellarator. J. Plasma Phys. 88 (5), 175880504.CrossRef Google Scholar

Kappel, J., Landreman, M. & Malhotra, D. 2023 The magnetic gradient scale length explains why certain plasmas require close external magnetic coils. Plasma Phys. Control. Fusion 66, 025018.CrossRef Google Scholar

Kim, P., Buller, S., Conlin, R., Dorland, W., Dudt, D.W., Gaur, R., Jorge, R., Kolemen, E., Landreman, M., Mandell, N.R., et al. 2024 Optimization of nonlinear turbulence in stellarators. J. Plasma Phys. 90 (2), 905900210.CrossRef Google Scholar

Kuczynski, M.D., Kleiber, R., Smith, H.M., Beidler, C.D., Borchardt, M., Geiger, J. & Helander, P. 2024 Self-consistent, global, neoclassical radial-electric-field calculations of electron-ion-root transitions in the W7-X stellarator. Nucl. Fusion 64 (4), 046023.CrossRef Google Scholar

Landreman, M., Buller, S. & Drevlak, M. 2022 Optimization of quasisymmetric stellarators with self-consistent bootstrap current and energetic particle confinement. Phys. Plasmas 29, 082501.CrossRef Google Scholar

Landreman, M., Medasani, B., Wechsung, F., Giuliani, A., Jorge, R. & Zhu, C. 2021 SIMSOPT: a flexible framework for stellarator optimization. J. Open Source Softw. 6 (65), 3525.CrossRef Google Scholar

Landreman, M. & Paul, E. 2022 Magnetic fields with precise quasisymmetry for plasma confinement. Phys. Rev. Lett. 128 (3), 035001.CrossRef Google Scholar PubMed

Landreman, M., Smith, H.M., Mollén, A. & Helander, P. 2014 Comparison of particle trajectories and collision operators for collisional transport in nonaxisymmetric plasmas. Phys. Plasmas 21 (4), 042503.CrossRef Google Scholar

Lascas Neto, E., Jorge, R., Beidler, C.D. & Lion, J. 2024 Dataset from: electron root optimisation for stellarator reactor designs. https://zenodo.org/doi/10.5281/zenodo.11302191.Google Scholar

Lee, B.F., Lazerson, S.A., Smith, H.M., Beidler, C.D. & Pablant, N.A. 2024 Direct optimization of neoclassical ion transport in stellarator reactors. Nucl. Fusion 64 (10), 106054.CrossRef Google Scholar

Maassberg, H., Brakel, R., Burhenn, R., Gasparino, U., Grigull, P., Kick, M., Kuhner, G., Ringler, H., Sardei, F., Stroth, U., et al. 1993 Transport in stellarators. Plasma Physics and Controlled Fusion 35 (SB), B319.CrossRef Google Scholar

Najmabadi, F., Raffray, A.R., Abdel-Khalik, S.I., Bromberg, L., Crosatti, L., El-Guebaly, L., Garabedian, P.R., Grossman, A.A., Henderson, D., Ibrahim, A., et al. 2008 The aries-cs compact stellarator fusion power plant. Fusion Sci. Technol. 54 (3), 655–672.CrossRef Google Scholar

Naujoks, D., Dhard, C., Feng, Y., Gao, Y., Stange, T., Buttenschön, B., Bozhenkov, S.A., Brezinsek, S., Brunner, K.J., Cseh, G., et al. 2023 Performance of tungsten plasma facing components in the stellarator experiment W7-X: recent results from the first OP2 campaign. Nucl. Mater. Energy 37, 101514.CrossRef Google Scholar

Nemov, V.V., Kasilov, S.V., Kernbichler, W. & Heyn, M.F. 1999 Evaluation of 1/v neoclassical transport in stellarators. Phys. Plasmas 6 (12), 4622.CrossRef Google Scholar

Pablant, N., Langenberg, A., Alonso, A., Beidler, C., Bitter, M., Bozhenkov, S., Burhenn, R., Beurskens, M., Delgado-Aparicio, L.F., Dinklage, A., et al. 2018 Core radial electric field and transport in wendelstein 7-X plasmas. Phys. Plasmas 25, 022508.CrossRef Google Scholar

Pedersen, T.S., König, R., Krychowiak, M., Jakubowski, M., Baldzuhn, J., Bozhenkov, S., Fuchert, G., Langenberg, A., Niemann, H., Zhang, D., et al. 2018 First results from divertor operation in Wendelstein 7-X. Plasma Phys. Control. Fusion 61 (1), 014035.CrossRef Google Scholar

Pütterich, T., Dux, R., Neu, R., Bernert, M., Beurskens, M.N.A., Bobkov, V., Brezinsek, S., Challis, C., Coenen, J.W., Coffey, I., et al. 2013 Observations on the W-transport in the core plasma of JET and ASDEX Upgrade. Plasma Phys. Control. Fusion 55 (12), 124036.CrossRef Google Scholar

Roberg-Clark, G.T., Plunk, G.G., Xanthopoulos, P., Nührenberg, C., Henneberg, S.A. & Smith, H.M. 2023 Critical gradient turbulence optimization toward a compact stellarator reactor concept. Phys. Rev. Res. 5, L032030.CrossRef Google Scholar

Rodríguez, E. & Plunk, G.G. 2023 Higher order theory of quasi-isodynamicity near the magnetic axis of stellarators. Phys. Plasmas 30 (6), 062507.CrossRef Google Scholar

Sanchez, R., Hirshman, S.P., Ware, A.S., Berry, L.A. & Spong, D.A. 2000 Ballooning stability optimization of low-aspect-ratio stellarators. Plasma Phys. Control. Fusion 42 (6), 641.CrossRef Google Scholar

Shaing, K.C. 1984 Stability of the radial electric field in a nonaxisymmetric torus. Phys. Fluids 27 (7), 1567–1569.CrossRef Google Scholar

Shaing, K.C., Ida, K. & Sabbagh, S.A. 2015 Neoclassical plasma viscosity and transport processes in non-axisymmetric tori. Nucl. Fusion 55 (12), 125001.CrossRef Google Scholar

Sugama, H., Okamoto, M., Horton, W. & Wakatani, M. 1996 Transport processes and entropy production in toroidal plasmas with gyrokinetic electromagnetic turbulence. Phys. Plasmas 3 (6), 2379–2394.CrossRef Google Scholar

Tamura, N., Sudo, S., Suzuki, C., Funaba, H., Nakamura, Y., Tanaka, K., Yoshinuma, M., Ida, K. & The LHD Experiment Group 2016 Mitigation of the tracer impurity accumulation by EC heating in the LHD. Plasma Phys. Control. Fusion 58 (11), 114003.CrossRef Google Scholar

Tokitani, M., Masuzaki, S. & Murase, T. 2019 Demonstration of suppression of dust generation and partial reduction of the hydrogen retention by tungsten coated graphite divertor tiles in lhd. Nucl. Mater. Energy 18, 23–28.CrossRef Google Scholar

Turkin, Y., Beidler, C.D., Maaberg, H., Murakami, S., Tribaldos, V. & Wakasa, A. 2011 Neoclassical transport simulations for stellarators. Phys. Plasmas 18 (2), 022505.CrossRef Google Scholar

van Rij, W.I. & Hirshman, S.P. 1989 Variational bounds for transport coefficients in three-dimensional toroidal plasmas. Phys. Fluids B 1 (3), 563–569.CrossRef Google Scholar

Velasco, J.L., Calvo, I., Satake, S., Alonso, A., Nunami, M., Yokoyama, M., Sato, M., Estrada, T., Fontdecaba, J.M., Liniers, M., et al. 2016 Moderation of neoclassical impurity accumulation in high temperature plasmas of helical devices. Nucl. Fusion 57 (1), 016016.CrossRef Google Scholar

Wechsung, F., Landreman, M., Giuliani, A., Cerfon, A. & Stadler, G. 2022 Precise stellarator quasi-symmetry can be achieved with electromagnetic coils. Proc. Natl Acad. Sci. USA 119 (13), e2202084119.CrossRef Google Scholar PubMed

Yamada, H., Harris, J.H., Dinklage, A., Ascasibar, E., Sano, F., Okamura, S., Talmadge, J., Stroth, U., Kus, A., Murakami, S., et al. 2005 Characterization of energy confinement in net-current free plasmas using the extended international stellarator database. Nucl. Fusion 45 (12), 1684.CrossRef Google Scholar

Yokoyama, M., Maaßberg, H., Beidler, C.D., Tribaldos, V., Ida, K., Estrada, T., Castejon, F., Fujisawa, A., Minami, T., Shimozuma, T., et al. 2007 Core electron-root confinement (CERC) in helical plasmas. Nucl. Fusion 47 (9), 1213.CrossRef Google Scholar

Figure 1. Comparison between SFINCS (crosses) and DKES (circles) monoenergetic radial transport coefficients (normalised to DKES conventions) at $r/a=0.12247$ and $r/a=0.5$ for the Hydra-Np04-20190108 configuration presented in Beidler et al. (2024). Values are shown for the electric fields away from resonant values. A SFINCS resolution of $N_{\theta }=39$, $N_{\phi }=59$ and $N_{\xi }=116$ is used.

Figure 2. Comparison between SFINCS (crosses) and DKES (circles) monoenergetic radial transport coefficients at $r/a=0.12247$ and $r/a=0.5$ for the Hydra-Np04-20190108 configuration presented in Beidler et al. (2024). Comparison is shown for large electric fields near resonant values. A SFINCS resolution of $N_{\theta }=39$, $N_{\phi }=59$ and $N_{\xi }=116$ is used.

Figure 3. Comparison between SFINCS (crosses) and DKES (circles) monoenergetic radial transport coefficients at $r/a=0.12247$ and $r/a=0.5$ for the Hydra-Np04-20190108 configuration presented in Beidler et al. (2024). Comparison is shown for large electric field values near resonance. A SFINCS resolution of $N_{\theta }=39$, $N_{\phi }=431$ and $N_{\xi }=351$ is used.

Figure 5. Contours of the magnetic field strength of the initial configuration at $s=0.1025$ in Boozer coordinates and the same quantity at the plasma boundary.

Figure 6. Scans of the monoenergetic radial transport coefficients at $r/a=0.12247$ and $r/a=0.5$ for the initial magnetic configuration.

Figure 9. Scans of the monoenergetic radial transport coefficients at $r/a=0.12247$ and $r/a=0.5$ for the magnetic configuration optimised only for the QI cost function.

Figure 12. Scans of the monoenergetic radial transport coefficients at $r/a=0.12247$ and $r/a=0.5$ for the magnetic configuration optimised only for the electron root cost function.

Figure 15. Scans of the monoenergetic radial transport coefficients at $r/a=0.12247$ and $r/a=0.5$ for the magnetic configuration optimised for both the electron root and QI.

Figure 21. Scans of the monoenergetic radial transport coefficients at $r/a=0.12247$ and $r/a=0.5$ for the magnetic configuration optimised for both electron root and QI with finite-$\beta$ effects.

Figure 23. Epsilon effective and loss fraction of fusion alpha particles for the vacuum standard W7-X configuration, for the configuration optimised only for QI and for the vacuum and finite-$\beta$ configurations optimised for both QI and electron root. For the loss fraction of alpha particles, the volume-averaged magnetic field and minor radius of the configurations were scaled to $5.7\ \textrm {T}$ and $1.7\ \textrm {m}$, respectively, to mimic the ARIES-CS scaling (Najmabadi et al.2008). The alpha particles are initialised at $s=0.25$.

Figure 24. Density (a) and temperature (b) profiles obtained with the NTSS code by performing a power balance simulation for the configuration optimised only for QI.

Figure 25. Radial electric field profile obtained from the NTSS code by performing a power balance simulation for the configuration optimised only for QI.

Figure 27. Density (a) and temperature (b) profiles obtained from a power balance simulation, performed by the NTSS code for the vacuum configuration optimised for both QI and electron root.

Figure 28. Radial electric field profile obtained from a power balance simulation, performed with the NTSS code for the vacuum configuration optimised for both QI and electron root.

Figure 30. Density (a) and temperature (b) profiles obtained with the NTSS code by performing a power balance simulation for the finite-$\beta$ configuration optimised for both QI and electron root.

Figure 31. Radial electric field profile obtained with the NTSS code by performing a power balance simulation for the finite-$\beta$ configuration optimised for both QI and electron root.

Figure 36. Poloidal cross-section at $\phi =0$ of the three-dimensional Poincaré plots for the optimised set of coils obtained for the vacuum configuration optimised for both QI and electron root.

Article contents

Electron root optimisation for stellarator reactor designs

Abstract

Keywords

1 Introduction

2 Theoretical background

3 Optimisation method

4 Validation method

5 Results

5.1 The initial configuration

5.2 The QI optimisation

5.3 Electron root optimisation

5.4 The QI and electron root optimisation

5.5 Finite-$\beta$ optimisation

5.6 Confinement properties

5.7 Power balance considerations

5.8 Impurity transport properties

5.9 Coil configurations

5.10 Intuitive picture of electron root optimisation

6 Conclusion

Acknowledgments

Funding

Declaration of interests

Appendix A. The DKES–SFINCS normalisations

References

Save article to Kindle

Save article to Dropbox

Save article to Google Drive

Reply to: Submit a response

Your details

You have entered the maximum number of contributors

Conflicting interests