2. Model description

The M3D-C1 code is a high-fidelity extended-MHD code (Jardin et al. Reference Jardin, Ferraro, Breslau and Chen2012a ). It employs a split-implicit time scheme, which allows for time steps that extend beyond the Alfvénic time scale, enabling stable simulations on the transport time scale (Jardin Reference Jardin2012). The code utilizes high-order finite elements with $C^1$ continuity, constructed on an axisymmetric mesh. Recent developments have enabled M3D-C1 to model the MHD evolution of stellarator plasmas by accommodating non-axisymmetric domain geometries (Zhou et al. Reference Zhou, Ferraro, Jardin and Strauss2021). The governing equations for both tokamak and stellarator simulations are fully detailed in Jardin et al. (Reference Jardin, Ferraro, Breslau and Chen2012a ) and Zhou et al. (Reference Zhou, Ferraro, Jardin and Strauss2021), while the single-fluid model equations, relevant to the analysis of plasma behavior in this work, are reproduced in Appendix C.

In the present work, we extend M3D-C1 to include a non-inductive current source, specifically, the bootstrap current, by modifying Ohm’s law as follows:

(2.1) \begin{equation} \boldsymbol{E}=-\boldsymbol{v}{\times} \boldsymbol{B}+\eta [\boldsymbol{J}-\boldsymbol{J}_x]. \end{equation}

Where, $\boldsymbol{E}$ is the electric field, $\boldsymbol{v}$ is fluid velocity, $\eta$ is the electrical resistivity, and $\boldsymbol{J}$ is the current density. This is equivalent to adding a force on the electrons $\boldsymbol{F}_x^{e}$ such that $\boldsymbol{F}_x^{e}=-\eta \, n_e\, e\,\boldsymbol{J}_x$ , which corresponds to removing the contribution of $\boldsymbol{J}_x$ to the friction force on the electrons due to collisions with ions (Braginskii Reference Braginskii1965). Thus, collisional drag will cause the current density to decay towards $\boldsymbol{J}_x$ . Even if the force is applied instantly, the actual current in the plasma $\boldsymbol{J} = \boldsymbol{\nabla } {\times} \boldsymbol{B} / \mu _0$ will evolve on resistive time scales. How this force appears in the ion momentum equation depends on whether this is an internal force (e.g. the bootstrap current) or an external force (here, $\mu_0$ is the permeability of free space). For an internal force, there must be an equal and opposite force in the ion momentum equation, which results in no net force added to the MHD force balance equation.

For a force consistent with a purely parallel, divergence-free current ( $\boldsymbol{J}_x=J_{\parallel } ({\boldsymbol{B}}/{B})$ ), the quantity $J_{\parallel }/B$ must remain constant on a magnetic surface. Using this condition, we can relate the local parallel current to the magnetic field via the expression

(2.2) \begin{equation} J_{\parallel }=\frac {\langle \boldsymbol{J}_x \boldsymbol{\cdot }\boldsymbol{B}\rangle }{\langle B^2\rangle }{B}, \end{equation}

where $\langle \boldsymbol{\cdot }\rangle$ denotes the magnetic surface average. To calculate $\langle \boldsymbol{J_x \boldsymbol{\cdot }B} \rangle$ , the generalized neoclassical models of Sauter et al. (Reference Sauter, Angioni and Lin-Liu1999) and Redl et al. (Reference Redl2021) are employed, which provide analytical expressions for this quantity. The specific forms of these models are given in Appendix A.

For stellarators, the bootstrap current is computed using the method introduced by Landreman et al. (Reference Landreman, Buller and Drevlak2022) (shown in (2.3)), which exploits the isomorphism between axisymmetric and quasi-symmetric geometries

(2.3) \begin{align} \langle \boldsymbol{J_x \boldsymbol{\cdot }B}\rangle &= \frac {\tilde {G} }{\iota - N} \Bigg [ L_{31} \left ( p_e \frac {\partial \ln n_e}{\partial \psi _t} + p_i \frac {\partial \ln n_i}{\partial \psi _t} \right ) \notag \\ &\quad + p_e (L_{31} + L_{32}) \frac {\partial \ln T_e}{\partial \psi _t} + p_i (L_{31} + \alpha L_{34}) \frac {\partial \ln T_i}{\partial \psi _t} \Bigg ]. \end{align}

