Nonlinear saturation of resistive tearing modes in a cylindrical tokamak with and without solving the dynamics

Abstract

We show that the saturation of resistive tearing modes in a cylindrical tokamak, as well as the corresponding island width, can be directly calculated with a magnetohydrodynamics (MHD) equilibrium code without solving the dynamics and without considering resistivity. The results are compared with initial-value resistive MHD simulations and to an analytical nonlinear theory. For small enough islands, the agreement is remarkable. For sufficiently large islands, the equilibrium calculations, which assume a flat current profile inside the island, overestimate the saturation amplitude. On the other hand, excellent agreement between nonlinear resistive MHD simulations and nonlinear theory is observed for all the considered tearing unstable equilibria.

1 Introduction

In magnetic fusion, tearing modes are an important class of instabilities that can often be described by resistive magnetohydrodynamics (MHD) (White 2014). Tearing modes are current-driven internal instabilities that drive spontaneous magnetic reconnection at rational surfaces in toroidally confined plasmas. Magnetic islands grow in time, thus potentially affecting confinement, and eventually saturate nonlinearly. If their saturated amplitude is too large, they can trigger disruptions, sometimes even inducing the complete termination of the plasma discharge (La Haye 2006). While resistivity is required to allow these modes to grow, their saturated states are three-dimensional MHD equilibria that are often independent of resistivity, at least for sufficiently large Lundquist numbers (Poyé et al. 2014). This suggests that a direct MHD equilibrium computation might be able to describe and even predict these saturated states, as was shown to be the case for the saturation of ideal instabilities in tokamaks (Cooper et al. 2010, 2011; Kleiner et al. 2019) or the formation of saturated double-axis states in reversed field pinches (Dennis et al. 2013b). The possibility of finding these nonlinearly saturated tearing modes directly, namely, without solving the dynamics and without free parameters, was investigated and successfully shown for the first time in a slab geometry (Loizu et al. 2020). More recently, a similar type of calculation was carried out in a cylindrical geometry for the case of externally forced reconnection in a cylinder (Wright et al. 2022), although the volume available for reconnection remained unconstrained. Here, we investigate the possibility of directly predicting the saturated state of tearing modes in a cylindrical tokamak geometry at zero pressure. In this configuration, we can exploit an exact nonlinear saturation theory (Arcis, Escande & Ottaviani 2006) that predicts the saturated island width for a given initially axisymmetric equilibrium that is unstable to a tearing mode. The results are compared with those obtained with the SpeCyl code (Cappello & Biskamp 1996; Bonfiglio, Chacón & Cappello 2010), which solves the time-dependent nonlinear visco-resistive MHD equations in a cylinder. These predictions serve as a benchmark for the direct equilibrium approach, where we use the stepped-pressure equilibrium code (SPEC) (Hudson et al. 2012, 2020) to describe, without free parameters, the spontaneous formation of these nonlinearly saturated tearing islands. In general, the direct calculation of saturated states with an equilibrium approach is interesting for two reasons: on the one hand, it provides insights into the nature of the saturated state itself; on the other hand, the approach is potentially much faster than the initial-value approach, especially as we consider larger systems (such as magnetic fusion reactors) or more complex geometries (such as stellarators), where analytical theories are not available and initial-value approaches are very demanding.

The remainder of this paper is organized as follows. In § 2 we construct and analyze the linear stability of a family of MHD equilibria both analytically and with SPEC. Section 3 describes a nonlinear theory that predicts the saturated island width of these unstable equilibria. In § 4 we carry out nonlinear simulations with the SpeCyl code and compare the obtained saturated tearing modes with the theoretical predictions. Section 5 describes how SPEC can also be used to retrieve the saturated states corresponding to each unstable equilibrium, and a discussion on the level of agreement is provided. Conclusions follow in § 6.

2 Family of equilibria and linear stability analysis

We construct a family of non-rotating, cylindrical MHD equilibria with zero pressure and parametrized by the on-axis safety factor $q_0$

(2.1)\begin{equation} q(r) = q_0\left[1+\left(\frac{r}{r_0}\right)^2\right], \end{equation}

for $0\leq r\leq a$ and where a and $r_0$ are fixed and correspond, respectively, to the radius of the plasma boundary and to the current channel width. As we will see, the value of $q_0$ controls the stability of the equilibrium, the presence and location of resonant rational surfaces as well as the width of the saturated islands. In the strong guide field limit ($B_z\gg B_{\theta}$ where $\boldsymbol{B}$ is the magnetic field), which for a tokamak with $q\sim 1$ also corresponds to the large aspect ratio limit ($R \gg a$ where R is the tokamak major radius), the axial current density profile associated with this $q$-profile is given by

(2.2)\begin{equation} j_z(r) = \frac{j_0}{\left[1+\left(\dfrac{r}{r_0}\right)^2\right]^2}, \end{equation}

where $j_0 = 2B_0/(\mu_0q_0R)$ is the on-axis current density and $\mu_0$ the vacuum magnetic permeability, $B_0=B_z(0)$ and $2{\rm \pi} R$ is the length of the cylinder, which is considered to be $2{\rm \pi}$-periodic in both the polar angle $\theta$ and $\varphi =z/R$. Figure 1 shows examples of profiles corresponding to (2.1)–(2.2). Given the fixed, circular plasma boundary at $r=a$, and given $q_0$, $r_0$, $B_0$ and $R$, the ideal MHD equilibrium magnetic field $\boldsymbol {B}$ is fully determined by the force-balance equation $\boldsymbol {j}\times \boldsymbol {B}=\boldsymbol {\nabla } p=0$, which reduces to a one-dimensional, first-order differential equation for $B_z(r)$. These equilibria are ideally stable as long as $q_0>1$ (Freidberg 2014), but as we will see in the next section, they can be resistive unstable (Furth, Rutherford & Selberg 1973).

