Direct prediction of saturated neoclassical tearing modes in slab using an equilibrium approach

Abstract

We demonstrate for the first time that the nonlinear saturation of neoclassical tearing modes (NTMs) can be found directly using a variational principle based on Taylor relaxation, without needing to simulate the intermediate, resistivity-dependent dynamics. As in previous investigations of classical tearing mode saturation (Loizu et al., 2020; Loizu & Bonfiglio, 2023), we make use of SPEC (Hudson et al., 2012), an equilibrium solver based on the variational principle of the Multi-Region relaxed MHD, featuring stepped pressure profiles and arbitrary magnetic topology. We work in slab geometry and employ a simple bootstrap current model $J_{bs}=C\>\bnabla p$ to study the bootstrap-driven tearing modes, scanning over the asymptotic matching parameter $\Delta^{\prime}$ and bootstrap current strength. Saturated island widths produced by SPEC agree well with the predictions of an initial value resistive MHD code (Huang & Bhattacharjee, 2016) while being orders of magnitude faster to calculate. Additionally, we observe good agreement with a simple analytical Modified Rutherford Equation, without requiring any fitting coefficients. The match is obtained for both linearly unstable classical tearing modes in the presence of bootstrap current, and neoclassical tearing modes, which are linearly stable but nonlinear-unstable due to the effects of the bootstrap current.

1 Introduction

Tearing modes convert magnetic field energy into kinetic energy of the plasma through spontaneous magnetic reconnection (Goldston, 2020). Generally present in toroidally-confined magnetic fusion plasmas, these internal modes violate the frozen-in constraint of ideal MHD, modifying the flux surface topology and creating magnetic islands. In many cases, tearing modes do not lead to the complete collapse of the plasma, but rather grow into non-linearly saturated states (Kong, 2020) that can strongly increase global transport and lower performance (Günter et al., 1999). These instabilities often appear in the form of neoclassical tearing modes (NTM), the behavior of which is driven by the bootstrap current. The latter is generated by the collisional interaction of trapped and passing particles in a toroidal plasma and is typically proportional to pressure gradients (Helander & Sigmar, 2005). As the transport of particles and heat along the magnetic field dominates that across the field, islands tend to locally flatten pressure (Fitzpatrick, 1995), which leads to the removal of bootstrap current inside the island, driving the tearing mode unstable. Effects of bootstrap current on tearing modes include an enhanced growth of linearly unstable tearing modes and, more importantly, a nonlinear destabilization of configurations that are linearly stable to tearing modes.

As a first approximation, tearing modes are studied using linear theory which involves tearing-like perturbations of an initial equilibrium. These are constructed by solving the linearized resistive MHD equations in a narrow layer around the resonant surface, where the safety factor $q$ is rational, and performing an asymptotic matching with a solution of the linearized ideal MHD equations outside the layer. If a mode can be found that lowers the MHD potential energy, $\mathcal{W}=\int dV\>[B^{2}/2\mu_{0}+p/(\gamma-1)]$, then the equilibrium is deemed unstable to this mode. In the large aspect ratio limit, an important quantity is $\Delta^{\prime}$ which comes from the discontinuity of the derivative of the perturbed flux across the resistive layer for a given perturbation. Neglecting the stabilizing effects of pressure (Glasser et al., 1976), if $\Delta^{\prime}>0$, then a particular tearing mode is unstable and will lower $\mathcal{W}$, and if $\Delta^{\prime}<0$ the mode will increase $\mathcal{W}$ as it takes more energy to bend the fieldlines than what is released by reconnection. The linear approach provides a good start for uncovering the physics of tearing modes, but it only yields information on linear growth and stability. Therefore, it is of limited use for understanding nonlinear growth and saturation (Goldston, 2020), as is measured in experiments (Kong, 2020).

To study the nonlinear dynamics, Rutherford (1973) started from the linear theory and extended it into a nonlinear regime, deriving an ODE for the evolution of the island width $w$, which is a good measure of the mode amplitude and also relates to the level of confinement degradation of the underlying equilibrium. Fitzpatrick (1995) extended Rutherford’s model by including the effects of the bootstrap current. An example of such ODE is:

$$\frac{\mu_{0}}{\eta}\frac{\partial w}{\partial t}=\Delta^{\prime}_{0}+\Delta^{\prime}_{\textrm{stab}}(w)+\Delta^{\prime}_{\textrm{bs}}(w)$$ (1)

where $\Delta^{\prime}_{0}$ is the value of $\Delta^{\prime}$ for perturbations of an initial equilibrium with $w=0$, $\Delta^{\prime}_{\textrm{stab}}(w)\leq 0$ includes stabilizing terms and ensures saturation of the island, and $\Delta^{\prime}_{\textrm{bs}}(w)$ accounts for the effects of the bootstrap current. The bootstrap term is destabilizing and can be found analytically in the large aspect ratio limit, $\Delta^{\prime}_{\textrm{bs}}(w)\sim 1/w$ (Fitzpatrick, 1995). This term diverges for small islands and so is modified to $\Delta^{\prime}_{\textrm{bs}}(w)\sim w^{2}/(w^{2}+w_{0}^{2})$, thus introducing a regularizing scale, $w_{0}\sim({\chi_{\perp}}/{\chi_{\parallel}})^{1/4}$, which describes the effects of incomplete pressure flattening in an island. Here $\chi_{\perp}$ and $\chi_{\parallel}$ are the heat diffusion coefficients across and along the magnetic field, respectively. The typical evolution of an island for a linearly stable equilibrium ($\Delta^{\prime}_{0}<0$) using either of the two bootstrap terms is shown in Figure 1. When $w_{0}$ isn’t included, the mode is always unstable and will grow regardless of the initial island size. However, when $w_{0}$ is included and incomplete flattening is accounted for, then the mode will not grow unless it has been perturbed with a sufficient seed, i.e. if the initial island size is sufficiently large. Seeding dynamics are complicated (Muraglia et al., 2021), but it is important to note that, in the large island limit $w\gg w_{0}$ , the saturated island is essentially the same, regardless of whether the seeding mechanism is included or not. For the purposes of this paper, we study the $w\gg w_{0}$ regime, which is one of greater relevance to magnetic fusion experiments (in a typical tokamak $w_{0}\approx 1\textrm{cm}$ while $w_{\textrm{sat}}\approx 4\textrm{cm}$) (Sauter et al., 2002). Non-linear equations such as Eq. 1 are known as Modified Rutherford equations (MRE) and are often used in tokamak experiments (Kong, 2020) as they are quick to evaluate and give good estimates for the saturated island widths. However, MRE often includes other $\Delta^{\prime}(w)$ ad-hoc terms (Kong et al., 2022) which are device-specific and can feature many fitting coefficients.