Here, $\tilde {G} (\psi _{t})=G+NI$ , $G$ is $\mu _0/2\pi$ times the poloidal current outside the flux-surface $(\psi _{t})$ , $I$ is $\mu _0/2\pi$ the toroidal current inside the flux-surface $(\psi _{t})$ , $\psi _{t}$ is toroidal flux per radian i.e. $\varPsi _t = 2\pi \,\psi _t$ , $\iota$ is the rotational transform, $p_{e/i}, n_{e/i}, T_{e/i}$ are the electron (ion) pressures, densities and temperatures and $\alpha$ , $L_{31}$ , $L_{32}$ , $L_{34}$ are the bootstrap coefficients (see Appendix B for further definitions). The Sauter–Redl–Landreman formulation requires information about the global equilibrium, quantities such as $G,I$ and $ \psi _t$ . This information is generally not known within M3D-C1, as its treatment of magnetic field does not assume the presence of magnetic surfaces. To address this challenge, a separate calculation is performed in which an approximate magnetic coordinate system is constructed from M3D-C1 output using Fusion-IO (Ferraro Reference Ferraro2025), taking isotherms of $T_e$ as proxies for magnetic surfaces. Due to the strongly anisotropic conduction of $T_e$ , these isotherms closely coincide with magnetic surfaces when surfaces exist and the dynamics is sufficiently slow; when surfaces do not exist, the isotherms are related to quadratic flux minimizing surfaces (Dewar, Hudson & Price Reference Dewar, Hudson and Price1994). However, it is worth reiterating that the Sauter–Redl–Landreman formula is valid only in quasi-symmetric geometries. Limitations introduced by deviations from quasi-symmetry are discussed in § 3.2.

To evaluate the neoclassical bootstrap current density and its related coefficients correctly, it is necessary to compute quantities such as the trapped particle fraction $f_t$ and geometric factors like the inverse aspect ratio $\epsilon$ and $qR$ (as defined in Appendix B). These parameters are integral to the expressions for the bootstrap current coefficients $\alpha$ , $L_{31}$ , $L_{32}$ and $L_{34}$ (see Redl et al. Reference Redl2021 for detailed expressions). To facilitate these calculations, we define a coordinate system based on the electron temperature isotherms, $\hat {T}_e=1-T_e / T_e^{max}$ with $T_e^{max}$ being the maximum electron temperature within the plasma domain. This is roughly analogous to the normalized flux $\varPsi _{t_N}=\psi _t/\psi _{LCFS}$ . Within this framework, the global equilibrium quantities $\langle I(\hat {T_e}) \rangle$ , $\langle G(\hat {T_e})\rangle$ , $\langle \hat {T_e} \rangle$ , along with $\langle f_t(\hat {T_e})\rangle$ , $\langle \epsilon (\hat {T_e})\rangle$ , $\langle qR(\hat {T_e})\rangle$ and $\iota (\hat {T_e}) = ({\rm d}\varPsi _p/{\rm d}\hat {T_e}) / ({\rm d}\varPsi _t/{\rm d}\hat {T_e})$ are calculated externally to M3D-C1 as flux-surface averages (all flux-surface averaging is performed externally). They are then read into M3D-C1 at the start of the simulation and extrapolated to local mesh coordinates. During the dynamical simulations, the bootstrap coefficients $\alpha$ , $L_{31}$ , $L_{32}$ and $L_{34}$ are evaluated locally at each time step using the evolving profiles of $\hat {T}_e$ . The global equilibrium quantities are re-extrapolated at the beginning of each time step to reflect the updated $\hat {T_e}$ distribution. Using $\hat {T}_e$ rather than $T_e$ enables a consistent treatment of cases where the temperature profile evolves in amplitude but maintains its shape, thereby avoiding unnecessary recomputation of the global quantities and other variables required for the calculation of the bootstrap coefficients.

To validate the accuracy of the bootstrap current calculations in M3D-C1 simulations, the results are compared with predictions from well-established neoclassical codes, namely: (i) XGCa, a global total-f gyrokinetic neoclassical code, where f is the 5D phase space distribution function (see Hager & Chang Reference Hager and Chang2016 for details), (ii) NEO, a drift-kinetic neoclassical steady-state solver, a detailed description of which can be found in Belli & Candy (Reference Belli and Candy2008, Reference Belli and Candy2011) and (iii) the stellarator neoclassical code SFINCS (the stellarator Fokker–Planck iterative neoclassical conservative solver), which solves the radially local four-dimensional drift-kinetic equation without assuming quasi-symmetry. See Landreman et al. (Reference Landreman, Smith, Mollén and Helander2014) for further details. These comparisons are performed for both tokamak (NEO, XGCa, SFINCS) and quasi-symmetric stellarator (SFINCS) geometries, ensuring consistent and reliable results across a range of magnetic confinement configurations.