Figure 1. Equilibrium profiles for the safety factor (a) and current density (b) as given by (2.1)–(2.2) (solid black lines) and as obtained from SPEC (red circles and lines) with $N_v=41$. Here, $q_0=1.2$ and $r_0/a=0.81$.

We can also construct these equilibria numerically by using SPEC (Hudson et al. 2012, 2020), with the purpose of seeking nearby, relaxed (lower potential energy) equilibria, which can be accessed via the possible opening of magnetic islands. SPEC considers a finite number $N_v$ of nested volumes, hereafter referred to as relaxation volumes, that are separated by magnetic surfaces, hereafter referred to as ideal interfaces. In a cylindrical geometry, the position of each ideal interface $l$ is generally described by $\boldsymbol{x}_l(\theta,\varphi) = r_l(\theta,\varphi)\cos{\theta} \hat{i} + r_l(\theta,\varphi)\sin{\theta} \hat{j} + R\varphi \hat{k}$, with $(\hat{i},\hat{j},\hat{k})$ the cartesian unit vectors, and thus can be determined by a set of Fourier coefficients, since we can write $r_l(\theta,\varphi )=\sum _{m,n} r_{l,mn}\cos {(m\theta -n\varphi )}$. In each relaxation volume, the magnetic field solves a Beltrami equation $\boldsymbol {\nabla }\times \boldsymbol {B}=\mu \boldsymbol {B}$ with $\mu$ a constant that is related to the parallel current density. The solution in the volumes is then determined by the shape of the boundary interfaces where $\boldsymbol {B}\boldsymbol {\cdot }\boldsymbol {n}=0$, the value of $\mu$ and two scalars $\varPsi _p$ and $\varPsi _t$ that correspond to the ‘poloidal’ (constant-$\theta$) and ‘toroidal’ (constant-$\varphi$) magnetic fluxes enclosed by each volume (in the innermost volume, only $\varPsi _t$ is required). The position and shape of each interface is unknown a priori and is determined by a jump condition across each interface, $[[B^2]]=B^2_+ - B^2_- = 0$, where $+$ and $-$ refers to each side of the interface. This condition is the local equivalent of force balance and ensures that the inner and outer pressures are in balance everywhere on an interface. We remark that equilibria with finite plasma pressure can also be described in SPEC but here we are considering zero-pressure equilibria. The magnetic fields satisfying these equations are extrema of the plasma potential energy, $W=\int _{V_p} B^2 / (2\mu _0)\,{\rm d}V$, when considering arbitrary variations in the magnetic field inside the relaxation volumes and ideal variations of the interface geometries, with the constraint of a conserved magnetic helicity, $K_l=\int _{V_l} \boldsymbol {A}\boldsymbol {\cdot }\boldsymbol {B} \, {\rm d}V$, and magnetic fluxes in each relaxation volume (Hudson et al. 2012). Here, $V_p=\sum _{l=1}^{N_v} V_l$ is the total plasma volume and $\boldsymbol {A}$ is the magnetic vector potential. Furthermore, it has been shown that the magnetic field satisfying these equations retrieves exactly ideal MHD in the limit $N_v\rightarrow \infty$ (Dennis et al. 2013a). Hence, SPEC can be used to construct the family of ideal MHD equilibria described by (2.1) if a large value of $N_v$ is chosen and if the appropriate constraints are provided. In fact, the constants ($\mu$, $\varPsi _p$, $\varPsi _t$) in each volume, which determine the Beltrami fields, can be related to other physical quantities such as the safety factor $q$ on each side of each interface, $q_+$ and $q_-$, or the net toroidal current flowing in each volume, $I_{\mathrm {vol}}$, and on each interface, $I_{\mathrm {surf}}$ (Baillod et al. 2021). One can thus for example constrain in each volume $(q_+, q_-, \varPsi _t)$ if the safety factor is to be imposed, or $(I_{\mathrm {vol}} , I_{\mathrm {surf}}, \varPsi _t)$ if the current profile is to be imposed. Figure 1 shows an example of an equilibrium reconstructed with SPEC, where $N_v=41$ and the safety factor was constrained on each interface by following (2.1). After the SPEC calculation, one can look at the global behaviour of $q(r)$ as well as that of $j_z(r)$, which indeed reproduce well the analytical profiles.

2.1 Linear stability from classical tearing theory

The linear resistive stability of a cylindrical tokamak with zero pressure is uniquely determined by the sign of $\varDelta '$ (Furth et al. 1973)

(2.3)\begin{equation} \varDelta' = \lim_{\epsilon\rightarrow 0^+}\frac{\hat{\psi}'(r_s+\epsilon)-\hat{\psi}'(r_s-\epsilon)}{\hat{\psi}(r_s)}, \end{equation}