Figure 1: Sketch of the width growth rate as a function of width, as described by Equation 1. The black curve shows the behavior with the diverging $1/w$ term, while the red curve demonstrates the effect of the $w_{0}$ correction to $\Delta^{\prime}_{\textrm{bs}}$.

A more accurate approach to studying nonlinear tearing modes is the time integration of resistive MHD equations, evolving a perturbed equilibrium into a steady state featuring a saturated island. Such a method gives a detailed description of tearing mode dynamics and can be expanded to add various non-MHD effects, but it is numerically very expensive (Hoelzl et al., 2021; Teng et al., 2018). The cost comes from the large separation between ideal and resistive timescales (measured by Lundquist number $S=\tau_{\textrm{res}}/\tau_{\textrm{Alfven}}\gg 1$), which requires long simulation times, and the separation between system size and tearing spatial scales, which requires fine meshes. In summary, given existing approaches, it is difficult to predict the saturation of tearing modes in a manner that is simultaneously fast, geometry-agnostic, and rooted in solid theory.

However, there may exist another approach based on energy principles that can remedy some of these stated drawbacks. The method was pioneered by Cooper et al. (2010), who used the ideal MHD equilibrium solver VMEC (Hirshman & Whitson, 1983) to predict the saturated state of internal kink modes in a tokamak. The idea is that, in addition to the initial MHD equilibrium, one can think of the saturated state of an instability as yet another (lower energy) MHD equilibrium. Thus, if an equilibrium code exists where both the initial and saturated states are in the space of possible solutions, then by clever perturbations in the parameter space it may be possible to directly jump from the initial to the final state. In addition to ideal modes such as the internal kink, this approach proved valuable for finding nonlinear states of non-resonant modes (Brunetti et al., 2014), where $q=m/n$ surface is not present. Here, nonlinear states from VMEC (Hirshman & Whitson, 1983) agreed well with a resistive MHD code when a mode was nonresonant, but led to a poor agreement otherwise.

Loizu extended the method for resonant resistive instabilities. His work showed that such an approach can be used for the prediction of the saturated state of classical tearing modes, both in the slab (Loizu et al., 2020) and cylindrical geometry (Loizu & Bonfiglio, 2023). Tearing mode prediction required using an equilibrium code such as the Stepped Pressure Equilibrium Code (SPEC) (Hudson et al., 2012), which does not constrain the magnetic topology around the resistive layer. In this paper, we present an extension of the work done in slab (Loizu et al., 2020) where the effects of bootstrap current are included and lead to the growth and saturation of neoclassical tearing modes. The paper is organized as follows. In section 2 we introduce the set of tearing-unstable initial equilibria, including the geometry and starting profiles. Section 3 describes how the saturation of tearing modes is predicted using resistive MHD initial value simulations. Finally, section 4 presents the study of nonlinear saturation using SPEC, along with a comparison of the variational approach to the resistive MHD and MRE.

2 Family of tearing unstable equilibria