3. Numerical results: bootstrap current calculations

3.1. Tokamak case verification

For the verification study, the low aspect ratio ‘CIRC1’ case from Hager & Chang (Reference Hager and Chang2016) is considered here. The configuration features a circular cross-section with inverse aspect ratio $\epsilon = 0.84$ at the outer boundary. In this analysis, following Hager & Chang (Reference Hager and Chang2016), the inverse aspect ratio $\epsilon$ is defined as the ratio of the mean minor radius $a_{\textit{mean}} = (R_{\textit{max}} - R_{\textit{min}})/2$ to the geometrical center $R_c= (R_{\textit{max}} + R_{\textit{min}})/2$ of a flux surface, where $R_{\textit{max}}$ and $R_{\textit{min}}$ are the maximum and minimum major radii on a flux surface. The density and temperature profiles used in the simulations are tanh-type pedestals, as shown in figure 1(a).

Figure 1(b) shows the bootstrap current profile for the CIRC1 configuration. Maximum bootstrap current occurs at $\psi _N \sim 0.89$ , where $\epsilon = 0.63$ and the electron collisionality $\nu _e^*$ is 0.92. A comparison between the M3D-C1 code using the Sauter et al. (Reference Sauter, Angioni and Lin-Liu1999) model and the NEO code’s Sauter model reveals a 2.18 % difference at the peak of the bootstrap current profile. Furthermore, the bootstrap current profiles from M3D-C1 and XGCa are nearly identical, with a difference of only 0.04 $\%$ at the peak. The M3D-C1 profile calculated using the Redl et al. (Reference Redl2021) model shows a 2.12 % difference from XGCa and a 1.01 $\%$ difference from SFINCS profiles at the peak. It is important to highlight that the electron collisionality in this configuration is less than unity ( $\nu _e^* \lt 1$ ), which is within the regime where the Sauter and Redl models are expected to provide reasonable approximations. While no single code serves as a definitive reference, NEO and SFINCS solve the drift-kinetic equation directly and provide high-fidelity neoclassical transport calculations. Given the complexity of bootstrap current calculations, agreement within a few percent at the profile peak is generally considered acceptable. Accordingly, the strong agreement among M3D-C1, NEO-3D, SFINCS and XGCa simulations supports the validity of the bootstrap current calculations for this case.

Figure 1. Set-up and results for tokamak verification: (a) prescribed density and temperature distributions; (b) bootstrap current profiles showing that M3D-C1 results closely approximate those from the drift-kinetic calculations. $^1$ Hager & Chang (Reference Hager and Chang2016), $^2$ Belli & Candy (Reference Belli and Candy2008, Reference Belli and Candy2011), $^3$ Landreman et al. (Reference Landreman, Smith, Mollén and Helander2014).

3.2. Stellarator verification

The implementation of Redl et al. (Reference Redl2021) formulae with the isomorphism as defined by Landreman et al. (Reference Landreman, Buller and Drevlak2022) are verified on two QA stellarator configurations from Landreman & Paul (Reference Landreman and Paul2022) and Landreman et al. (Reference Landreman, Buller and Drevlak2022). These QA configurations both have a minor radius of 1.70 m and volume-averaged ${{B}}=5.86\,{\rm T}$ . The density ( $n$ ) and temperature ( $T$ ) profiles for both cases are specified using

(3.1a) \begin{align} n(\psi _{t_N})&=n_0\big(1-\psi _{t_N}^5\big),\\[-10pt]\nonumber \end{align}

(3.1b) \begin{align} T(\psi _{t_N})&=T_0\big(1-\psi _{t_N}\big) ,\end{align}