with $\varDelta '>0$ implying instability. Here, $\hat {\psi }(r)$ is the amplitude of the perturbed flux function $\psi (r,\theta,\varphi )=\hat {\psi }(r)\cos {(m\theta -n\varphi )}$ for the linear ideal mode $(m,n)$ associated with the rational surface $q(r_s)=q_s=m/n$, where $m$ and $n$ are, respectively, the poloidal (along $\theta$) and toroidal (along $\varphi$) mode numbers. The corresponding magnetic field perturbation is $\delta \boldsymbol {B}=\boldsymbol {\nabla }\psi \times \boldsymbol {\hat {z}}$, and thus the radial magnetic field perturbation, $\delta B_r$, for this mode is proportional to $\hat {\psi }$. Notice that $\varDelta '$ depends on equilibrium profiles but is independent of resistivity. It can be calculated from the linearized ideal MHD equations, by matching the so-called ‘outer solutions’. A shooting code provides a numerical solution for $\hat {\psi }(r)$ and a value of $\varDelta '$ for a given resonant surface. For the family of equilibria described by (2.1), choosing $r_0/a=0.81$, $R/a=10$ and assuming the presence of a perfectly conducting shell at the plasma boundary, we find instability for the $m=2$, $n=1$ mode whenever $1< q_0<2$. In that range of $q_0$ values, we also find no other instability, hence, we can focus on the $q_s=2$ resonant surface. Figure 2 shows an example of solution for $\hat {\psi }(r)$ for $q_0=1.2$. Figure 3 shows the obtained values of $\varDelta '$ as a function of $q_0$. For larger values of $q_0$, the value of $\varDelta '$ keeps increasing but of course for $q_0>2$ there is no instability because there is no $q_s=2$ surface anymore.

Figure 2. Radial profile of the perturbed flux function $\hat {\psi }$ for the unstable $m=2,n=1$ mode obtained from the linearized ideal MHD equations (black solid line). The same function can be inferred from the linear radial displacement $\xi (r)$ obtained with SPEC using the relation $\hat {\psi }(r)=rB_{z0}(1/q(r)-1/q_s)\xi (r)/R$. The initial equilibrium is the same as for figure 1. The dashed line indicates $r=r_s$.

Figure 3. Value of $\varDelta '$ vs $q_0$ as obtained from a shooting code (black solid line) for the $m=2,n=1$ mode. The equilibrium $q$ profile is given in (2.1) with $r_0/a=0.81$ and $R/a=10$. The smallest eigenvalue of the force-gradient matrix from SPEC, $\lambda _{{\rm SPEC}}$, is also shown in arbitrary units for each value of $q_0$ (red circles).

2.2 Linear stability from SPEC

The stability of each of these equilibria can also be assessed with the SPEC code by evaluating the eigenvalues of the Hessian matrix $\mathcal {H}$, whose coefficients measure the second variation of the plasma potential energy with respect to perturbations in the geometry of the ideal interfaces, with the constraint of fixed magnetic helicity and fluxes in each volume (Loizu & Hudson 2019; Kumar et al. 2021, 2022). In a cylindrical geometry, however, one can also directly evaluate the eigenvalues of the force-gradient matrix $\mathcal {G}$, whose coefficients measure the variation in the force on each interface, $f_l = [[B^2/2\mu _0]]_l$, with respect to perturbations in the geometry of the ideal interfaces, $r_{l,mn}$, with the constraint of fixed magnetic helicity and fluxes in each volume (see Appendix A for more details). Let us call $\lambda$ the smallest eigenvalue of $\mathcal {G}$. If $\lambda \geq 0$, the equilibrium is stable. If $\lambda <0$, then the equilibrium is unstable and the eigenmode $\boldsymbol {u}_{\lambda }\equiv \{r_{l,mn}\}$ corresponding to the eigenvalue $\lambda$ provides information about the mode number $(m,n)$ of the instability as well as its radial structure. For the family of equilibria described by (2.1), we find that, for $1< q_0<2$, there is a negative eigenvalue (thus instability) that corresponds to the $m=2$, $n=1$ mode (figure 3). Outside this interval of $q_0$-values, all eigenvalues are positive (thus stability). In the case of instability, the radial structure of the unstable eigenmode agrees well with the one obtained from the shooting code, as shown in figure 2. Indeed, this comparison is carried out by relating the radial displacement of the flux surfaces to the perturbed flux function, $\delta \boldsymbol {B} = \boldsymbol {\nabla } \times ( \xi \times \boldsymbol {B}_0) = \boldsymbol {\nabla } \psi \times \boldsymbol {\hat {z}}$, and evaluating the radial displacement of the SPEC ideal interfaces, which has an amplitude $\xi (r_l)=r_{l,mn}-r_{l,00}$ for each mode $(m,n)$.

3 Nonlinear saturation theory

The growth of the tearing mode promotes the growth of a magnetic island. The rate $\gamma$ at which the island grows depends on the plasma resistivity $\eta$, both in the linear (Furth et al. 1973) and nonlinear (Rutherford 1973) phases. The saturation of the tearing mode, which in itself proves the existence of a nearby, generally three-dimensional MHD equilibrium, can be characterized by the width of the saturated island, $w_{{\rm sat}}$, and that value becomes independent of resistivity for sufficiently small resistivity (Poyé et al. 2014). For a cylindrical tokamak, Arcis et al. (2006) derived a nonlinear theory for the time evolution of the island width $w(t)$

(3.1)\begin{equation} \frac{\mu_0}{\eta_{{\rm eq}}(r_s)}\frac{{\rm d}w}{{\rm d}t} = 1.22\varDelta' + w\left\{ \frac{\mathcal{A}^2}{2}\ln{\frac{w}{w_0} - 2.21\mathcal{A}^2 + 0.40\frac{\mathcal{A}}{r_s} + \frac{\mathcal{B}}{2}} + 0.17\sigma\frac{\mathcal{A}^2s}{2-s}\right\}, \end{equation}

that is exact to first order in $w$. The coefficients in (3.1) can all be calculated from the equilibrium profiles, namely,

(3.2a–c)\begin{equation} \mathcal{A} = \frac{j'_{{\rm eq}}(r_s)}{j_{{\rm eq}}(r_s)}\left(1-\frac{2}{s}\right) ;\quad \mathcal{B} = \frac{j''_{{\rm eq}}(r_s)}{j_{{\rm eq}}(r_s)}\left(1-\frac{2}{s}\right); \quad w_0 = a \exp{\left(-\frac{\varSigma'}{2\mathcal{A}}\right)}, \end{equation}