To demonstrate the new approach for the prediction of tearing mode saturation with bootstrap current, we work in slab geometry. This simplifies the analysis since the slab does not have curvature or a coordinate singularity. In addition, we can set up an equilibrium that is stable to ideal modes and that has only one resonant surface where a tearing mode can grow unstable. Furthermore, analytical linear and nonlinear theories are simpler in slab geometry. We use a cartesian coordinate system where $x$ is a radial-like coordinate, $y$ is the poloidal-like direction, $z$ is the toroidal-like direction, and the slab is doubly periodic in $y$ and $z$. The size of the slab is fixed to $2\pi$ along $x$ and $z$, while we consider different values $L$ of the system size along $y$, thus allowing us to scan different stability regimes, with effects similar to varying the safety factor $q$ in toroidal geometry. An initial equilibrium magnetic field $\boldsymbol{B}_{0}=B_{z0}\hat{z}+\hat{z}\times\bnabla\psi_{0}$ is constructed with a flux function $\psi_{0}(x)=1/\cosh^{2}{(x)}+\alpha x^{2}$, similar to that in Loureiro et al. (2005). Associated to this flux function is the poloidal field, $B_{y0}(x)=\psi_{0}^{\prime}(x)$, and the toroidal current, $J_{z0}(x)=\psi_{0}^{\prime\prime}(x)$, which we call the ohmic current. The initial pressure $p_{0}(x)$ is simply a linear function of $x$, with the scale being given by the condition on the average plasma beta $\beta=\langle 2\mu_{0}p/B^{2}\rangle\approx 3\%$. The equilibrium, profiles of which are shown in Fig. 2, is ideally stable and has a single resonant surface at $x=0$, and $B_{z0}$ is an essentially constant strong guide field, such that the mode is approximately incompressible (Goldston, 2020), with minor variations to set up an initial force balance, $B_{z0}^{2}+|\bnabla\psi_{0}|^{2}+p_{0}=\textrm{const}$. The field at the resonant surface is solely in the $z$ direction, which means that the mode, in order to have the resonance condition $\boldsymbol{k}\cdot\boldsymbol{B}=0$, is symmetric in the $z$ direction. The bootstrap current, which in a torus originates from the field curvature, is not inherently present in a slab equilibrium. We overcome this by employing a simple model for the bootstrap current $\boldsymbol{J}_{\textrm{bs}}=\hat{C}\>C_{0}(\partial p/\partial x)\>\hat{z}$, where $C_{0}=J_{z0}(0)/(\partial p_{0}(0)/\partial x)$ and $\hat{C}$ is a free constant which specifies the relative strength of the bootstrap and ohmic current in the initial equilibrium. This term has been successfully used in previous studies of slab NTMs (Muraglia et al., 2021). Thanks to the POEM theory (Escande & Ottaviani, 2004; Militello & Porcelli, 2004), the stabilization term of MRE, Eq. 1, can be expressed analytically as $\Delta^{\prime}_{\textrm{stab}}(w)=-(1/2.44d^{2})w$, where $d$ is the length scale of the initial equilibrium current, $d=\sqrt{J_{z0}(0)/J_{z0}^{\prime\prime}(0)}\approx 0.35$. Using the POEM stabilization term and a simplified bootstrap contribution $\Delta^{\prime}_{bs}=\hat{C}\>C_{0}\>p^{\prime}_{r}\>b\>(1/w)$, where $p^{\prime}_{r}$ is the pressure gradient across the resonant surface and $b$ is a fitting coefficient (which can be estimated analytically, see Appendix B), we arrive at the following expression

$$w_{sat}=1.22d^{2}\left(\Delta^{\prime}_{0}+\sqrt{\Delta_{0}^{\prime 2}+1.64a^{-2}\>\hat{C}\>C_{0}\>p^{\prime}_{r}\>b}\right)$$ (2)

This equation provides an estimate of the island width for tearing modes with small constant $J_{\textrm{bs}}$ and is valid only in the large island limit, $w\gg w_{0}$.

3 Resistive MHD

The canonical model equations describing the nonlinear evolution of tearing modes are the resistive MHD equations:

$$\frac{\partial\rho}{\partial t}+\bnabla\boldsymbol{\cdot}\left(\rho\boldsymbol{u}\right)=0$$ (3)