where $\psi _{t_N}$ is the normalized toroidal flux. For the first case (QA_Case1), a pure hydrogen plasma with $n_0=n_{H,0}=4.13 \times 10^{20}$ m $^{-3}$ and $T_0=T_{H,0}=12$ keV (Landreman et al. Reference Landreman, Buller and Drevlak2022; Landreman et al. Reference Landreman, Buller and Drevlak2022 – § IV) is considered. This case uses a VMEC (Hirshman & Whitson, Reference Hirshman and Whitson1983) equilibrium with zero plasma pressure. VMEC is an equilibrium solver commonly used in stellarator optimization, which uses a variational method to find the minimum total energy of a system. The bootstrap current is calculated using the temperature and density profiles shown in figure 2(a), however, it is not self-consistent with the equilibrium. The second case (QA_Case2) is a QA configuration with $n_0=n_{e,0}=2.38 \times 10^{20}$ m $^{-3}$ , $T_0=T_{e,0}=9.45$ keV, that was optimized for volume-averaged $\beta =$ 2.5 $\%$ (see Sect. VIC in Landreman et al. Reference Landreman, Buller and Drevlak2022)). Here, $\beta$ is the ratio of the plasma pressure to the magnetic pressure. Figures 2 and 3 provide further details of each configuration including the equilibrium profiles, cross-sections at various toroidal angles and three-dimensional views.

Figure 2. Quasi-axisymmetric configuration (QA_Case1): (a) density and temperature equilibrium profiles, (b) toroidal cross-sections of the plasma boundary and (c) three-dimensional view.

Figure 3. Optimized QA configuration with a volume-averaged $\beta = 2.5\,\%$ (QA_Case2): (a) density and temperature equilibrium profiles, (b) toroidal cross-sections of the plasma boundary and (c) three-dimensional view.

Figure 4. Bootstrap current profiles for (a) quasi-axisymmetric configuration (QA_Case1) and (b) optimized QA configuration with volume-averaged $\beta = 2.5 \%$ . (QA_Case2). $^1$ Landreman et al. (Reference Landreman, Smith, Mollén and Helander2014), $^2$ Redl et al. (Reference Redl2021).

Figure 4 compares M3D-C1’s bootstrap current profile with those of the SFINCS and Redl formulae from Landreman et al. (Reference Landreman, Buller and Drevlak2022). It is evident from figure 4 that M3D-C1’s implementation of the modified Redl et al. (Reference Redl2021) formulae, as presented in Landreman et al. (Reference Landreman, Buller and Drevlak2022), is in close agreement with the SFINCS and Redl based on calculations from Landreman et al. (Reference Landreman, Buller and Drevlak2022). Minor discrepancies between the profiles may be attributed to numerical treatments between the codes or to temperature isotherms not being exactly aligned with the magnetic surfaces, although the exact sources remain unclear. Nevertheless, the results strongly support the robustness of M3D-C1’s implementation in these complex configurations.

3.3. Stellarator: nonlinear evolution

In this section, the impact of the bootstrap model on the nonlinear evolution of a QA stellarator equilibrium is analyzed. The simulations conducted in M3D-C1 are initialized using the QA equilibrium profile from section VIC of Landreman et al. (Reference Landreman, Buller and Drevlak2022) (QA_Case2). The nonlinear evolution of the equilibrium are begun from an initial state with nested magnetic surfaces. This equilibrium, developed using SIMSOPT optimization software (Landreman et al. Reference Landreman, Medasani, Wechsung, Giuliani, Jorge and Zhu2021), was optimized for low quasi-symmetry error, good energetic particle confinement and a self-consistent bootstrap current. The initial equilibrium profiles used in the simulations in this section are shown in figure 3(a). The computational domain is aligned with the shape of the equilibrium plasma, bounded by the last closed flux surface. The viscosity coefficient is fixed at ( $\nu = 3.65 \times 10^{-4}$ kg (ms)⁻¹). The heat transport is modeled according to (C5), which accounts for thermal energy exchange and diffusion across the plasma. To explore a range of physical scenarios, several simulation configurations are examined. Four resistivity profiles are considered, using the Spitzer resistivity model with varying scaling factors of $\eta _0=$ {1, 10, 1000, 10 000}. The general resistivity form is expressed as

(3.2) \begin{equation} \eta (R,\phi ,Z)=\eta _{\textit{norm}}\,\eta _0\,T^{-3/2}_e(R,\phi ,Z), \end{equation}

where $\eta _{\textit{norm}}=2.74$ $\varOmega$ m is the normalization factor for the classical resistivity value, $\eta _0$ is an artificial scale factor and $T_e$ is the electron temperature at the position $(R,\phi ,Z)$ . In this set of simulations, the perpendicular thermal conductivity ( $\kappa$ ) is held constant at $2.18\times 10^{-20}$ m $^{-1}$ s $^{-1}$ . To assess the effect of the bootstrap model, all configurations are evaluated with the bootstrap model both enabled and disabled.