where $j_{{\rm eq}}(r)$ is the equilibrium current density profile, which in our case is given by (2.2), $s = r_sq'(r_s)/q(r_s)$ is the magnetic shear at the resonant surface,

(3.3)\begin{equation} \varSigma' = \lim_{\epsilon\rightarrow 0^+}\left[ \frac{\hat{\psi}'(r_s+\epsilon)+\hat{\psi}'(r_s-\epsilon)}{\hat{\psi}(r_s)} - 2\mathcal{A}\left(1+\ln{\left(\frac{\epsilon}{a}\right)}\right) \right] \end{equation}

is a quantity that, similar to $\varDelta '$, can be obtained from the ideal outer solution and $\sigma$ is a constant that defines whether the equilibrium resistivity is uniform ($\sigma =0$) or the equilibrium electric field is uniform ($\sigma =1$). Equation (3.1) can be used to predict the saturated island width, $w_{{\rm sat}}$, as a function only of the initial equilibrium profiles, by setting the left-hand side to zero and solving for $w$ (notice that, as expected, the solution will not depend on $\eta$). This involves solving a single nonlinear algebraic equation and that can be easily done numerically. Figure 4 shows the predicted values of $w_{{\rm sat}}$ for the same family of equilibria as in § 2, namely for different values of $q_0$. We observe that $w_{{\rm sat}}$ increases with $q_0$ up to a maximum, $w_{{\rm sat}}/a \approx 0.15$ at $q_0\approx 1.3$, and then decreases again towards zero as $q_0$ approaches $q_s=2$. It is worth noticing that the value of $\varDelta '$, which gives a measure of the available energy for reconnection (Furth et al. 1973), is monotonically increasing with $q_0$ and thus the final amplitude of the mode is not entirely reflected by the amplitude of $\varDelta '$. We also observe that the choice of $\sigma$ (0 or 1) only introduces minor changes.

Figure 4. Value of $w_{{\rm sat}}$ vs $q_0$ as obtained from (3.1) with $\sigma =1$ (solid black line) and $\sigma =0$ (dashed grey line). Also shown are the predictions from the initial-value code SpeCyl (blue stars) and from the equilibrium code SPEC (red circles). The initially unstable equilibria are the same as for figure 3.

4 Nonlinear resistive MHD simulations

The dynamics of the tearing mode, from the initially unstable equilibrium up until the nonlinear saturation of the island, can be calculated by using an initial-value approach. Here, we use the SpeCyl code (Cappello 2004; Bonfiglio et al. 2010) to solve the nonlinear visco-resistive MHD equations in a cylinder, in the constant-pressure, constant-density approximation. The model equations are, in dimensionless units,

(4.1)\begin{equation} \left. \begin{aligned} \frac{\partial \boldsymbol{B}}{\partial t} & ={-}\boldsymbol{\nabla}\times\boldsymbol{E}, \\ \boldsymbol{E} & = \eta \boldsymbol{j} - \boldsymbol{v}\times\boldsymbol{B},\\ \boldsymbol{j} & = \boldsymbol{\nabla}\times\boldsymbol{B} ,\\ \rho\frac{{\rm d} \boldsymbol{v}}{{\rm d} t} & = \boldsymbol{j}\times\boldsymbol{B} + \rho\nu{\nabla}^2\boldsymbol{v} ,\\ \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{B} & =0, \end{aligned} \right\} \end{equation}

where $\boldsymbol {E}$ and $\boldsymbol {B}$ are the electric and magnetic fields, $\boldsymbol {j}$ is the current density, $\boldsymbol {v}$ is the plasma velocity, ${\textrm {d} \boldsymbol {v}}/{\textrm {d} t}=[{\partial \boldsymbol {v}}/{\partial t} + (\boldsymbol {v}\boldsymbol {\cdot }\boldsymbol {\nabla })\boldsymbol {v} ]$ and $\rho =m_in_0$ is the (constant) mass density. Lengths are normalized to the minor radius $a$, the plasma particle density to the uniform density $n_0$, the magnetic field to the initial value $B_0$ of the axial magnetic field on axis, the velocity to the Alfvén velocity $v_A=B_0/(\mu _0m_i n_0)^{1/2}$ and time to the Alfvén time $\tau _A=a/v_A$. In these units, the resistivity $\eta$ corresponds to the inverse Lundquist number $S^{-1}=\tau _A/\tau _{\eta }$, and the viscosity $\nu$ to the inverse viscous Lundquist number (for a scalar kinematic viscosity), $M^{-1}=\tau _A/\tau _{\nu }$, where $\tau _{\eta }$ and $\tau _{\nu }$ are the resistive time scale and viscous time scale, respectively. The equations are solved in cylindrical geometry $(r,\theta,z)$ and periodic boundary conditions are used for both the azimuthal coordinate $\theta$ and the axial coordinate $z$, which has periodicity $2{\rm \pi} R$. The code uses finite differences in the radial coordinate $r$ and a spectral formulation in the two periodic coordinates $\theta$ and $z$. In the calculations presented in this paper, the Fourier harmonics $(m,n)=(2j,j)$ are retained, with $j=0,\ldots,16$. Such large number of helical harmonics is required to properly capture the sharp flow patterns developing around the island separatrix, in particular for high Lundquist numbers. The boundary conditions at the plasma wall, $r=a$, are $B_r=E_{\theta }=0$ (perfect conductor), $E_z=E_0$ (externally imposed, constant, toroidal electric field) and $v_{\theta }=v_z=0$ (no-slip boundary conditions). The rest of the boundary conditions can be derived from Ohm's law. More details about the SpeCyl code and its numerical algorithm can be found in Cappello & Biskamp (1996) and Cappello (2004).

