Constant-$\psi$ Magnetic Island Evolution

Let us consider the nonlinear evolution of a chain of constant-$\psi$ magnetic islands whose thickness is much greater than the linear layer thickness, but much less than the thickness of the equilibrium current sheet. In other words,

$$\hat{\delta} \ll \hat{W}\ll 1,$$ (9.88)

where $\hat{\delta} = S^{-2/5}$.

Writing $\psi = -B_0\,a\,\hat{\psi}$, $J=-(B_0/a)\,(1+\hat{J})$, $\phi=(a/k\,\tau_R)\,\skew{3}\hat{\phi}$, $U=(1/a\,k\,\tau_R)\,\hat{U}$, and $t=\tau_R\,\hat{t}$, the reduced-MHD equations, (9.11)–(9.14), become

$$\frac{\partial\hat{\psi}}{\partial\hat{t}} =\{\skew{3}\hat{\phi},\hat{\psi}\}+\hat{J},$$ (9.89)

$$\frac{\partial\hat{U}}{\partial\hat{t}} = \{\skew{3}\hat{\phi},\hat{U}\}+S^{2}\,\{\hat{J},\hat{\psi}\},$$ (9.90)

$$\frac{\partial^2\hat{\psi}}{\partial \hat{x}^2} = 1 + \hat{J},$$ (9.91)

$$\hat{U} =\frac{\partial^2 \skew{3}\hat{\phi}}{\partial \hat{x}^2}$$ (9.92)

in the limit $\vert\hat{x}\vert\ll 1$, where

$$\{A, B\} \equiv \frac{\partial A}{\partial\hat{x}}\,\frac{\partia... ...rtial\xi} - \frac{\partial A}{\partial\xi}\,\frac{\partial B}{\partial\hat{x}}.$$ (9.93)

In the island region, $\hat{x}\sim \hat{W}$, $\xi\sim 1$, $\hat{\psi}\sim \hat{W}^{\,2}$, $\hat{J}\sim \hat{W}$ (this is consistent with ${\mit\Delta}'\sim 1$), $\hat{t}\sim \hat{W}$ (as will become apparent), $\skew{3}\hat{\phi}\sim 1$ (as will also become apparent), and $\hat{U}\sim 1/\hat{W}^{\,2}$. It follows that all terms in Equations (9.89) and (9.92) are of the same order of magnitude. On the other hand, the terms involving $\hat{U}$ in Equation (9.90) are smaller than the other term by a factor $(\hat{\delta}/\hat{W})^5$. Finally, the term involving $\hat{J}$ in Equation (9.91) is smaller than the other terms by a factor $\hat{W}$. Thus, to lowest order, Equations (9.89)–(9.92) reduce to

$$\frac{\partial\hat{\psi}}{\partial\hat{t}} =\{\skew{3}\hat{\phi},\hat{\psi}\}+\hat{J},$$ (9.94)

$$\{\hat{J},\hat{\psi}\} =0,$$ (9.95)

$$\frac{\partial^2\hat{\psi}}{\partial \hat{x}^2} = 1.$$ (9.96)

It is clear that the terms involving plasma inertia (i.e., the terms involving $\hat{U}$) in the vorticity evolution equation (9.90) become negligible as soon as the island width exceeds the linear layer width, leaving (9.95), which corresponds to the force balance criterion $\nabla\times({\bf j}\times {\bf B} +\nabla p)=0$. Thus, a nonlinear magnetic island chain is essentially a $y$-dependent magnetic equilibrium.

Equation (9.96) can be integrated to give

$$\hat{\psi}(\hat x,\xi,\hat{t})=\frac{\hat{x}^2}{2}+\skew{3}\hat{\mit\Psi}(\hat{t})\,\cos\xi,$$ (9.97)

which is identical to the constant-$\psi$ result (9.83). It follows that the constant-$\psi$ approximation is valid provided that $\vert\hat{J}\vert\ll 1$. In other words, provided that the perturbed current density in the island region is small compared to the equilibrium current density.

Equation (9.95) implies that

$$\hat{J} = \hat{J}(\hat{\psi}).$$ (9.98)

In other words, the current density in the island region is a flux-surface function (i.e., it is constant on magnetic field-lines). Equations (9.94), (9.97), and (9.98) can be combined to give

$$\frac{d\skew{3}\hat{\mit\Psi}}{d\hat{t}}\,\cos\xi = \{\skew{3}\hat{\phi},\hat{\psi}\} + \hat{J}(\hat{\psi}).$$ (9.99)

It is helpful to define $X=4\,\hat{x}/\hat{W}$. Writing $\hat{\psi}(\hat{x},\xi,\hat{t}) = \skew{3}\hat{\mit\Psi}(\hat{t})\,{\mit\Omega}(X,\xi)$, we find that

$${\mit\Omega}(X,\xi) = \frac{X^2}{2} + \cos\xi.$$ (9.100)

Thus, the contours of ${\mit\Omega}(X,\xi)$ map out the magnetic field-lines in the vicinity of the island chain. The region inside the magnetic separatrix corresponds to $-1\leq{\mit\Omega}\leq 1$, whereas the region outside the separatrix corresponds to ${\mit\Omega}>1$. Let us transform into the new coordinate system $s\equiv{\rm sgn}(X)$, ${\mit\Omega}$, and $\xi$. It follows that $X=s\sqrt{2\,({\mit\Omega}-\cos\xi)}$. When written in terms of the new coordinates, the plasma Ohm's law, (9.99), becomes

$$\frac{d\skew{3}\hat{\mit\Psi}}{d\hat{t}}\,\cos\xi =-\frac{\hat{W}... ...ew{3}\hat{\phi}}{\partial\xi}\right\vert _{\mit\Omega} + \hat{J}({\mit\Omega}).$$ (9.101)

It is helpful to define the flux-surface average operator,

$$\langle A(s,{\mit\Omega},\xi)\rangle =\left\{ \begin{array}{l... ...os\xi)}}\,\frac{d\xi}{2\pi}&&{\mit\Omega}>1 \end{array}\right.,$$ (9.102)