The simulations are performed on a semi-structured grid consisting of 36 toroidal planes, resulting in a total of $1.88 \times 10^5$ three-dimensional elements within the computational domain. The grid is designed to provide sufficient resolution to accurately capture the nonlinear dynamics of the plasma.

Figure 5 shows the toroidal current density profiles from simulations with $\eta _0=10000$ at two times, $t=0$ and $t=250\tau _A$ , comparing the cases with and without the bootstrap current model. In the absence of the bootstrap model, the current decays over time. In contrast, when the bootstrap model is enabled, the current is maintained, emphasizing its role in sustaining plasma currents. However, despite these differences, the simulations remain MHD unstable, and the current continues to evolve. The plasma instability is further evidenced by the Poincaré plots in figure 6. At $t=250.0\tau _A$ in figure 6, these plots reveal that the chaotic region at the plasma boundary is significantly smaller when the bootstrap model is enabled, indicating that enabling the bootstrap model modifies the evolution of instabilities by slowing the breakup of magnetic surfaces.

Figure 5. Toroidal current density profiles for simulations with resistivity scale factor $\eta _0=10^4$ , comparing simulations with and without the bootstrap (BS) model at $t = 0$ and $250\tau _A$ , highlighting the effect of the bootstrap current model on the profile evolution.

Figure 6. Poincaré sections of the magnetic field starting from a QA stellarator equilibrium optimized at 2.5 % plasma beta for the case with resistivity scale factor $\eta _0=10^4$ (QA_Case2). The sections are shown as a function of time for $0 \leq t \leq 350\tau _A$ , comparing simulations with the bootstrap model disabled (top row) and enabled (bottom row).

It was found in Wright & Ferraro (Reference Wright and Ferraro2024a ) that the growth of instabilities in equilibria similar to the ones under consideration here exist in the limit of low resistivity (e.g. are essentially ideal), but are accelerated at high resistivity. Because the instabilities progressed more rapidly than the resistive decay time, it was concluded that this result was not a consequence of neglecting the bootstrap current drive; for example, through the spurious resistive decay of equilibrium currents that are actually driven by the bootstrap effect. Now that we are in a position to evaluate this with a bootstrap model included, we vary the resistivity scaling factor $\eta _0$ here, running simulations with $\eta _0=$ 1, 10, 1000 and 10 000. Figure 7 shows the kinetic energy grows more rapidly with increasing resistivity, supporting the conclusions of Wright & Ferraro (Reference Wright and Ferraro2024a ) that resistivity enhances the growth rate, even when the bootstrap model is included.

Figure 7. Time evolution of total kinetic energy for varying resistivity scaling factors $\eta _0 =$ 1, 10, 1000 and 10 000 demonstrating faster kinetic energy growth at higher resistivity.

To further investigate the stability of the equilibrium over time, the deviation from quasi-symmetry is quantified as the system evolves. This is achieved by calculating the two-term quasi-symmetry error ( $F_{QS}$ ) (Helander & Simakov Reference Helander and Simakov2008; Helander Reference Helander2014; Paul, Antonsen & Cooper Reference Paul, Antonsen, Landreman and Cooper2020), as shown in the following equation:

(3.3) \begin{equation} F_{QS}= \left\langle \left\{ \frac {1}{B^3}[(N-\iota M)\boldsymbol{ B } {\times} \boldsymbol{\nabla}B \boldsymbol{\cdot }\boldsymbol{\nabla }\psi - (MG+NI)\boldsymbol{ B }\boldsymbol{\cdot }\boldsymbol{\nabla}B] \right\}^2 \right\rangle \end{equation}

where $B$ is the magnetic field, $M=1,N=0$ for quasi-axisymmetry. The notation $\langle \boldsymbol{\cdot }\rangle$ denotes the magnetic surface average, with isotherms of $T_e$ serving as proxies for magnetic surfaces, consistent with the approach used in the bootstrap current calculation. Figure 8 compares $F_{QS}$ obtained from M3D-C1 outputs with the bootstrap model enabled, at two distinct times $t=0.0$ , $t=250.0\tau _A$ , as well as from SFINCS and VMEC using the $t=0$ equilibrium. The M3D-C1 $F_{QS}$ metric closely resembles the VMEC quasi-symmetry error. While the equilibrium evolves over time, as observed in the Poincaré plots, the $F_{QS}$ value does not significantly change in the range $0\lt t\lt 250$ $\tau _A$ , indicating that quasi-symmetry is maintained throughout this period. However, beyond $t\gt 250$ , the chaotic region grows, making it difficult to generate well-defined isothermal surfaces globally (at which point neoclassical transport is likely strongly subdominant to parallel transport along chaotic field lines anyway), so analysis is concluded there.