$$\frac{\partial(\rho\boldsymbol{u})}{\partial t}+\bnabla\boldsymbol{\cdot}\left(\rho\boldsymbol{u}\boldsymbol{u}\right)=-\bnabla\left(p+\frac{1}{2}B^{2}\right)+\bnabla\boldsymbol{\cdot}\left(\boldsymbol{B}\boldsymbol{B}\right)+\bnabla\boldsymbol{\cdot}\left(\mu\rho\boldsymbol{\mathsfbi{\epsilon}}$$

$$\frac{\partial p}{\partial t}+\bnabla\boldsymbol{\cdot}\left(p\boldsymbol{u}\right)=(1-\gamma)\left(p\bnabla\boldsymbol{\cdot}\boldsymbol{u}-\frac{\chi_{\parallel}}{2}\bnabla\boldsymbol{\cdot}\left(\rho\boldsymbol{b}\boldsymbol{b}\boldsymbol{\cdot}\bnabla(p/\rho)\right)\right)$$ (4)

$$\frac{\partial\psi}{\partial t}=(\boldsymbol{u}\times\boldsymbol{B})_{z}+\eta(J_{z0}-J_{z})+\eta J_{\textrm{bs}}$$ (5)

These equations are solved, in a dimensionless form, using the code HMHD (Huang & Bhattacharjee, 2016). Viscous forces are introduced in the momentum balance, Eq. 3.2, through the strain tensor $\boldsymbol{\mathsfbi{\epsilon}}=(\bnabla\boldsymbol{u}+\bnabla^{T}\boldsymbol{u})/2$. In Ohm’s law, Eq. 3.4, $\eta J_{z0}$ is the loop voltage which maintains a non-zero steady state, and $\eta J_{\textrm{bs}}$ acts like a secondary driving term, allowing for the inclusion of an effective bootstrap current, as used by Muraglia et al. (2021). The variables are normalized as follows (real $\rightarrow$ normalized): length $x,y,z\rightarrow x,y,z/L_{0}$, density $\rho\rightarrow\rho/\rho_{0}$, magnetic field $\boldsymbol{B}\rightarrow\boldsymbol{B}/B_{0}$, current $\boldsymbol{J}\rightarrow\boldsymbol{J}/(B_{0}/L_{0}\mu_{0})$, velocities $\boldsymbol{u}\rightarrow\boldsymbol{u}/u_{\textrm{A}}$ (with Alfven speed $u_{\textrm{A}}=B_{0}/\sqrt{\mu_{0}\rho_{0}}$), time $t\rightarrow t/\tau_{\textrm{A}}$ (with Alfven time $\tau_{\textrm{A}}=L_{0}/u_{\textrm{A}}$), and pressure $p\rightarrow p/(B_{0}^{2}\mu_{0})$. The parameters of the model are the adiabatic coefficient $\gamma$, which we fix to that of the ideal gas $\gamma=5/3$, parallel diffusion coefficient, $\chi_{\parallel}\approx 10$, used for numerical stability, dynamic viscosity $\mu$, and resistivity $\eta$. The latter two are specified in terms of their respective dimensionless numbers, the Lundquist number $S=L_{0}u_{\textrm{A}}/\eta$, and the magnetic Prandtl number $Pr_{m}=\mu/\eta$. For the simulations shown below, the dimensionless numbers are $S\in[10^{3},10^{4}]$ and $Pr_{m}=[10^{-1},10^{1}]$. The code relies on the 5-point finite difference scheme for the spatial discretization, and the trapezoidal leapfrog method for the temporal discretization. The simulations are performed using dynamic time steps based on the CFL condition (Courant et al., 1967) for Alfvén dynamics. The mesh comprises approximately 400x400 points, periodic boundary conditions are applied in the $y$ direction, and wall boundary conditions are applied in the $x$ direction.

Simulating a tearing mode entails setting up the initial equilibrium, Figure 2, according to the profiles derived from the initial flux function $\psi_{0}$ and pressure $p_{0}$. The initial density is constant and set to $\rho\in[2,8]$ in normalized units. Note that the initial equilibrium does not include bootstrap current. The next step involves perturbing the initial flux with $\psi_{p}=\hat{\psi}_{p}\cos{(y/L)\cos{(x/2\pi)}}$, where $\hat{\psi}_{p}$ is the perturbation magnitude (the exact value does not affect saturation, see Fig. 3). The mode then grows spontaneously and reaches, on an approximately resistive time scale, $\tau_{R}={\mu_{0}L^{2}}/{\eta}=S\tau_{\textrm{Alfven}}$, a state where the profiles and island size stabilize. The time trace of the width of a typical island is depicted in Figure 3. The final step is to fit a decaying exponential of the form $w(t)=c_{0}+c_{1}\exp(c_{2}t)$ to the final stage of the island evolution. The best fit for the constant coefficient gives the saturated island width, $w_{\textrm{sat}}:=w(t\rightarrow\infty)=c_{0}$. Fig. 3 also shows that the value of $w_{\textrm{sat}}$ is independent of $\eta$, $\mu$, as well as the mesh resolution. The structure of the saturated island is shown in Figure 4a, where an $m=1$ mode can be seen to dominate. This supports an underlying assumption in MRE, which is that the nonlinear profiles have the same spatial dependence as the linear ones and are of single harmonicity.

Figure 2: Initial profiles ($\alpha=-0.002$) of the poloidal field, toroidal current, and pressure, as used in HMHD (red) and SPEC (blue). SPEC equilibrium consists of $N_{vol}=151$ volumes, the largest of which is the central volume, located around $x=0$.

Figure 3: Time evolution of the island width from a seed to a saturated island, as a function of the resistive time. Here the equilibrium has $L=10.7d$ and $\Delta^{\prime}_{0}\>d=0.97$. Different curves correspond to (1) a reference case (solid black), (2) a smaller seed (blue triangles), (3) smaller resistivity (orange circles), (4) finer mesh (green crosses), and (5) smaller viscosity (red stars).

Figure 4: Separatrices (red) and flux contours (black) of saturated islands in HMHD (left) and SPEC (right) for an equilibrium with $L=7.8d$, which corresponds to $\Delta^{\prime}_{0}\>d=-0.23$. The SPEC plot also depicts a selection of the ideal interfaces (green), the shape of which resembles the flux contours in HMHD.

The islands shown in Figs. 3,4 are for a slab of a certain size $L$, to which corresponds a value of $\Delta^{\prime}_{0}$. It is interesting to study the dependence of $w_{\textrm{sat}}$ on $\Delta^{\prime}_{0}$. Such a study is shown in Figure 5, where the island width is depicted as a function of the linear stability parameter $\Delta^{\prime}_{0}$. The black triangles reproduce the original study by Loizu et al. (2020), showing that classical modes ($J_{\textrm{bs}}=0$) do not have a saturated state if linearly stable ($\Delta^{\prime}_{0}<0$) and do have a saturated state, with $w_{\textrm{sat}}$ approximately linear in $\Delta^{\prime}_{0}$, if linearly unstable ($\Delta^{\prime}_{0}>0$). However, this is no longer true for tearing modes with bootstrap current, simulations of which are shown as black crosses (for $\hat{C}=0.048$). In this case, linearly stable modes ($\Delta^{\prime}_{0}<0$) can in fact grow into saturated islands. These are usually referred to as neoclassical tearing modes (NTM) and are commonly observed in tokamaks (Kong, 2020). As depicted in Figure 1, NTMs need to be initiated with a sufficiently large seed, otherwise islands collapse back into the initial equilibrium. The precise size of this seed, however, does not impact the saturated island (see Fig. 3). For linearly unstable equilibria ($\Delta^{\prime}_{0}>0$), the bootstrap current again acts in a destabilizing manner, resulting in an enhanced island width as compared to the cases with $J_{\textrm{bs}}=0$.

Figure 5: Scan of island width as a function of the linear stability $\Delta^{\prime}_{0}$. Classical modes are shown from SPEC (green circle), HMHD (black triangle). Bootstrap-enhanced classical modes and linearly stable neoclassical tearing modes are shown from SPEC (blue diamond) and HMHD (black cross). The analytical prediction from MRE is also displayed (red dash). The bootstrap coefficient used is $\hat{C}=0.048$.

4 Variational approach

For directly predicting the saturation of tearing modes, we use SPEC (Hudson et al., 2012), a code that finds Multi Region Relaxed MHD (MRxMHD) equilibrium states using a variational principle. Taylor (1974) hypothesized that resistive plasmas often evolve in a way that conserves global magnetic helicity $\mathcal{H}=\int_{V_{p}}dV\boldsymbol{A}\boldsymbol{\cdot}\boldsymbol{B}$ while minimizing plasma energy $\mathcal{W}$. Here $\boldsymbol{A}$ is the vector potential and $V_{p}$ is the plasma volume. Hole et al. (2006) extended this idea to build the MRxMHD model that subdivides a plasma into $N_{\textrm{vol}}$ constant-pressure Taylor-relaxed volumes $\mathcal{V}_{l}$ which are separated by ideal interfaces $\mathcal{I}_{l}$. Thus SPEC solves MRxMHD by finding the extrema of the energy functional

$$\mathcal{F}=\mathcal{W}+\mu(\mathcal{H}-\mathcal{H}^{0})=\sum_{l}\int_{\mathcal{V}_{l}}dV\>\left(\frac{B^{2}}{2\mu_{0}}+\frac{p_{l}}{\gamma-1}\right)+\mu_{l}\left(\int_{\mathcal{V}_{l}}dV\>\boldsymbol{A}\boldsymbol{\cdot}\boldsymbol{B}-\mathcal{H}^{0}_{l}\right)$$

where $\mathcal{H}^{0}_{l}$ are the target helicities and $\mu_{l}$ are the Lagrange multipliers. In slab geometry, the interfaces $\mathcal{I}_{l}$ are represented as

$$\boldsymbol{x}_{l}(\theta,\zeta)=\sum_{m}x_{l,m}\cos{(m\theta)}\>\hat{x}+L\frac{\theta}{2\pi}\>\hat{y}+\zeta\>\hat{z}$$

where $x_{l,m}$ are the coefficients of a Fourier expansion that is numerically truncated at $m=m_{\textrm{pol}}$. The field within each subvolume is found through minimization with respect to the vector potential $\boldsymbol{A}$, and force balance across the interfaces is achieved by varying the interface harmonics $x_{l,m}$. SPEC finds the states that extremize $\mathcal{F}$, which satisfy (Hudson et al., 2012)

$$\bnabla\times\boldsymbol{B}=\mu_{l}\boldsymbol{B}$$ (6)

$$\boldsymbol{F}_{l}=[p+B^{2}/2\mu_{0}]_{l}=0$$ (7)

where the first equation describes linear force-free fields inside each $\mathcal{V}_{l}$, and the second equation represents the force balance across the two sides of each $\mathcal{I}_{l}$. The position and shape of each interface is in fact determined by Eq. 4.2, which is the local equivalent of the MHD force balance equation $\boldsymbol{J}\times\boldsymbol{B}=\bnabla p$. The standard approach for finding minimal MRxMHD states uses a Newton method for finding the zero of $\boldsymbol{F}_{l}$ and the alternative is a gradient descent algorithm for minimizing the functional $\mathcal{F}$. We found that the Newton method alone is not robust enough for equilibria with islands. Instead, a more robust approach involves first minimizing the energy functional using a descent method, which is slower but more stable, and then using the Newton method to make fine adjustments to interface harmonics, resulting in an equilibrium with good force balance. The setup for solving MRxMHD in slab geometry with SPEC for a case with $N_{\textrm{vol}}=3$ is shown in Figure 6.

Figure 6: The setup for SPEC in slab geometry for a case of $N_{\textrm{vol}}=3$. SPEC solves for Beltrami fields within each volume and seeks force balance across ideal interfaces by perturbing the shapes of these interfaces.

Inputs to SPEC include four physical constraints per volume that are fixed during functional minimization. These are the pressure $p_{l}$, the integrated toroidal and poloidal fluxes, $\Delta\Psi^{t}_{l}$ and $\Delta\Psi^{p}_{l}$, and the helicities $\mathcal{H}_{l}$ inside each volume. Alternatively, SPEC can directly solve Eqs. 4.1 and 4.2 provided a different set of four constraints, ($p_{l}$, $\Delta\Psi^{t}_{l}$, $\Delta\Psi^{p}_{l}$, $\mu_{l}$) or ($p_{l}$, $\Delta\Psi^{t}_{l}$, $I^{\textrm{vol}}_{l}$, $I^{\textrm{surf}}_{l}$), where $I^{\textrm{vol}}_{l}$ is the integrated toroidal current in each volume and $I^{\textrm{surf}}_{l}$ is the sheet current on each ideal interface (Baillod et al., 2021).

While the linear growth of tearing modes occurs on a time scale that is intermediate between Alfvenic, $\tau_{\textrm{Alfven}}$, and resistive, $\tau_{R}$ (for slab $\tau=\tau_{\textrm{Alfven}}^{2/5}\tau_{R}^{3/5}$ (Goldston, 2020)), the nonlinear saturation occurs on an approximately resistive time scale, as is demonstrated in Figure 3. Therefore, Taylor’s hypothesis of global helicity conservation, which works well for faster relaxation events (Berger, 1999), is not valid for slow tearing modes. Instead, the POEM theory (Militello & Porcelli, 2004) reveals that a good invariant is the flux surface average of the toroidal current $J_{z}(\psi)=\langle J_{z0}\rangle_{\psi}$, with $\psi$ being the nonlinear flux function. See Loizu et al. (2020) for a numerical demonstration of this invariant and its comparison with helicity. To construct equilibria where an integrated version of $J_{z}(\psi)=\langle J_{z0}\rangle_{\psi}$ is fixed, we run SPEC with the constraints ($p_{l}$, $\Delta\Psi^{t}_{l}$, $I^{\textrm{vol}}_{l}$, $I^{\textrm{surf}}_{l}$).

As with HMHD, finding tearing modes in SPEC involves first setting up the equilibrium. The initial domain is covered along the $x$ direction with ideal interfaces, with a number $N_{\textrm{vol}}$ that is sufficient to reproduce the equilibrium profiles. Special care needs to be given to the placement of the two interfaces bounding the resonant surface $x=0$. If these are placed too closely to the resonant surface then they might limit the growth of the island. Conversely, if the spacing between them is too large then the analytical profiles are not well represented around $x=0$, and can be artificially stabilized. The arbitrariness in interface placement is removed by placing the two interfaces symmetrically around $x=0$ with a spacing of $\Psi_{w}$ (defined as the normalized toroidal flux enclosed by these two interfaces) which maximizes the width of the resulting saturated island. For more details on interface placement, see Loizu et al. (2020). From the interface positions and the analytical equilibrium profiles, the SPEC input profiles including $\Delta\Psi^{t}_{l}$, $I_{l}^{\textrm{vol}}$, and $p_{l}$ are calculated. Since the pressure in SPEC is stepped, the bootstrap model $J_{\textrm{bs}}=\hat{C}\>C_{0}\>\partial p/\partial x$ is implemented as a set of sheet currents on the interfaces $I_{l}^{\textrm{surf}}=\hat{C}\>C_{0}\>\Delta p_{i}$. The integrated bootstrap current in SPEC is the same as in HMHD, and, as we will show, equivalent saturated islands can be obtained despite the discreteness of the current sheets. Note that the bootstrap current in SPEC, unlike in HMHD, exists already in the initial equilibrium, which explains why the initial toroidal volume currents, shown in Fig. 2, are not exactly the same between the two codes. SPEC is run with these inputs to find the initial MRxMHD state, profiles of which are shown in Fig. 2.

The next step is to evaluate linear stability and determine the unstable eigenmode with which the initial equilibrium can be perturbed, before a secondary equilibrium with an island is sought. This is accomplished by analyzing the Hessian matrix (Kumar, 2022) defined as

$$\mathsfbi{H}_{ij}=\frac{\partial\mathcal{F}}{\partial\boldsymbol{X}_{i}\>\partial\boldsymbol{X}_{j}}$$

where $i,j$ index the packed degrees of freedom of the interfaces $\boldsymbol{X}=\{x_{l,m}\>|\>l\!=\!1,...,N_{\textrm{vol}};\>m\!=\!0,...,m_{\textrm{pol}}\}$. Eigendecomposition of $\mathsfbi{H}$ gives information on linear stability, where a negative eigenvalue implies instability and the associated eigenmode describes the shape of the mode in the space of interface Fourier coefficients. Since the initial equilibrium is ideally stable, depending on $L$, there is either zero or one negative eigenvalue, the eigenmode of which has the shape of an odd-parity tearing-like mode. Linear stability in SPEC has been verified for classical tearing modes (Loizu et al., 2020), while the stability of modes with bootstrap current is more subtle. In fact, as discussed in Sec. 1 and illustrated in Fig. 1, a system for which $w_{0}\rightarrow 0$ (which corresponds to perfect pressure flattening in the island) is unconditionally unstable. Since the pressure profile in SPEC is, by definition, flat in the relaxation region, we expect that tearing modes will always be unstable when bootstrap effects are included. This is in fact what we observe when looking at the Hessian (see Appendix A). Nevertheless, as discussed in Sec. 1 and seen in Fig. 1, the saturation of an island is essentially independent of $w_{0}$ in the limit $w\gg w_{0}$ studied here. After finding the eigendecomposition of the Hessian, the next step involves perturbing the interfaces of the SPEC equilibrium, creating the seed island. The perturbation is constructed from the unstable eigenmode which is scaled to give the largest possible seed that does not overlap the ideal interfaces. Finally, SPEC is run for the second time, iterating on the interface harmonics to reach a secondary MRxMHD equilibrium with a lower energy and a topology of an island.

The results of running SPEC for a range of equilibria can be seen in Figure 5. The green circles correspond to classical tearing modes without bootstrap current, as reported in Loizu et al. (2020). In this case, SPEC recovers the island width obtained using HMHD as well as the prediction by the POEM theory (Militello & Porcelli, 2004). As is shown, SPEC correctly predicts the linear stability threshold such that stable equilibria for which $\Delta^{\prime}_{0}<0$ are also nonlinearly stable. The blue diamonds in Fig. 5 show the results from running SPEC with the inclusion of bootstrap current ($\hat{C}=0.048$). SPEC again manages to predict islands both in the case of classical modes (enhanced with bootstrap) and neoclassical tearing modes (which are only destabilized by the bootstrap current effect). A slight difference exists in the width between SPEC and HMHD for the larger islands, which we believe may be due to a violation of the current invariant with the increasing island size. Also shown in the figure is a prediction of an MRE, Eq. 2.1, which includes the bootstrap contribution. The results related to the MRE rely on a constant that we fit numerically, $b=4.25$, as is standard for these models in general (Kong, 2020), instead of using the less-accurate, analytical value of $b=8$, as derived in Appendix B. The fit is obtained by matching the island width for the equilibrium with $\Delta^{\prime}_{0}=0$.

The analysis in Figure 5 was done using a fixed bootstrap strength, $\hat{C}=0.048$. To verify that the prediction from SPEC is independent of the exact value of $\hat{C}$, we chose two equilibria, one linearly stable, $\Delta^{\prime}_{0}\>d=-0.26$, and one linearly unstable, $\Delta^{\prime}_{0}\>d=0.26$, and we vary the constant $\hat{C}$ from $0$ to $0.3$. The resulting island widths are shown in Figure 7, where it can be seen that the value of $\hat{C}$ does not significantly affect the agreement between SPEC, HMHD, and the analytical MRE. In addition to the island width, one can look at the structure of the island in SPEC and compare it with the island generated in HMHD. As shown in Figure 4, the island from SPEC matches that of HMHD closely. Its shape is dominantly given by the $m=1$ harmonic, but it seems to have a slightly larger content of $m>1$ harmonics. It is interesting to note that the two interfaces around the island, shown in green in Fig. 4b, are effectively coincident on the island separatrix, such that the entire volume between them is filled by the island. This is a consequence of the choice of $\Psi_{w}$ that maximizes the island width.

The saturated profiles in SPEC also agree well with the profiles from HMHD, as can be seen in Figure 8. The pressure in HMHD flattens across the island, which in this case happens primarily due to the generated parallel plasma flow, and not the parallel diffusion. In SPEC, the flattening is a consequence of the MRxMHD model. The radial profile of $B_{y}$ over the entire domain is not much affected by the tearing mode, meaning that the initial agreement between SPEC and HMHD is sustained in the saturated state.

Two important parameters in SPEC are $m_{\textrm{pol}}$ and $N_{\textrm{vol}}$. The first one determines the number of poloidal harmonics used for representing the interfaces and force imbalances, while the second one is the number of relaxed volumes in SPEC and it determines how well the outer ideal solution is represented. In Figure 9, the saturation width of a classical tearing mode obtained with SPEC is shown as a function of both $m_{\textrm{pol}}$ and $N_{\textrm{vol}}$. It can be seen that the saturated island width is converged, in the sense that increasing the two parameters beyond $m_{\textrm{pol}}=3$ and $N_{\textrm{vol}}=151$, which are used for the simulations presented in this paper, leads to less than $1\%$ difference in the calculated widths.

Figure 7: Scan of island width as a function of the bootstrap current strength for SPEC (green and blue) and HMHD (black), in the case of a classical tearing mode ($\Delta^{\prime}_{0}\>d=0.26$) and an NTM ($\Delta^{\prime}_{0}\>d=-0.26$).

Figure 8: Profiles of a) pressure and b) poloidal field for a saturated tearing mode from SPEC and HMHD. Flattened pressure in HMHD matches well with SPEC where the pressure is, ipso facto, flat inside the island. The poloidal field is mostly unchanged from the initial equilibrium.