We run SpeCyl for a number of initially unstable equilibria that, as described in § 2, are characterized by $q_0$. For each value of $q_0$, we consider different values for the Lundquist number from $S=10^6$ to $S=10^9$ at fixed magnetic Prandtl number $P=S/M=1$, we evolve the system until saturation and measure the island width. Figure 5 illustrates examples of such saturated states by showing a Poincaré section of the magnetic field obtained for $q_0=1.3$ (top left) and $q_0=1.6$ (bottom left) with $S=10^8$. For sufficiently high Lundquist number $S \geq 10^8$, the radial profiles of the perturbed flux function $\hat {\psi }$ show good agreement in the linear growth phase with those produced by the shooting code, and the nonlinear saturated island width becomes independent of $S$ (data not shown). This is consistent with previous analytical and numerical studies (see e.g. Poyé et al. 2014) showing that the saturated tearing mode island width is independent of resistivity and viscosity. We also find very good agreement with the analytical predictions from the nonlinear theory of Arcis for all the considered tearing unstable equilibria (figure 4). In particular, the best agreement with nonlinear theory is found assuming $\sigma =1$ (i.e. a uniform equilibrium electric field) consistently with the assumption of the resistive MHD code. From the resistive simulations, we also observe that the saturated island is not symmetric with respect to the X-point. More precisely, we can define the asymmetry parameter

(4.2)\begin{equation} A_{{\rm sym}} = \frac{r_X-r_-}{r_+{-}r_X} - 1, \end{equation}

where $r_X$ is the radial position of the X-point of the island and $r_-$, $r_+$ are respectively the smallest and largest radial positions of the island separatrix (which occur at the poloidal position of the O-point). For a symmetric island, $r_X-r_- = r_+ - r_X$ and thus $A_{\textrm {sym}}=0$. We can evaluate (4.2) from the SpeCyl simulation results and display $A_{\textrm {sym}}$ as a function of $q_0$ (figure 6). Figure 7 also shows how the individual values of $r_X$, $r_-$ and $r_+$ change with $q_0$. As soon as $q_0>1$, the islands become quickly asymmetrical, $A_{\textrm {sym}}>0$, and become more and more asymmetric as $q_0$ approaches 2. We remark that the asymmetric feature of tearing mode islands has been observed in previous investigations and has been suggested to play an important role in determining the saturation amplitude from a modified Rutherford equation approach (Teng et al. 2016).

Figure 5. Poincaré sections of the magnetic field at $\varphi =0$ showing the saturated island as obtained from SpeCyl (a,c) and SPEC (b,d) for $q_0=1.3$ (a,b) and $q_0=1.6$ (c,d). The violet lines on the SPEC Poincaré sections (b,d) indicate the two interfaces encapsulating the resonant volume.

Figure 6. Value of $A_{\textrm {sym}}$ as a function of $q_0$ as obtained from SpeCyl (blue stars) and SPEC (red circles). The dashed line indicates the expected island asymmetry $A_{\textrm {max}}$ obtained by solving by solving (5.4) and that is used to initialize the SPEC equilibrium.

Figure 7. Radial positions $r_X$, $r_-$ and $r_+$ characterizing the island, as a function of $q_0$ as obtained from SpeCyl (blue) and SPEC (red). The dashed line indicates the position of the resonant surface.

5 Nonlinear equilibrium calculations

We now use the SPEC code to seek stable equilibrium solutions by starting from the unstable equilibria described in § 2. Indeed, in § 2.2 we obtained, for each value of $q_0$, the unstable eigenmode $\boldsymbol {u}_{\lambda }\equiv \{r_{l,mn}\}$ that is associated with the unstable (negative) eigenvalue $\lambda$ of the force-gradient matrix $\mathcal {G}$. Here, $m=2$, $n=1$ and the index $l=1,2,\ldots,N_v-1$ labels the interfaces. Figure 8 shows an example of such an eigenmode. We now perturb the Fourier geometry of the initial equilibrium interfaces by adding a displacement proportional to the eigenmode, namely, we add $\alpha \ \boldsymbol {u}_{\lambda }$ with $\alpha >0$. Notice that, given a certain volume (or equivalently, magnetic flux) enclosed by the two interfaces surrounding the resonant surface, there is a maximum perturbation amplitude beyond which the interfaces would overlap – but this is never exceeded because the perturbation is typically small enough. SPEC then seeks a new equilibrium solution. If $\alpha$ is too small, SPEC snaps back to the initial equilibrium (figure 9). If $\alpha$ is large enough, however, SPEC finds another equilibrium solution (independent of $\alpha$, see figure 9) that has a helical structure and displays a magnetic island inside the volume containing the resonant surface and thus where reconnection happens (two examples are shown in figure 5). This approach was successful in slab geometry and showed that the saturation amplitude of the tearing mode is correctly retrieved (Loizu et al. 2020). In that study, the toroidal magnetic flux $\varPsi _w$ enclosed by the resonant volume, which constrains the maximum achievable island width, was chosen by following a recipe without free parameters that is summarized in § 5.1 for convenience. If the island is expected to be perfectly symmetric around the resonant surface, the value of $\varPsi _w$ uniquely determines the initial positions $r_+$ and $r_-$ of the two interfaces defining the resonant volume, since $\varPsi _w = B_z {\rm \pi}(r_+^2 - r_-^2)$ and the symmetry argument forces the flux to be equally distributed on each side of the resonant radius $r_s$. In general, however, there are multiple ways of placing the two interfaces given the same total flux $\varPsi _w$ (see figure 10 for an illustration). For this reason, we need an additional constraint that is related to the expected island asymmetry, and that is the subject of § 5.2. Finally, the search for a nearby, lower energy equilibrium state can be done with different constraints, and an appropriate choice must be made by identifying the good quasi-invariants in the process at play. In the case of tearing mode reconnection, it is the net-toroidal-current profile that shall be constrained, and that is the subject of § 5.3. Section 5.4 summarizes the results of nonlinear simulations carried out with SPEC to reproduce the saturated tearing mode island width as a function of $q_0$.