where $\xi_0=\cos^{-1}({\mit\Omega})$. The flux-surface average operator is designed to annihilate the first term on the right-hand side of Equation (9.101). Thus, the flux-surface average of this equation yields

$$\hat{J}({\mit\Omega}) = \frac{d\skew{3}\hat{\mit\Psi}}{d\hat{t}}\,\frac{\langle\cos\xi\rangle}{\langle 1\rangle}.$$ (9.103)

Let

$$\widetilde{\cos\xi} \equiv \cos\xi - \frac{\langle \cos\xi\rangle}{\langle 1\rangle},$$ (9.104)

which implies that $\langle\widetilde{\cos\xi}\rangle = 0$. Equations (9.101) and (9.103) can be combined to give

$$\skew{3}\hat{\phi}(s,{\mit\Omega},\xi,\hat{t}) = -\frac{4}{\hat{W... ...nt_0^\xi \frac{\widetilde{\cos\xi'}}{\sqrt{2\,({\mit\Omega}-\cos\xi')}}\,d\xi'.$$ (9.105)

According to Equations (9.37), (9.91), and (9.102), asymptotic matching between the island solution and the solution in the outer region (i.e., the region $\vert\hat{x}\vert\gg \hat{W}$) yields

$${\mit\Delta}' = \frac{2}{\skew{3}\hat{\mit\Psi}}\int_{-\infty}^\infty \oint \fr... ...infty}^\infty \oint \hat{J}({\mit\Omega})\,\cos\xi\,d\hat{x}\,\frac{d\xi}{2\pi}$$

$$= \frac{\hat{W}}{\skew{3}\hat{\mit\Psi}}\int_{-1}^\infty \hat{J}({\mit\Omega})\,\langle\cos\xi\rangle\,d{\mit\Omega}.$$ (9.106)

Note that it is necessary to specifically project the $\cos\xi$ component out of $\hat{J}({\mit\Omega})$ because $\hat{J}({\mit\Omega})$ is a nonlinear function that possesses many $\cos(m\,\xi)$ components, where $m$ ranges from $1$ to $\infty$. (See later.) Equations (9.103) and (9.106) can be combined to give

$$I_1\,\frac{d\hat{W}}{d\hat{t}} = {\mit\Delta}',$$ (9.107)

where

$$I_1 =2\int_0^\infty \frac{\langle \cos\xi\rangle^2}{\langle 1\rangle}\,d{\mit\Omega} = 0.8227.$$ (9.108)

(See Table 9.1.) Equation (9.107) confirms that $\hat{t}\sim \hat{W}$ (provided that ${\mit\Delta}'\sim 1$), as was previously assumed. Moreover, Equations (9.84), (9.103) and (9.107) demonstrate that the constant-$\psi$ approximation (which requires $\vert\hat{J}\vert\ll 1$) is valid provided that

$${\mit\Delta}'\,\hat{W}\ll 1.$$ (9.109)

Note that this criterion becomes harder to satisfy as the island width grows, which implies that if the tearing mode is in the constant-$\psi$ regime when it enters the nonlinear regime (i.e., if ${\mit\Delta}'\,\hat{\delta}\ll 1$) then it may spontaneously leave the constant-$\psi$ regime as it subsequently evolves in time. (However, if ${\mit\Delta}'\,\hat{W}$ approaches unity when ${\mit\Delta}'\sim 1$ then this indicates a breakdown of asymptotic matching, due to the fact that the island width is no longer small compared to the current sheet thickness, rather than a breakdown of the constant-$\psi$ approximation.) When written in unnormalized form, the island width evolution equation, (9.107), yields the Rutherford island width evolution equation  (Rutherford 1973),

$$I_1\,\tau_R\,\frac{d(W/a)}{dt} = {\mit\Delta}',$$ (9.110)

where $W=a\,\hat{W}$ is the island width in $x$. According to the Rutherford equation, the growth of a constant-$\psi$ tearing mode slows down as it enters the nonlinear regime (i.e., as the island width exceeds the linear layer width). Indeed, the tearing mode transitions from growing exponentially in time on the hybrid timescale $\tau_H^{2/5}\,\tau_R^{3/5}$ to growing algebraically in time on the much longer timescale $\tau_R$.

Figure: 9.6 Contours of the normalized perturbed current density distribution, $\langle\cos\xi\rangle/\langle 1\rangle$, in the vicinity of a constant-$\psi$ magnetic island chain. Positive/negative values are indicated by solid/dashed contours.

The current density in the island region that is specified in Equation (9.103) can be written in the form

$$\hat{J}(\hat{x},\xi) = \sum_{m=1,\infty}J_m(\hat{x})\,\cos(m\,\xi),$$ (9.111)

where the $J_m(\hat{x})$ are even functions of $\hat{x}$ that are similar in magnitude to one another. (Note that there is no $m=0$ harmonic.) This is clear from Figure 9.6, which shows contours of the normalized perturbed current density distribution, $\langle\cos\xi\rangle/\langle 1\rangle$. [See Equation (9.103).] It can be seen that the current density is mostly confined to the interior of the magnetic separatrix, and becomes particularly large on the separatrix itself. (In fact, the current density blows up logarithmically on the separatrix.) Clearly, such a current distribution cannot be represented as $J_1(\hat{x})\,\cos\xi$. In other words, the current density distribution is multi-harmonic (i.e., it is not dominated by the $m=1$ harmonic, but contains substantial contributions from the $m>1$ harmonics). It would seem reasonable, therefore, to write the solution to (9.96) in the form

$$\hat{\psi}(\hat{x},\xi,\hat{t})=\frac{\hat{x}^2}{2}+\sum_{m=1,\infty}\skew{3}\hat{\mit\Psi}_m(\hat{t})\,\cos(m\,\xi),$$ (9.112)

where the $\skew{3}\hat{\mit\Psi}_m(\hat{t})$ are independent of $\hat{x}$. However, when we actually wrote the solution to this equation, in Equation (9.97), we omitted the higher harmonics (i.e., the $\skew{3}\hat{\mit\Psi}_m$ for $m>1$). Let us now investigate under which circumstances this approximation can be justified. Let us, first of all, assume that the allowed wavenumbers are quantized: that is, $k=k_m$, for $m=1,\infty$, where $k_m=m\,k_1$. (Here, $k_1$ is what we previously referred to as $k$.) Obviously, the easiest way to justify the quantization of the allowed wavenumbers is to assume that the equilibrium current sheet has a finite length $2\pi/k_1$ in the $y$-direction. Asymptotic matching between the island solution and the solution in the outer region yields

$$I_m\,\frac{d\hat{W}}{d\hat{t}} = {\mit\Delta}'_m\,\frac{\skew{3}\hat{\mit\Psi}_m}{\skew{3}\hat{\mit\Psi}_1},$$ (9.113)

where

$$I_m = 2\int_{-1}^\infty \frac{\langle\cos(m\,\xi)\rangle\,\langle \cos\xi\rangle}{\langle 1\rangle}\,d{\mit\Omega}.$$ (9.114)

Note that Equation (9.113) is the multi-harmonic generalization of Equation (9.107), Here, $\skew{3}\hat{\mit\Psi}_1$ is what we previously referred to as $\skew{3}\hat{\mit\Psi}$, and $\hat{W}= \skew{3}\hat{\mit\Psi}_1^{1/2}$. Moreover, ${\mit\Delta}_m'$ is the tearing stability index calculated with the wavenumber $k_m$. (So ${\mit\Delta}_1'$ is what we previously referred to as ${\mit \Delta }'$.) Equation (9.113) yields

$$\frac{\skew{3}\hat{\mit\Psi}_m}{\skew{3}\hat{\mit\Psi}_1}=\frac{{\mit\Delta}_1'}{{\mit\Delta}_m'}\,\frac{I_m}{I_1}.$$ (9.115)

Figure: 9.7 Integrands for $I_n$ integrals.

It is helpful to define $k =\sqrt{(1+{\mit\Omega})/2}$. Thus, $k=0$ at the centers of the magnetic islands, $k=1$ on the magnetic separatrix, and $k>1$ in the region outside the separatrix. It can be demonstrated that

$$\langle 1\rangle \begin{align*}= \left\{ \begin{array}{lll} K(\pi/2,k)/\pi&~~~~&0\leq k < 1\\ [2ex] K(\pi/2,1/k)/k\,\pi&&k>1 \end{array}\right.,\end{align*}$$ (9.116)

$$\langle\cos\xi\rangle = \begin{align*}\left\{ \begin{array}{lll} [K(\pi/2,k)-2\,E(\pi/2,k)]/\pi &~~~~&0\... ...\,K(\pi/2,1/k)-2\,k^2\,E(\pi/2,1/k)]/k\,\pi&&k>1 \end{array}\right.,\end{align*}$$ (9.117)

$$\langle \cos(m\,\xi)\rangle= \begin{align*}\left\{ \begin{array}{lll} \int_0^{\pi/2}\frac{\cos[2\,m\,\cos^{-1... ...t{k^2-\sin^2\varphi}}\,\frac{d\varphi}{\pi}&&k>1 \end{array}\right.,\end{align*}$$ (9.118)

where

$$K(\varphi,k) = \int_0^\varphi \frac{du}{\sqrt{1-k^2\sin^2 u}},$$ (9.119)

$$E(\varphi,k) = \int_0^\varphi \sqrt{1-k^2\sin^2 u}\,\,du$$ (9.120)

are elliptic integrals (Abramowitz and Stegun 1965). Hence, we can write

$$I_m = \int_0^\infty K_m(k)\,dk,$$ (9.121)

where

$$K_m(k) = 8\,k\, \frac{\langle\cos(m\,\xi)\rangle\,\langle \cos\xi\rangle}{\langle 1\rangle}.$$ (9.122)

Figure 9.7 shows the integands $K_m(k)$. Note that all integrands are singular at the magnetic separatrix ($k=1$). However, the singularities are logarithmic in nature, and, therefore, integrable. Note, further, that the higher order (i.e., $m>1$) integrands are of similar magnitude to the $m=1$ integrand, which confirms that the $J_m(\hat{x})$, for $m>1$, appearing in Equation (9.111), are of similar magnitude to $J_1(\hat{x})$. In other words, the current density distribution is truly multi-harmonic. However, it can be seen, from Figure 9.7, that the $m=1$ integrand is always positive, whereas the $m>1$ integrands oscillate about zero. It is not surprising, therefore, that the $I_m$ for $m>1$ are much smaller in magnitude than $I_1$. In fact, as is shown in Table 9.1, the $I_m$ for $m>1$ are, at least, 10 times smaller than $I_1$. It follows from Equation (9.115) that the single-harmonic approximation for $\hat{\psi}$ (i.e., the neglect of the $\skew{3}\hat{\mit\Psi}_m$ for $m>1$) is justified as long as $\vert{\mit\Delta}_1'/{\mit\Delta}_m'\vert$ does not exceed the critical value $\vert I_1/I_m\vert$, for $m>1$. As shown in Table 9.1, this is a comparatively easy criterion to meet, especially if the $m=1$ harmonic is close to marginal stability (i.e., ${\mit\Delta}_1'>0$ is small compared to unity).

Table: 9.1 The $I_m$ integrals for $m=1$ to $5$.

$m$	1	2	3	4	5
$I_m$	$8.227\times 10^{-1}$	$-8.529\times 10^{-2}$	$-1.659\times 10^{-2}$	$-6.410\times 10^{-3}$	$-3.228\times 10^{-3}$
$\vert I_1/I_m\vert$	$1.000\times 10^{0}$	$9.645\times 10^{0}$	$4.959\times 10^1$	$1.283\times 10^2$	$2.548\times 10^2$

Figure: 9.8 Equally-spaced contours of the normalized stream-function, $H(\hat{x},\xi)$, in the vicinity of a constant-$\psi$ magnetic island chain. Positive/negative values are indicated by solid/dashed contours.

Equations (9.84), (9.105), and (9.107) yield

$$\skew{3}\hat{\phi}(s,k,\xi) = \frac{{\mit\Delta}'}{I_1}\,H(s,k,\xi),$$ (9.123)

where

$$H(s,k,\xi)=s\left\{ \begin{array}{lll} \frac{F(\varphi,k)\,E(... ...,1/k)\,E(\varphi',1/k)]}{F(\pi/2,1/k)}&&k>1 \end{array}\right.,$$ (9.124)

and

$$\varphi = \frac{\pi}{2}-{\rm sgn}(\xi)\,\cos^{-1}[\cos(\xi/2)/k],$$ (9.125)

$$\varphi' = \frac{\pi}{2}-\frac{\xi}{2}.$$ (9.126)

Note that Equation (9.123) confirms that $\skew{3}\hat{\phi}\sim 1$ in the island region (assuming that ${\mit\Delta}'\sim 1$), as was previously assumed. Figure 9.8 shows contours of the normalized stream-function, $H$, in the $\hat{x}$-$\xi$ plane. It can be seen that the flow pattern associated with magnetic reconnection, which is clearly strongly multi-harmonic, is concentrated at the magnetic separatrix, being particularly large at the magnetic X-points (Biskamp 1993). In fact, the flow velocity has a logarithmic singular on the magnetic sepatratrix, indicating that plasma inertia cannot be neglected in the immediate vicinity of the sepatatrix. Consequently, an inertial layer [i.e., a layer in which plasma inertia cannot be neglected in Equation (9.90)], whose width is of order the linear layer width, develops on the separatrix. Fortunately, the development of such a layer does not invalidate the results obtained in this section (Edery, et al., 1983). It is interesting to note that the Rutherford island width evolution equation, (9.110), can be derived without explicitly calculating the flow pattern. Nevertheless, the flow pattern is implicitly specified in the analysis.

The island width evolution equation (9.107) seems to indicate that if ${\mit\Delta}'> 0$ then the island width, $\hat{W}$, grows without limit. In fact, this is not the case. If we perform the asymptotic matching between the solution in the island region and the solution in the outer region more carefully, taking the finite width of the magnetic island chain into account, then the tearing stability index, ${\mit \Delta }'$, which is defined in Equation (9.34), is replaced by (White, et al. 1977; Biskamp 1993)

$${\mit\Delta}'(\hat{W}) =\frac{1}{\hat{\psi}(0)}\left[\frac{d\hat{\psi}(\bar{W}/2)}{d\hat{x}}- \frac{d\hat{\psi}(-\bar{W}/2)}{d\hat{x}}\right].$$ (9.127)

Here, $\bar{W}=2\,\hat{W}/\pi$ is the average (over $\xi$) width of the magnetic separatrix, and $\hat{\psi}(\hat{x})$ is a solution of the tearing mode equation, (9.31). Making use of Equation (9.36), we deduce that

$${\mit\Delta}'(\hat{W}) = {\mit\Delta}'(0) - \alpha\,\hat{W} + {\cal O}(\hat{W}^2),$$ (9.128)

where

$$\alpha = \frac{4}{\pi}\,(2-\hat{k}^2),$$ (9.129)

for the specific plasma equilibrium discussed in Section 9.3. Here, ${\mit\Delta}'(0)$ is specified in Equation (9.39). Thus, the modified island width evolution equation takes the form

$$I_1\,\frac{d\hat{W}}{d\hat{t}} \simeq {\mit\Delta}'(0)-\alpha\,\hat{W}.$$ (9.130)

It is clear from this equation that nonlinear growth of the tearing mode slows down, as the island width increases, and eventually stops when the island width attains the value

$$\hat{W}_s \simeq \frac{{\mit\Delta}'(0)}{\alpha}.$$ (9.131)

For the specific equilibrium discussed in Section 9.3, the so-called saturated island width is

$$\hat{W}_s = \frac{\pi}{2}\,\frac{1-\hat{k}^2}{\hat{k}\,(2-\hat{k}^2)}.$$ (9.132)

The previous expression is only accurate when $\hat{W}_s\ll 1$ (i.e., when the saturated island width is much smaller than the width of the current sheet). This is only the case when ${\mit\Delta}'\ll 1$ (i.e., when the tearing mode is close to marginal stability). However, it seems reasonable to deduce that if ${\mit\Delta}'\sim 1$ then the tearing mode eventually attains a steady-state with a saturated island width that is comparable to the width of the equilibrium current sheet (i.e., $\hat{W}\sim 1$). In other words, the tearing mode completely changes the topology of the current sheet's magnetic field on a timescale that is of order $\tau_R$. Equation (9.127) is only approximate. However, more rigorous calculations give essentially the same result (Thyagararja, 1981; Escande & Ottaviani 2004; Militello & Porcelli 2004; Hastie, et al. 2005).

At first sight, the time evolution of a constant-$\psi$ tearing mode in the nonlinear regime seems completely different to that in the linear regime. However, it turns out to be comparatively easy to formulate a theory that takes both regimes into account. The time evolution of the tearing mode in the linear regime is specified by

$$\frac{d\skew{3}\hat{\mit\Psi}}{d\hat{t}} = \gamma\,\tau_R\, \skew{3}\hat{\mit\Psi}.$$ (9.133)

Making use of Equation (9.60), this equation reduces to

$$\hat{\delta}\,\frac{d\skew{3}\hat{\mit\Psi}}{d\hat{t}} = {\mit\Delta}'\, \skew{3}\hat{\mit\Psi},$$ (9.134)

where

$$\hat{\delta}=\left[\frac{2\pi\,\Gamma(3/4)}{\Gamma(1/4)}\right]^{4/5}{\mit\Delta}'^{1/5}\,S^{-2/5}$$ (9.135)

is the exact (normalized) linear layer width. The normalized Rutherford island width evolution equation, (9.107), can be written

$$\frac{I_1\,\hat{W}}{2}\,\frac{d\skew{3}\hat{\mit\Psi}}{d\hat{t}} = {\mit\Delta}'\, \skew{3}\hat{\mit\Psi},$$ (9.136)

where use has been made of Equation (9.84). It can be seen, by comparison with Equation (9.134), that a nonlinear magnetic island evolves in time in an analogous manner to a constant-$\psi$ linear layer whose (normalized) thickness is $I_1\,\hat{W}/2$. Thus, the essential nonlinearity in the nonlinear regime arises because the effective layer width scales as the square root of the mode amplitude. [See Equation (9.84).] Recall that Equation (9.134) is valid when $\hat{\delta}\gg\hat{W}$, whereas Equation (9.136) is valid when $\hat{\delta}\ll \hat{W}$. Thus, Equations (9.134) and (9.136) can be combined to give the composite evolution equation

$$\left(\hat{\delta} + \frac{I_1\,\hat{W}}{2}\right)\frac{d\skew{3}\hat{\mit\Psi}}{d\hat{t}} = {\mit\Delta}'\, \skew{3}\hat{\mit\Psi}$$ (9.137)

that interpolates between the linear and nonlinear regimes. As is clear, from the previous equation, the magnetic reconnection rate decelerates in a smooth fashion as the island width exceeds the linear layer width.

Note, finally, that the type of slow magnetic reconnection, mediated by tearing modes that grow and eventually saturate, described in this section is observed on a routine basis in tokamak plasmas (Wesson 2011).