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 Reference Freidberg2014), but as we will see in the next section, they can be resistive unstable (Furth, Rutherford & Selberg Reference Furth, Rutherford and Selberg1973).

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. Reference Hudson, Dewar, Dennis, Hole, McGann, Von Nessi and Lazerson2012, Reference Hudson, Loizu, Zhu, Qu, Nührenberg, Lazerson, Smiet and Hole2020), 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. Reference Hudson, Dewar, Dennis, Hole, McGann, Von Nessi and Lazerson2012). 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. Reference Dennis, Hudson, Dewar and Hole2013a). 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. Reference Baillod, Loizu, Qu, Kumar and Graves2021). 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. Reference Furth, Rutherford and Selberg1973)

(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).

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. Reference Furth, Rutherford and Selberg1973) and nonlinear (Rutherford Reference Rutherford1973) 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. Reference Poyé, Agullo, Smolyakov, Benkadda and Garbet2014). For a cylindrical tokamak, Arcis et al. (Reference Arcis, Escande and Ottaviani2006) 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. Reference Furth, Rutherford and Selberg1973), 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 Reference Cappello2004; Bonfiglio et al. Reference Bonfiglio, Chacón and Cappello2010) 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 (Reference Cappello and Biskamp1996) and Cappello (Reference Cappello2004).

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. Reference Poyé, Agullo, Smolyakov, Benkadda and Garbet2014) 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. Reference Teng, Ferraro, Gates, Jardin and White2016).

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. Reference Loizu, Huang, Hudson, Baillod, Kumar and Qu2020). 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. Reference White, Monticello, Rosenbluth and Waddell1977), 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 Reference Fitzpatrick1995)

(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 Reference Wesson2004) 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 Reference Berger1999), and this happens for example during sawteeth in tokamaks (Heidbrink & Dang Reference Heidbrink and Dang2000), 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 Reference Escande and Ottaviani2004; Militello & Porcelli Reference Militello and Porcelli2004)

(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. Reference Hudson, Loizu, Zhu, Qu, Nührenberg, Lazerson, Smiet and Hole2020), 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).