Figure 9: Convergence of the width of a classical tearing mode with the bootstrap current in SPEC with respect to the number of a) poloidal harmonics and b) relaxed volumes. Reaching an equilibrium with a saturated island requires at least 3 harmonics.

5 Discussion

We have demonstrated that an "equilibrium approach", based on the solver SPEC (Hudson et al., 2012), can be used to directly predict the saturated state of tearing modes that evolve under the influence of the neoclassical bootstrap current. Good agreement with the resistive MHD equations, solved using the code HMHD (Huang & Bhattacharjee, 2016), was found for linearly unstable classical tearing modes, where the bootstrap current enhances the width of the saturated island. More importantly, we have shown that the approach can be applied to neoclassical tearing modes, which are linearly stable and for which bootstrap current is the primary driving mechanism. In addition to demonstrating the agreement in the width of saturated islands, a good match was also shown for the saturated island shapes and profiles of the pressure and the magnetic field. The independence of the SPEC results with respect to resolution parameters was established, as shown by the scans of the island width as a function of the number of poloidal harmonics and relaxed volumes. Having demonstrated that SPEC can correctly reproduce the saturated state of tearing modes, it is important to stress the relative benefits of this approach. First and foremost, the new method is significantly faster than resistive MHD simulations. For the scan of tearing modes with the bootstrap current in Figure 5, it took approximately 300 CPU hours to simulate with HMHD the saturation of just a single tearing mode. Conversely, it took a mere 2 CPU hours to simulate the entire dataset for SPEC, which includes 30 different saturated islands. An additional benefit is that the SPEC-based approach does not involve fitting coefficients. Such is the case with the analytical MRE model which, despite being fast, has several coefficients that need to be fit to a specific device ($\sim 10$ coefficients for the TCV tokamak) (Kong, 2020). The last notable benefit of the approach is that it is general, both in terms of geometry and particular instability. While initial work (Loizu et al., 2020) and its extension presented here were done in a slab, this approach could be applied to other configurations. The method has already been applied to classical tearing modes in a cylindrical tokamak (Loizu & Bonfiglio, 2023), and work has started in toroidal geometry. We are interested in predicting 2-1 NTMs observed in the TCV tokamak (Kong, 2020) and comparison with experimental measurements and MRE models. More broadly, this work promotes the development of a fast and simple tool for the general prediction of saturated instabilities in arbitrary geometries. In particular, this method may eventually allow for quick prediction of nonlinear saturation in stellarators, which have been observed (de Aguilera et al., 2015; Watanabe et al., 2005; Geiger J. E., 2004) to operate stably in regimes deemed unstable by linear theory.

