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

E. Balkovic1  J. Loizu1  J. P. Graves1  Y.-M. Huang2  C. Smiet1  1 École
Polytechnique Fédérale de Lausanne, Swiss Plasma Center, CH-1015 Lausanne,
Switzerland
2 Princeton University, Princeton NJ 08543, USA

###### 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:

$\displaystyle\frac{\partial\rho}{\partial
t}+\bnabla\boldsymbol{\cdot}\left(\rho\boldsymbol{u}\right)=0$ (3)
$\displaystyle\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}}$
$\displaystyle\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) $\displaystyle\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)

$\displaystyle\bnabla\times\boldsymbol{B}=\mu_{l}\boldsymbol{B}$ (6)
$\displaystyle\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)

## References

  * de Aguilera et al. (2015) de Aguilera, Adriana M., Castejón, Francisco, Ascasíbar, Enrique, Blanco, Emilio, la Cal, Eduardo De, Hidalgo, Carlos, Liu, Bing, López-Fraguas, Antonio, Medina, Francisco, Ochando, María Antonia, Pastor, Ignacio, Ángeles Pedrosa, María, Milligen, Boudewijn Van, Velasco, José Luis & the TJ-II Team 2015 55 (11), 113014.
  * Baillod et al. (2021) Baillod, A., Loizu, J., Qu, Z.S., Kumar, A. & Graves, J.P. 2021 Computation of multi-region, relaxed magnetohydrodynamic equilibria with prescribed toroidal current profile. Journal of Plasma Physics 87 (4), 905870403.
  * Berger (1999) Berger, Mitchell A 1999 Introduction to magnetic helicity. Plasma Physics and Controlled Fusion 41 (12B), B167.
  * Brunetti et al. (2014) Brunetti, D., Graves, J.P., Cooper, W.A. & Terranova, D. 2014 Ideal saturated mhd helical structures in axisymmetric hybrid plasmas. Nuclear Fusion 54 (6), 064017.
  * Cooper et al. (2010) Cooper, W. A., Graves, J. P., Pochelon, A., Sauter, O. & Villard, L. 2010 Tokamak magnetohydrodynamic equilibrium states with axisymmetric boundary and a 3d helical core. Phys. Rev. Lett. 105, 035003.
  * Courant et al. (1967) Courant, R., Friedrichs, K. & Lewy, H. 1967 On the partial difference equations of mathematical physics 11 (2), 215–234.
  * Escande & Ottaviani (2004) Escande, D.F. & Ottaviani, M. 2004 Simple and rigorous solution for the nonlinear tearing mode. Physics Letters A 323 (3), 278–284.
  * Fitzpatrick (1995) Fitzpatrick, Richard 1995 Helical temperature perturbations associated with tearing modes in tokamak plasmas. Physics of Plasmas 2, 825–838.
  * Geiger J. E. (2004) Geiger J. E., Weller, A. Zarnstorff M. C. Nührenberg C. Werner A. H. F. Kolesnichenko Y. I. W7-AS Team and Neutral Beam Injection Group 2004 Equilibrium and stability of high-b plasmas in wendelstein 7-as. Fusion Science and Technology 46 (1), 13–23.
  * Glasser et al. (1976) Glasser, A. H., Greene, J. M. & Johnson, J. L. 1976 Resistive instabilities in a tokamak. The Physics of Fluids 19 (4), 567–574.
  * Goldston (2020) Goldston, R.J 2020 Introduction to Plasma Physics. CRC Press.
  * Günter et al. (1999) Günter, S, Gude, A, Maraschek, M, Yu, Q & the ASDEX Upgrade Team 1999 Influence of neoclassical tearing modes on energy confinement. Plasma Physics and Controlled Fusion 41 (6), 767.
  * Helander & Sigmar (2005) Helander, Per & Sigmar, D. J. 2005 Collisional transport in magnetized plasmas. Cambridge University Press.
  * Hirshman & Whitson (1983) Hirshman, S. P. & Whitson, J. C. 1983 Steepest-descent moment method for three-dimensional magnetohydrodynamic equilibria. The Physics of Fluids 26 (12), 3553–3568.
  * Hoelzl et al. (2021) Hoelzl, M., Huijsmans, G.T.A., Pamela, S.J.P., Bécoulet, M., Nardon, E., Artola, F.J., Nkonga, B., Atanasiu, C.V., Bandaru, V., Bhole, A., Bonfiglio, D., Cathey, A., Czarny, O., Dvornova, A., Fehér, T., Fil, A., Franck, E., Futatani, S., Gruca, M., Guillard, H., Haverkort, J.W., Holod, I., Hu, D., Kim, S.K., Korving, S.Q., Kos, L., Krebs, I., Kripner, L., Latu, G., Liu, F., Merkel, P., Meshcheriakov, D., Mitterauer, V., Mochalskyy, S., Morales, J.A., Nies, R., Nikulsin, N., Orain, F., Pratt, J., Ramasamy, R., Ramet, P., Reux, C., Särkimäki, K., Schwarz, N., Verma, P. Singh, Smith, S.F., Sommariva, C., Strumberger, E., van Vugt, D.C., Verbeek, M., Westerhof, E., Wieschollek, F. & Zielinski, J. 2021 The jorek non-linear extended mhd code and applications to large-scale instabilities and their control in magnetically confined fusion plasmas. Nuclear Fusion 61 (6), 065001.
  * Hole et al. (2006) Hole, M. J., Hudson, S. R. & Dewar, R. L. 2006 Stepped pressure profile equilibria in cylindrical plasmas via partial taylor relaxation. Journal of Plasma Physics 72 (6), 1167–1171.
  * Huang & Bhattacharjee (2016) Huang, Yi-Min & Bhattacharjee, A. 2016 Turbulent magnetohydrodynamic reconnection mediated by the plasmoid instability. The Astrophysical Journal 818 (1), 20.
  * Hudson et al. (2012) Hudson, S. R., Dewar, R. L., Dennis, G., Hole, M. J., McGann, M., von Nessi, G. & Lazerson, S. 2012 Computation of multi-region relaxed magnetohydrodynamic equilibria. Physics of Plasmas 19.
  * Kong (2020) Kong, Mengdi 2020 Towards integrated control of tokamak plasmas: physics-based control of neoclassical tearing modes in the tcv tokamak p. 228.
  * Kong et al. (2022) Kong, M, Felici, F, Sauter, O, Galperti, C, Vu, T, Ham, C J, Hender, T C, Maraschek, M, Reich, M, the TCV team, the ASDEX Upgrade team, the MAST team & the EUROfusion MST1 Team 2022 Physics-based control of neoclassical tearing modes on tcv. Plasma Physics and Controlled Fusion 64 (4), 044008.
  * Kumar (2022) Kumar, A. 2022 Plasma Physics and Controlled Fusion 64 (6), 065001.
  * Loizu & Bonfiglio (2023) Loizu, J. & Bonfiglio, D. 2023 Nonlinear saturation of resistive tearing modes in a cylindrical tokamak with and without solving the dynamics. Journal of Plasma Physics 89 (5), 905890507.
  * Loizu et al. (2020) Loizu, J., Huang, Y. M., Hudson, S. R., Baillod, A., Kumar, A. & Qu, Z. S. 2020 Direct prediction of nonlinear tearing mode saturation using a variational principle. Physics of Plasmas 27.
  * Loureiro et al. (2005) Loureiro, N. F., Cowley, S. C., Dorland, W. D., Haines, M. G. & Schekochihin, A. A. 2005 $x$-point collapse and saturation in the nonlinear tearing mode reconnection. Phys. Rev. Lett. 95, 235003.
  * Militello & Porcelli (2004) Militello, Fulvio & Porcelli, F. 2004 Simple analysis of the nonlinear saturation of the tearing mode. Physics of Plasmas 11, L13.
  * Muraglia et al. (2021) Muraglia, M., Poyé, A., Agullo, O., Dubuit, N. & Garbet, X. 2021 Nonlinear dynamics of ntm seeding by turbulence. Plasma Physics and Controlled Fusion 63.
  * Rutherford (1973) Rutherford, P. H. 1973 Nonlinear growth of the tearing mode. The Physics of Fluids 16 (11), 1903–1908.
  * Sauter et al. (2002) Sauter, O, Buttery, R J, Felton, R, Hender, T C, Howell, D F & contributors to the EFDA-JET Workprogramme 2002 Marginal beta-limit for neoclassical tearing modes in jet h-mode discharges. Plasma Physics and Controlled Fusion 44 (9), 1999.
  * Taylor (1974) Taylor, J. B. 1974 Relaxation of toroidal plasma and generation of reverse magnetic fields. Phys. Rev. Lett. 33, 1139–1141.
  * Teng et al. (2018) Teng, Q., Ferraro, N., Gates, D.A. & White, R.B. 2018 Nonlinear simulations of thermo-resistive tearing mode formalism of the density limit. Nuclear Fusion 58 (10), 106024.
  * Watanabe et al. (2005) Watanabe, K.Y., Sakakibara, S., Narushima, Y., Funaba, H., Narihara, K., Tanaka, K., Yamaguchi, T., Toi, K., Ohdachi, S., Kaneko, O., Yamada, H., Suzuki, Y., Cooper, W.A., Murakami, S., Nakajima, N., Yamada, I., Kawahata, K., Tokuzawa, T., Komori, A. & experimental group, LHD 2005 Effects of global mhd instability on operational high beta-regime in lhd. Nuclear Fusion 45 (11), 1247.