Figure 8. Two-term quasi-symmetry error as defined in (3.3) for the QA stellarator equilibrium optimized at 2.5 % plasma beta (QA_Case2) with a resistivity scale factor $\eta _0=10^4$ .

4. Summary and discussion

This work enhances the capabilities of the M3D-C1 extended-MHD code by incorporating self-consistent bootstrap current models for tokamak and quasi-symmetric stellarator geometries. Two models are implemented, namely, (a) the Sauter model (Sauter et al. Reference Sauter, Angioni and Lin-Liu1999) and (b) the Redl model (Redl et al. Reference Redl2021). For quasi-symmetric stellarators, the isomorphism outlined by Landreman et al. (Reference Landreman, Buller and Drevlak2022) is employed.

The numerical verification presented in this work demonstrates the accuracy of these new implementations. In § 3.1, the bootstrap current in a tokamak with a circular plasma boundary is tested using both bootstrap models. Comparisons between the flux-surface-averaged bootstrap current from the updated M3D-C1 code and results from established neoclassical codes, such as XGCa (Hager & Chang Reference Hager and Chang2016), NEO (Belli & Candy Reference Belli and Candy2008, Reference Belli and Candy2011) and SFINCS (Landreman et al. Reference Landreman, Smith, Mollén and Helander2014), show excellent agreement.

In § 3.2, the bootstrap current is evaluated in two quasi-symmetric stellarator configurations, QA_Case1 and QA_Case2. These results are compared with SFINCS calculations, yielding good agreement. Finally, § 3.3 explores the application of the updated M3D-C1 code to nonlinear MHD simulations of an optimized QA stellarator configuration, where the expected toroidal current sustainment is observed.

The accurate modeling of bootstrap currents is critical to understanding and improving plasma performance in magnetic confinement devices. These results not only verify the bootstrap models in the updated M3D-C1 code but also highlight its potential to improve the design and optimization of future fusion reactors through accurate neoclassical current predictions in non-axisymmetric geometries. The successful implementation and verification of these models provides a solid foundation for advancing nonlinear MHD simulations, offering crucial insights that can guide the design and operational strategies of next-generation fusion devices.

This capability will enable new investigations into nonlinear MHD physics where bootstrap current is expected to be important. In particular, this includes transport-time scale simulations where the global equilibrium evolves due to heating, as is done in nonlinear calculations of pressure limits (cf. Wright & Ferraro Reference Wright and Ferraro2024b ), sawtooth cycles (cf. Jardin et al. Reference Jardin, Ferraro, Breslau and Chen2012b ) and ELM cycles (cf. Futatani et al. Reference Futatani, Cathey, Hoelzl, Lang, Huijsmans and Dunne2021), for example. While the implementation described here is only strictly valid for quasi-symmetric magnetic geometries, the ability to calculate the neoclassical transport of generic non-axisymmetric output using NEO 3D (Sinha et al. Reference Sinha, Ferraro and Belli2022, Reference Sinha, Ferraro and Belli2023) and SFINCS (as described above) enable post hoc evaluation of the accuracy of the Sauter–Redl–Landreman model. These capabilities also lay the groundwork for coupled M3D-C1/NEO 3D or M3D-C1/SFINCS modeling, which would be applicable to a wider range of geometries. One fundamental challenge is that the treatment of stochastic regions is beyond the scope of existing neoclassical theory, and requires a kinetic-MHD theory to treat self-consistently. In cases where the accuracy of the dynamics of these regions is less critical, such as when the time evolution is slow or time-dependent effects are negligible, the flattening of temperature and density profiles by classical transport along stochastic field lines, which is naturally included in the extended-MHD model, is expected to dominate over neoclassical effects in any case. We also note that the Sauter–Redl–Landreman models implemented here, as well as NEO 3D and SFINCS, all calculate time-independent transport, and are therefore not able to treat time-dependent neoclassical response. In principle this limits the applicability of the model on Alfvénic time scales (e.g. kink modes) – compounded by the fact that these modes are not likely to maintain quasi-symmetry – but the neoclassical response is not expected to play a significant role in such phenomena anyway.