6 Acknowledgments

This work was supported in part by the Swiss National Science Foundation. And by a grant from the Simons Foundation (1013657, JL). This work has been carried out within the framework of the EUROfusion Consortium, via the Euratom Research and Training Programme (Grant Agreement No 101052200 — EUROfusion) and funded by the Swiss State Secretariat for Education, Research and Innovation (SERI). 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 SERI. Neither the European Union nor the European Commission nor SERI can be held responsible for them. We thank Antoine Baillod for his help with the codes and insightful discussions.

Appendix A

Linear stability in SPEC has been verified against numerical calculations of $\Delta^{\prime}_{0}$ for classical tearing modes (Loizu et al., 2020). The stability for equilibria with bootstrap current is more involved. In Figure 10, we measure the size of a marginally stable slab as a function of the number of relaxed volumes. Each set of points in the plot represents marginal stability, such that points above are unstable and below stable. The top curve in blue represents equilibria without bootstrap current, and it can be seen that it has good convergence as $N_{\textrm{vol}}$ is increased, in this case to a marginal size of $\sim 3.0$. However, the same is not true for the equilibria with bootstrap current shown in red, for which the limit $N_{\textrm{vol}}\rightarrow\infty$ does not seem to exist. This leads us to believe that such equilibria are linearly unconditionally unstable in SPEC, where stability in any particular equilibria is due to limited resolution. Such a conclusion is supported by the bootstrap term from the MRE (see Figure 1). In SPEC, the ratio of parallel to perpendicular heat conductivity is infinite, which means that the term that accounts for heat transport in the island becomes $w_{0}=0$. This leads to the bootstrap contribution becoming $\Delta^{\prime}_{\textrm{bs}}\sim 1/w$, which also leads to unconditional linear stability in the MRE model.