Figure 8. Unstable eigenmode profile showing the typical tearing parity around the resonant surface. Here, $q_0=1.1$ and $N_v=41$.

Figure 9. Saturated island width obtained from SPEC nonlinear calculations with $q_0=1.1$ as a function of the initial perturbation amplitude as measured by $\alpha$. The normalizing factor $\alpha _{\textrm {sat}}$ represents the actual saturation amplitude of the helical state obtained after the nonlinear calculation. The Fourier resolution is $(m,n)=(4,2)$.

Figure 10. An illustration of how the initial positions $r_+$ and $r_-$ of the interfaces defining the resonant volume can be symmetric (a) or asymmetric (b) with respect to the resonant surface $r=r_s$, while preserving the same enclosed toroidal flux $\varPsi _w$.

5.1 Determining $\varPsi _w$

In the limit $\varPsi _w\rightarrow 0$, the volume available for reconnection vanishes and therefore the maximum island width that is possibly achievable goes to zero. For increasingly large values of $\varPsi _w$, the island is allowed to grow over a larger volume, although it may well only occupy a fraction of it. However, if the value of $\varPsi _w$ exceeds a certain threshold, $\varPsi _w^*$, the initial equilibrium becomes stable and no island grows at all. This is illustrated in figure 11 for two cases. Obviously, for each $q_0$, the saturated island width obtained with SPEC depends on $\varPsi _w$. However, there is a distinct point, which we identify as $w_{\textrm {sat}}$, and that is the maximum achieved island width, which occurs at around $\varPsi _w\approx \varPsi _w^*$. We thus apply the following procedure: for each value of $q_0$, a linear stability analysis is first performed with SPEC (see § 2.2) and the value of $\varPsi _w^*$ is obtained by solving $\lambda (\varPsi _w)=0$. This criterion is reminiscent of the quasi-linear theory of R. White (White et al. 1977), in which the saturation of a tearing mode occurs at a marginal stability point evaluated as $\varDelta '(w)=0$. We remark that the value of $\varPsi _w^*$ obtained with this procedure is independent of the number of volumes for sufficiently large $N_v$ (see figure 12). Then, SPEC is run nonlinearly with $\varPsi _w \lesssim \varPsi _w^*$ and the obtained value of $w_{\textrm {sat}}$ is extracted from the corresponding saturated state.

Figure 11. Island width, $w$, obtained from SPEC nonlinear calculations as a function of the enclosed toroidal flux $\varPsi _w$ imposed in the resonant volume. Two values of $q_0$ are shown: $q_0=1.1$ (orange squares) and $q_0=1.3$ (purple diamonds). The vertical dashed lines indicate the value of $\varPsi _w^*$. Here, the flux is always equally distributed on each side of the resonant surface, and the toroidal-current profile is constrained. The Fourier resolution is low here ($m=2,n=1$).

Figure 12. Values of $\varPsi _w^*$ (obtained by solving $\lambda (\varPsi _w)=0$) as a function of the number of volumes. Here, the equilibrium has $q_0=1.1$.

5.2 Model for the island asymmetry

The SPEC calculations shown in figure 11 were done by assuming that the flux $\varPsi _w$ in the resonant volume is equally distributed on each side of the resonant surface. However, a different choice (as illustrated in figure 10) can substantially modify the nonlinear saturation amplitude and thus the value of $w_{\textrm {sat}}$. In the following, we describe how we can estimate the expected asymmetry of the island and thus constrain accordingly the initial position of the interfaces, $r_+$ and $r_-$, around the resonant surface.

We start by writing the total magnetic field as $\boldsymbol {B} = \boldsymbol {B}_0(r) + \boldsymbol {\nabla }\psi \times \hat {\boldsymbol {z}}$ and assume an arbitrarily large perturbation with single helicity, $\psi (r,\zeta )=\hat {\psi }(r)\cos {\zeta }$, where $\zeta =m\theta -n\varphi$. This perturbation generates an island around the rational surface with $q_s = q(r_s)=m/n$. The X-point of the island is located at $\zeta =0$ and the O-point of the island is located at $\zeta ={\rm \pi}$. However, as we show now, this island is not symmetric in the sense of (4.2), namely $A_{\textrm {sym}}>0$. The island itself consists of a set of nested magnetic surfaces that are described by a flux function $\chi (r,\zeta )$ such that $\boldsymbol {B}\boldsymbol {\cdot }\boldsymbol {\nabla }\chi =0$. One can show that (Fitzpatrick 1995)

(5.1)\begin{equation} \chi(r,\zeta)=\int_r^{r_s}\left( 1-\frac{q}{q_s} \right)B_{\theta}\,{\rm d}r + \hat{\psi}(r)\cos{\zeta} ,\end{equation}

and so the magnetic surfaces in the island, and in particular the island separatrix, are surfaces of constant $\chi$. Furthermore, the asymptotic behaviour of $\hat {\psi }(r)$ around the resonant surface can be obtained from the linear ideal outer solutions (see § 6.9 of Wesson 2004) and is given by