Figure 10: Size of a marginally stable slab in SPEC as a function of the number of volumes. The top curve (blue) is for the case with $J_{\textrm{bs}}=0$, which converges to a well-defined value. The bottom curve (red) has $J_{\textrm{bs}}\neq 0$ and it does not have a limit.

Appendix B

Here we derive the MRE for the width of tearing modes in slab with bootstrap current. Neglecting flows which are assumed to be weak, we start from Ohm’s law

$$\frac{\partial\psi}{\partial t}=\eta(J_{z0}-J_{z})+\eta J_{\textrm{bs}}$$ (8)

Expressing the flux function as a sum of the equilibrium and perturbed part, $\psi(x,y,t)=\psi_{0}(x)+\psi_{1}(x,y,t)$, and assuming gradients in $x$ dominate over those in $y$, $\partial/\partial x\gg\partial/\partial y$,

$$\frac{\partial\psi_{1}}{\partial t}=\eta\frac{\partial^{2}\psi_{1}}{\partial x^{2}}+\eta J_{\textrm{bs}}$$ (9)

Assuming all quantities are of a single harmonicity in the poloidal $y$ direction, $k=2\pi/L$, the functions are expressed as $f(x,y,t)\sim\hat{f}(x,t)\cos(ky)$. Next, we integrate both sides from $-w/2$ to $w/2$, where $w$ is the island width at time $t$