(5.2)\begin{equation} \hat{\psi}^\pm(r) = \psi_s \left[ 1 + \mathcal{A}(r-r_s)\ln{|r-r_s|} + A^{{\pm}}(r-r_s)) \right], \end{equation}

where $\pm$ refers to each side of the resonant surface, $\psi _s$ is the (unknown) amplitude of the perturbation at the resonant surface, $\mathcal {A}$ is given by (3.2a–c) and $A^{\pm }$ is a constant that can be determined numerically by solving for the radial profile of $\hat {\psi }(r)$ outside of the resonant layer. As a matter of fact, the constants $A^{\pm }$ are used to evaluate $\varDelta '$ and $\varSigma '$ as defined in (2.3) and (3.3), since $\varDelta '=A^+ - A^-$ and $\varSigma '=A^+ + A^-$. At the separatrix, we have

(5.3)\begin{equation} \chi(r_+,{\rm \pi})=\chi(r_-,{\rm \pi})\equiv \chi_{{\rm sep}} \approx \chi(r_s,0) = \psi_s, \end{equation}

and hence combining (5.1) and (5.2) into (5.3) we obtain a constraint for $(r_+, r_-)$

(5.4)\begin{equation} g(r_+) = g(r_-), \end{equation}

where

(5.5)\begin{equation} g(r) = \frac{\int_r^{r_s}\left( 1-\dfrac{q}{q_s} \right)B_{\theta}\,{\rm d}r }{2+\mathcal{A}(r-r_s)\ln{|r-r_s|} + A^{{\pm}}(r-r_s)}, \end{equation}

is a known function of the equilibrium. Equation (5.4) can be combined with the constraint for the toroidal flux, $\varPsi _w={\rm \pi} B_{z0}(r_+^2-r_-^2)$, to solve for the positions $r_+$ and $r_-$ of the interfaces that encapsulate the resonant volume. The values obtained for $r_+$ and $r_-$ can then be used to initialize the SPEC equilibrium interfaces. Given this initial geometry, one can show that the expected (maximum) value of the island asymmetry at saturation is

(5.6)\begin{equation} A_{{\rm max}} = 2 \left(\frac{r_s-r_-}{r_+{-}r_s}-1\right). \end{equation}

The value of $A_{\textrm {max}}$ obtained by solving (5.4) for different values of $q_0$ is shown in figure 6 and the agreement with the actual saturated island asymmetry obtained by SpeCyl is quite remarkable. Notice that this prediction is merely based on linear theory.

5.3 Constraining the current profile

Magnetic reconnection is a process through which magnetic energy is converted into kinetic and thermal energy. If the process is fast enough, the magnetic helicity is well conserved (Berger 1999), and this happens for example during sawteeth in tokamaks (Heidbrink & Dang 2000), which occur on quasi-ideal time scales. Resistive tearing modes, however, are instabilities that promote magnetic reconnection at a rate that is fast when compared with the resistive time scale $\tau _{\eta }$, but slow when compared with the Alfvénic time scale $\tau _A$. The reconnection process is then too slow for the magnetic helicity to be a good invariant. Nevertheless, it can be shown that the net toroidal current is a flux function at saturation and is given by (Escande & Ottaviani 2004; Militello & Porcelli 2004)

(5.7)\begin{equation} j_z(\psi) = \langle j_{{\rm eq}} \rangle_{\psi}, \end{equation}

where $\langle \boldsymbol {\cdot } \rangle _\psi$ is the average over a flux surface and $\psi (r,\theta,\varphi )$ is the saturated flux function. Equation (5.7) implies the quasi-invariance of the current profile, the invariance being exact in the limit of small island width. We can thus exploit this constraint on the current profile and run SPEC nonlinear equilibrium calculations by constraining the current in each volume, $I_{\mathrm {vol}}$, to remain constant, and the surface current on each interface, $I_{\mathrm {surf}}$, to remain zero, thereby ensuring that the current profile remains smooth.

5.4 Nonlinear simulation results

For each value of $q_0$, we thus proceed as follows: (i) perform a linear stability analysis with SPEC to find the unstable eigenmode $\boldsymbol {u}_{\lambda }$ and the value of $\varPsi _w^*$ (§ 5.1); (ii) use (5.4) to infer the values of $(r_+,r_-)$ (§ 5.2); (iii) perturb the equilibrium interfaces according to $\boldsymbol {u}_{\lambda }$ and use SPEC to seek a new equilibrium by constraining the current profile (§ 5.3); (iv) measure the width of the island in the resonant volume of the newly found equilibrium. We remark that this procedure has no free parameters. Figure 5 shows two examples of saturated states obtained with SPEC for $q_0=1.3$ (top right) and $q_0=1.6$ (bottom right), which qualitatively look very similar to the ones obtained with SpeCyl (top left and bottom left). Figure 4 shows the saturated island width $w_{\textrm {sat}}$ obtained with SPEC as a function of $q_0$. The qualitative shape and order of magnitude of the expected function $w_{\textrm {sat}}(q_0)$ is reproduced. For small enough islands, the agreement with the predictions from nonlinear theory and resistive MHD simulations is remarkable. The agreement remains good up to island sizes of the order of $10\,\%$ of the minor radius. For increasingly large islands, the equilibrium calculations overestimate the saturation amplitude from nonlinear theory and resistive MHD simulations (the maximum relative error is approximately $25\,\%$). This might be due to the fact that these calculations assume a flat current profile inside the island. Furthermore, the initial positioning of the interfaces around the rational surface, which is determined by solving (5.4), has some approximations, namely that the nonlinear perturbed flux function around the resonant surface has a single helicity and a radial structure given by the linear solution. As a matter of fact, we can measure a posteriori the positions $r_X$, $r_-$ and $r_+$ that characterize the island at saturation in SPEC, and compare them with those obtained with SpeCyl for different values of $q_0$. This is shown in figure 7. In particular, we see that the sign of the asymmetry is well retrieved, namely the property that $r_X - r_- > r_X - r_+$. However, the exact value of $A_{\textrm {sym}}$, (4.2), is quite sensitive to small errors in the values of $r_X$, $r_-$ and $r_+$, and a good agreement is only observed for $q_0\lesssim 1.3$ (see figure 6). We also compare the profiles of each component of $\boldsymbol {B}$ and $\boldsymbol {j}$ at saturation as obtained from SPEC and SpeCyl. Figures 13 and 14 show this comparison for the case $q_0=1.3$, for which the agreement in the island width is the least good. Even for this case, we still observe that SPEC is capable of reproducing (both qualitatively and quantitatively) the structure of the saturated tearing mode. We can also compare the axial current profile at saturation, $j_z(r)$, taken through the O-point and X-point of the island (figure 15). In particular, we indeed see that in SpeCyl the current significantly flattens, but not perfectly, across the island.

Figure 13. Components of $\boldsymbol {B}$ at saturation as obtained from SpeCyl (a,c,e) and SPEC (b,d, f) for $q_0=1.3$. For each component, the color scales are the same between SpeCyl and SPEC.

Figure 14. Components of $\boldsymbol {j}$ at saturation as obtained from SpeCyl (a,c,e) and SPEC (b,d, f) for $q_0=1.3$. For each component, the color scales are the same between SpeCyl and SPEC.

Figure 15. Comparison of the axial current profiles at saturation (red is SPEC, blue is SpeCyl) taken through the O-point (a) and X-point (b) of the island.

Finally, we would like to remark that a convergence analysis has been carried out for the SPEC simulations presented herein. Increasing further the number $N_v=41$ of volumes does not significantly change the values of $w_{\textrm {sat}}$ (data not shown). Similarly, no significant changes are obtained with a further increase in the order of the radial basis functions (Tschebyshev polynomials) used to represent the magnetic vector potential in the volumes. The Fourier resolution, however, seems to play a more important role. Even if the dominant Fourier mode in the saturated island is the $(m,n)=(2,1)$ mode and higher harmonics $(m,n)=(2j,j)$, with $j\in \mathbb {N}$, are present with exponentially decreasing amplitude, the values of $w_{\textrm {sat}}$ can strongly depend on the Fourier resolution used in SPEC, as illustrated in figure 16. At a resolution of $(m,n)=(10,5)$, which is the one used to obtain the results in figure 4 and already implies describing 116 different Fourier modes (Hudson et al. 2020), the relative changes in $w_{\textrm {sat}}$ due to a further increase in Fourier resolution become lower than $5\,\%$. It is also important to note that, for $q_0$ close to 2, the resonant surface is close to the magnetic (and coordinate) axis, and that makes it harder for SPEC to robustly find a solution, mainly because of the singularity in the coordinates at $r=0$. Furthermore, the SPEC code is still numerically fragile in the sense that the Newton-like method employed to find force balance might not be able to converge depending on the initial guess. The main result of this paper is thus the evidence that an equilibrium approach is capable of predicting, to a certain level of accuracy, the saturated state of tearing mode islands in a cylindrical tokamak. The potential of carrying out these predictions very fast is still to be shown by improving the numerical algorithm in SPEC. Nevertheless, when SPEC is able to find a solution, it does so in $\sim 3$ CPU hours, while the time-dependent resistive MHD calculations performed with SpeCyl can take hours to days, mainly depending on the value of the Lundquist number $S$, which strongly affects the simulation time required to reach saturation, and on the number of Fourier harmonics included in the calculation. The SpeCyl predictions are, however, more accurate, as was shown in figure 4.

Figure 16. Saturated island width (upper curve, left axis) obtained from SPEC nonlinear calculations with $q_0=1.05$ as a function of the poloidal Fourier resolution $m$. The toroidal Fourier resolution is always $n=m/2$. The relative error defined here as $(w_{\textrm {sat}}(m+1)-w_{\textrm {sat}}(m))/w_{\textrm {sat}}(m)$ is also shown (lower curve, right axis).

6 Conclusions

In this paper we have investigated the nonlinear saturation of resistive tearing modes in a zero-pressure cylindrical tokamak. We have validated the results of nonlinear visco-resistive MHD simulations against an exact nonlinear saturation theory, showing excellent agreement for a wide range of equilibria. Furthermore, we have shown that an equilibrium approach can be used to predict, without free parameters, the saturated state of these modes to a certain degree of accuracy, with a very good agreement for sufficiently small islands. The natural next step is to perform similar tearing mode saturation studies in toroidal geometries. In stellarators, classical tearing modes have been observed (Bartlett et al. 1980) and modelled (Nikulsin et al. 2022), and we plan on investigating them with the SPEC code. In tokamaks, neoclassical tearing modes play a crucial role in the potential trigger of plasma disruptions, and we plan on extending this approach with the addition of the bootstrap current effects, which are already present in the SPEC code (Baillod et al. 2021, 2023). There are also toroidal corrections to classical tearing modes in tokamaks (Glasser, Greene & Johnson 1975), which have already been reproduced linearly with the SPEC code (Kumar et al. 2023). We plan on verifying whether the effect of these corrections on the nonlinear saturation of classical tearing modes can be retrieved with an equilibrium approach as well. Finally, the effect of rotation on tearing modes could also be studied by exploiting recent extensions of the stepped-pressure equilibrium model that include the presence of flows (Qu et al. 2020).