$$\int_{-w/2}^{w/2}dx\>\frac{\partial\psi_{1}}{\partial t}=\int_{-w/2}^{w/2}dx\>\eta\frac{\partial^{2}\psi_{1}}{\partial x^{2}}+\int_{-w/2}^{w/2}dx\>\eta J_{\textrm{bs}}$$ (10)

For the term on the left hand side, we apply the constant-$\psi$ approximation to lift $\psi$ outside the integral, and for the bootstrap term, we assume that the pressure gradient, and hence $J_{\textrm{bs}}$, is slowly varying across the island. We obtain

$$\frac{\partial\psi_{1}}{\partial t}w=\eta\frac{\partial\psi_{1}}{\partial x}\Bigr{|}_{-w/2}^{w/2}+\eta J_{\textrm{bs}}w$$ (11)

To proceed further, we need an additional equation that relates the perturbed flux and instantaneous island width. Following an approach similar to Goldston (2020), the total flux is again composed of equilibrium and perturbed contributions. Expanding around the resonant surface, we use first-order expansion for the field $B_{y0}(x)=B^{\prime}_{y0}x$, which gives $\psi_{0}=-1/2\>B^{\prime}_{y0}x^{2}$, since we know that $B_{y0}(x)=-{\partial\psi_{0}}/{\partial x}$. For the perturbed flux, we use the zeroth-order expansion in $x$ and the fact that the island is of a single harmonicity $\psi_{1}(x,y,t)=-\hat{\psi}\cos(ky)$.

Evaluating the flux function at the X-point, we get $\psi(x=0,y=0)=-\hat{\psi}$, and the following expression traces the field line ($x_{l}$, $y_{l}$) corresponding to that flux

$$\hat{\psi}(1-\cos(ky_{l}))=1/2\>B^{\prime}_{y0}x^{2}_{l}$$ (12)

The island width is then obtained from the maximum excursion of the field line in $x$

$$w(t)=2\max(x_{l})=\left[\frac{16\hat{\psi}}{B^{\prime}_{y0}}\right]^{1/2}$$ (13)

Thus we obtain the relation between the reconnected flux and the island width

$$\psi=bw^{2}\>\>\textrm{ with }b=\frac{B^{\prime}_{y0}}{16}$$ (14)

Finally, this expression is inserted into the relation Eq. B4, leading to the MRE model

$$\frac{2}{\eta}\frac{\partial w}{\partial t}=\Delta^{\prime}+\frac{1}{b}J_{\textrm{bs}}\frac{1}{w}$$ (15)

$$\textrm{with }\Delta^{\prime}=\frac{1}{\psi_{1}}\frac{\partial\psi_{1}}{\partial x}\Bigr{|}_{-w/2}^{w/2}$$ (16)