# Large scale instabilities in two-dimensional magnetohydrodynamics

## I Introduction

Magnetic reconnection in two-dimensional magnetohydrodynamics (MHD) is one of the most intriguing problems in plasma physics, originally motivated by the observation of sudden and rapid release of energy occurring for instance in solar flares and tokamak disruptions. These phenomena are characterized by a change in the topology of the magnetic field lines and the formation of island-like structures. Resistive instabilities are often addressed as responsible of magnetic reconnection . Their development is a consequence of a resistive boundary layer formed nearby a magnetic field neutral line, where the ideal constraint of frozen magnetic flux is relaxed. Their growth rates typically depend on fractional powers of the dissipative coefficients (through the dimensionless inverse Lundquist number and Reynolds number). When nonlinear effects are taken into account, the growth of resistive instabilities slows down, eventually reaching saturation . Numerical studies  have shown that the nonlinear phase of evolution is characterized by the generation of large scale magnetic structures, by means of a coalescence process of small magnetic flux patches.

In this Paper we present a new kind of instability which sets in when the influence of dissipative terms cannot be restricted to a limited domain and the boundary layer approximation does not hold anymore. This instability can be viewed as the MHD counterpart of the large scale hydrodynamical instability which is known to develop in highly anisotropic flows . Large scale instabilities are formally associated to the existence of negative values of the eddy-viscosity and/or eddy-resistivity. The exact values of the turbulent dissipative coefficients can be explicitly derived from the knowledge of the equations of motion and of the basic equilibrium state, by means of multiple scale analysis. In the case of a sheared magnetic field embedded in a motionless fluid, a negative eddy-resistivity instability can develop for high enough Lundquist/Reynolds numbers, causing the formation of a chain of magnetic islands. As observed for resistive instabilities, also in this case, during the nonlinear evolution, the number of magnetic islands decreases, due to the merging of pairs of close islands. This coalescence process eventually leads to the formation of a single island, whose width is determined by the largest scale allowed in the system.

In Section II, we briefly introduce the fluid equations of magnetohydrodynamics and their basic equilibrium states. In Section III, the behavior of large scale perturbations is investigated making use of multiple-scale analysis and exact expressions for eddy-viscosity and eddy-resistivity are derived. The main result is the appearance of large scale transverse instabilities, associated to negative values of the renormalized dissipative coefficients. In Section IV, we focus on the case of marginal instability for which we obtain the effective equation for the large scale behavior and we show that the full nonlinear regime is characterized by the evolution towards a fixed point. In Section V, we present the results of direct numerical simulations (DNS), which display a clear quantitative agreement with the predictions of multiple-scale analysis, both in the linear and in the nonlinear phase.

## II MHD equations and basic equilibria

MHD equations are relevant in many different physical contests, such as astrophysics, laboratory plasma physics or magnetized liquid metal dynamics. The applicability of this model, which is a fluid description of plasma, relies on the assumption that all the lengthscales under consideration must largely exceed the ion Larmor radius. In a strong external magnetic field which is oriented along $`z`$, $`B_zB_{}`$, the motion becomes almost two-dimensional and the MHD model is well approximated by the 2D MHD equations for the magnetic flux function $`\psi `$ associated to the planar magnetic field ($`𝐁_{}=𝐞_z\times \psi `$) and for the stream function $`\phi `$ of the incompressible planar flow ($`𝐯_{}=𝐞_z\times \phi `$) :

$`{\displaystyle \frac{\psi }{t}}+[\phi ,\psi ]`$ $`=`$ $`\eta (^2\psi J_0),`$ (1)   
$`{\displaystyle \frac{^2\phi }{t}}+[\phi ,^2\phi ]`$ $`=`$ $`[\psi ,^2\psi ]+\nu ^4\phi ,`$ (2)   
where, following a standard notation, the convective terms are written as Jacobian operators($`[f,g]=_xf_yg_xg_yf`$). The forcing $`J_0=^2\psi _0`$ represents an input of magnetic energy preventing the decay due to resistive ($`\eta `$) and viscous ($`\nu `$) dissipation.

The above equations have been normalized with respect to the characteristic macroscopic length $`L`$, magnetic field $`\overline{B}`$ and Alfvèn time $`\tau _A=L/v_A`$, while $`\eta `$ and $`\nu `$ are, respectively, the inverse Lundquist number ($`\eta c^2/4\pi v_AL`$) and the inverse Reynolds number ($`\nu /v_AL`$). The Alfvèn velocity, $`v_A=\overline{B}/(4\pi mn)^{1/2}`$, is the velocity of small amplitude waves, propagating along the magnetic field $`\overline{B}`$ in a uniform plasma with density $`nm`$.

The equations (1,2) admit a basic equilibrium state $`\psi =\psi _0=(x),\phi =0`$, where $``$ is any function. We consider $`(x)=\mathrm{cos}x`$, in order to investigate the stability of a sheared magnetic field in a motionless conducting fluid, in a slab geometry with periodic boundary conditions. Such a configuration is widely used to study the evolution of reconnecting modes. It is known  that the slab geometry can be trusted as realistic and resistive instabilities can develop only if the slab aspect ratio ($`L_x/L_y`$) is less than one. For this reason it is of great interest the question related to the stability of the above configuration with respect to transverse perturbations on scales larger than the typical scale of the basic magnetic flux. These large scale instabilities are eventually responsible for the onset of the inverse cascade of the square magnetic potential $`A=\psi ^2`$ in fully developed 2D MHD turbulence .

## III Multiple-scale analysis

As a first step, using a standard multiple-scale technique, we derive analytically the large-scale magnetic flux $`\mathrm{\Psi }`$ and stream function $`\mathrm{\Phi }`$ equations. The multiscale method is based on the idea of exploiting the separation of scales as a perturbative parameter. We consider a periodic box width $`L_{box}`$ much larger than the basic magnetic flux typical length-scale ($`L_{box}L_x`$) in order to follow the dynamics on spatial scales of order $`O(1/\epsilon )L_x`$. Beside the fast variables ($`x,y,t`$) on which the basic flow evolves, a set of slow variables ($`X=\epsilon x,Y=\epsilon y,T=\epsilon ^2t`$) can be introduced. According to this choice the differential operators appearing in (1,2) are transformed to

$$_i_i+\epsilon _i,_t_t+\epsilon ^2_T.$$
(3)   
Expanding perturbatively in $`\epsilon `$ the fields we obtain:

$`\psi `$ $`=`$ $`\psi ^{(0)}(x,X,y,Y,t,T)+\epsilon \psi ^{(1)}(x,X,y,Y,t,T)+\epsilon ^2\psi ^{(2)}\mathrm{}`$ (4)   
$`\phi `$ $`=`$ $`\phi ^{(0)}(x,X,y,Y,t,T)+\epsilon \phi ^{(1)}(x,X,y,Y,t,T)+\epsilon ^2\phi ^{(2)}\mathrm{}`$ (5)   
The scaling of the slow time $`T`$ is suggested by physical hints: we are looking for a diffusive behavior of large scales which takes place on times $`O(\epsilon ^2)`$. It is worth noticing that in general the large-scale MHD dynamics is first order in time and space (the well-known $`\alpha `$ effect) , but it can be shown that this is not the case for parity-invariant basic configurations .

By substituting (3) and (5) in (1) and (2) and by equating the same power of $`\epsilon `$, one easily finds a hierarchy of equations in which perturbations belonging to different order of expansion appear coupled and depend on fast and slow variables. The dependence on the fast time variables can be discarded by observing that it reduces to a transient not affecting the long-time behavior, thanks to the fact that the forcing and the basic flux are time-independent (a rigorous proof needs the construction of a Poincarè inequality). We look for solutions with the same periodicities of the basic magnetic flux in $`L_{box}`$. At each step we can distinguish the pure large field contribution from the small scale oscillating part:

$$\psi ^{(k)}=\mathrm{\Psi }^{(k)}(X,Y,T)+\stackrel{~}{\psi }^{(k)}(X,Y,T,x),\phi ^{(k)}=\mathrm{\Phi }^{(k)}(X,Y,T)+\stackrel{~}{\phi }^{(k)}(X,Y,T,x).$$

Equations have to be solved recursively because solutions of lower order appear as coefficients in the following steps of the hierarchy. At each order one has to test the validity of the solvability condition. Indeed, the equation for the large scale magnetic flux $`\mathrm{\Psi }^{(0)}`$ is obtained as solvability condition at order $`\epsilon ^2`$, while the equation for the large scale vorticity $`^2\mathrm{\Phi }^{(0)}`$ comes out at order $`\epsilon ^4`$:

$`{\displaystyle \frac{\mathrm{\Psi }^{(0)}}{T}}+[\mathrm{\Phi }^{(0)},\mathrm{\Psi }^{(0)}]`$ $`=`$ $`\eta ^2\mathrm{\Psi }^{(0)}{\displaystyle \frac{1}{2\nu }}{\displaystyle \frac{^2\mathrm{\Psi }^{(0)}}{Y^2}},`$ (6)   
$`{\displaystyle \frac{^2\mathrm{\Phi }^{(0)}}{T}}+[\mathrm{\Phi }^{(0)},^2\mathrm{\Phi }^{(0)}]`$ $`=`$ $`[\mathrm{\Psi }^{(0)},^2\mathrm{\Psi }^{(0)}]+{\displaystyle \frac{1}{2}}{\displaystyle \frac{^2}{XY}}\left\{{\displaystyle \frac{1}{\nu ^2}}\left(1+2{\displaystyle \frac{\nu }{\eta }}\right)\left({\displaystyle \frac{\mathrm{\Psi }^{(0)}}{Y}}\right)^2+{\displaystyle \frac{1}{\eta ^2}}\left({\displaystyle \frac{\mathrm{\Phi }^{(0)}}{Y}}\right)^2\right\}`$ (8)   
$`+{\displaystyle \frac{1}{2\eta }}{\displaystyle \frac{^2}{Y^2}}\left({\displaystyle \frac{^2}{Y^2}}3{\displaystyle \frac{^2}{X^2}}\right)\mathrm{\Phi }^{(0)}+\nu ^4\mathrm{\Phi }^{(0)}.`$   
Let us focus our attention on diffusive terms in (6) and (8): as a consequence of the anisotropy of the basic small-scale flow the eddy-diffusivities are anisotropic too. For longitudinal perturbations ($`\mathrm{\Phi }^{(0)}=\mathrm{\Phi }^{(0)}(X,T)`$,$`\mathrm{\Psi }^{(0)}=\mathrm{\Psi }^{(0)}(X,T)`$) both viscosity and resistivity are left unchanged. On the other hand, for transverse perturbations ($`\mathrm{\Phi }^{(0)}=\mathrm{\Phi }^{(0)}(Y,T)`$,$`\mathrm{\Psi }^{(0)}=\mathrm{\Psi }^{(0)}(Y,T)`$) the renormalization of resistivity due to the small scale magnetic energy holds a negative term ($`1/2\nu `$), while molecular viscosity is increased by the eddy contribution ($`1/2\eta `$).

It is interesting to notice that the possibility of a negative eddy-resistivity due to small scale magnetic energy has already been presented in . That result, obtained in the framework of closure approaches to 2D MHD turbulence, was suggested as an explanation of the square magnetic potential inverse cascade. Analogous results are reported in , where the effective resistivity is shown to become negative in a small scale turbulent plasma, as long as the magnetic energy exceeds the kinetic one.

By inspection of equations (6) and (8), we stress that the large scale magnetic flux and stream function are linearly decoupled. If we assume large scale perturbations of the type $`\mathrm{\Psi }^{(0)}\mathrm{exp}(\mathrm{\Gamma }_\mathrm{\Psi }T+ıK_XX+ıK_YY)`$ and $`\mathrm{\Phi }^{(0)}\mathrm{exp}(\mathrm{\Gamma }_\mathrm{\Phi }T+ıK_XX+ıK_YY)`$, stability analysis leads to the following dispersion relationships

$`\mathrm{\Gamma }_\mathrm{\Psi }`$ $`=`$ $`\eta K_X^2(\eta {\displaystyle \frac{1}{2\nu }})K_Y^2,`$ (9)   
$`\mathrm{\Gamma }_\mathrm{\Phi }`$ $`=`$ $`{\displaystyle \frac{1}{K^2}}\left[(\nu +{\displaystyle \frac{1}{2\eta }})K_Y^4+2(\nu {\displaystyle \frac{3}{4\eta }})K_X^2K_Y^2+\nu K_X^4\right].`$ (10)   
The stability problem can be tackled more easily by introducing the parameters $`P=1/2\eta \nu `$ and $`T=(K_X/K_Y)^2`$. Marginal stability lines ($`\mathrm{\Gamma }_\mathrm{\Psi }=\mathrm{\Gamma }_\mathrm{\Phi }=0`$) are then given by

$`1P+T=0,`$   
$`1+P+2(1{\displaystyle \frac{3}{2}}P)T+T^2=0.`$   
and plotted in Figure 1.

For high enough values of molecular resistivity and viscosity ($`P<1`$), the basic flow is stable against any large scale perturbation. Increasing the Reynolds numbers, the first instability sets in at $`P=1`$ and $`T=0`$, i.e. for transverse perturbations. We notice that for $`1<P<16/9`$ the large scale vorticity is always stable ($`\mathrm{\Gamma }_\mathrm{\Phi }<0`$) and the magnetic potential growth rate $`\mathrm{\Gamma }_\mathrm{\Psi }`$ is maximum in the case $`T=0`$.

In the following, we investigate the linear and nonlinear evolution of large scale perturbations in the neighborhood of the critical point $`(P=1,T=0)`$, where we expect that there is no amplification of kinetic energy but that due to non-linear effects of coupling between velocity and magnetic field.

## IV Marginal instability

Special attention deserves the development of this large scale instability for Reynolds numbers close to the marginal stability threshold. In this regime, it is possible indeed to follow the full nonlinear evolution of the perturbation. We show that the instability eventually reaches a fixed point characterized by a magnetic island of the size of the box. The nonlinear evolution in the marginal regime is described by two coupled equations which generalize the Cahn-Hilliard equation, found for the hydrodynamical counterpart of this system, the so-called Kolmogorov flow .

Let us suppose to move the parameters just above the marginal stability line:

$$\eta =\eta _c(1\epsilon ^2),\nu =\nu _c(1\epsilon ^2)$$
(11)   
where $`\eta _c\nu _c=1/2`$ ($`P=1`$). The perturbative parameter $`\epsilon `$ is thus fixed by the distance between $`\eta ,\nu `$ and their critical values $`\eta _c,\nu _c`$. We will take into account only transverse perturbations, since, as shown in Figure 1, the large scale magnetic flux linear instability is mainly transverse, close to $`P=1`$. According to (6), the transverse eddy-resistivity in the neighborhood of the critical line defined by (11) is of order $`O(\epsilon ^2)`$, thus suggesting a scaling for the slow time $`T=\epsilon ^4t`$. The decomposition rules (3) become

$$_x_x_y_y+\epsilon _Y,_t_t+\epsilon ^4_T.$$
(12)   
The same multiscale technique described above can be adopted to solve perturbatively (1) and (2).

At first order in $`\epsilon `$ we obtain:

$`\psi `$ $`=`$ $`\mathrm{cos}x+\mathrm{\Psi }^{(0)}(Y,T)+\epsilon \mathrm{\Psi }^{(1)}(Y,T)`$ (13)   
$`\phi `$ $`=`$ $`2\eta _c\epsilon {\displaystyle \frac{\mathrm{\Psi }^{(0)}}{Y}}\mathrm{sin}x.`$ (14)   
We notice that, according to the conclusions drawn above, the large scale stream function is linearly stable and it is simply driven by the magnetic flux.

The evolution equations for $`\mathrm{\Psi }^{(0)}`$ and $`\mathrm{\Psi }^{(1)}`$ emerge as solvability conditions at order $`\epsilon ^4`$ and $`\epsilon ^5`$:

$`_T\mathrm{\Psi }^{(0)}`$ $`=`$ $`{\displaystyle \frac{27}{8}}\eta _c_{4Y}\mathrm{\Psi }^{(0)}2\eta _c_{YY}\mathrm{\Psi }^{(0)}+12\eta _c_{YY}\mathrm{\Psi }^{(0)}(_Y\mathrm{\Psi }^{(0)})^2`$ (15)   
$`_T\mathrm{\Psi }^{(1)}`$ $`=`$ $`{\displaystyle \frac{27}{8}}\eta _c_{4Y}\mathrm{\Psi }^{(1)}2\eta _c_{YY}\mathrm{\Psi }^{(1)}+12\eta _c_{YY}\mathrm{\Psi }^{(1)}(_Y\mathrm{\Psi }^{(0)})^2`$ (17)   
$`+24\eta _c_Y\mathrm{\Psi }^{(0)}_{YY}\mathrm{\Psi }^{(0)}_Y\mathrm{\Psi }^{(1)}`$   
In the first equation (15), one easily recognizes the renowned Cahn-Hilliard equation , which may be written in variational form:

$$\frac{\mathrm{\Psi }^{(0)}}{t}=\frac{\delta V\left[\mathrm{\Psi }^{(0)}\right]}{\delta \mathrm{\Psi }^{(0)}}.$$

The existence of the Lyapunov functional

$$V\left[\mathrm{\Psi }^{(0)}\right]=\eta _c𝑑Y\left[(_Y\mathrm{\Psi }^{(0)})^2+(_Y\mathrm{\Psi }^{(0)})^4+\frac{27}{16}(_{YY}\mathrm{\Psi }^{(0)})^2\right]$$
(18)   
indicates that asymptotically the solution of (15) in a bounded domain reaches a fixed point. This stationary solution is approached by a sequence of metastable states of decreasing dominating mode. We thus expect to observe a nonlinear evolution dominated by a magnetic island coalescence, analogous to the vortex pairing in 2D hydrodynamics .

Equation (17) is linear in $`\mathrm{\Psi }^{(1)}`$, with coefficients depending nonlinearly on $`\mathrm{\Psi }^{(0)}`$. It also can be written as a gradient flow, with a Lyapunov functional

$$V\left[\mathrm{\Psi }^{(1)}\right]=\eta _c𝑑Y\left[(_Y\mathrm{\Psi }^{(1)})^2+6(_Y\mathrm{\Psi }^{(0)})^2(_Y\mathrm{\Psi }^{(1)})^2+\frac{27}{16}(_{YY}\mathrm{\Psi }^{(1)})^2\right].$$
(19)   
We conclude this section by observing that the dispersion relation for $`\mathrm{\Psi }^{(0)}\mathrm{exp}(\mathrm{\Gamma }T+ıKY)`$ now reads:

$$\mathrm{\Gamma }=\frac{27}{8}\eta _cK^4+2\eta _cK^2.$$
(20)   
It implies instability ($`\mathrm{\Gamma }<0`$) for any ($`K\stackrel{}{<}0.77`$). We notice that information about the characteristic scale of unstable modes were totally absent in the general treatment presented in the previous section (see equation (9).

## V Numerical results

The results obtained in the previous section have been checked by extensive direct numerical simulations of MHD equations (1,2). In order to force a transverse perturbation, we integrate the equations on a rectangular slab with $`L_{box}=L_x=2\pi `$ and $`L_yL_x`$. In this way large scale instability can only develop on the $`y`$ direction for an aspect ratio $`r=L_x/L_y<1`$.

Given the numerical values of parameters $`\eta `$ and $`\nu `$, from (11) and the condition $`\eta _c\nu _c=1/2`$ we have

$$\epsilon =\sqrt{1\sqrt{2\eta \nu }},\eta _c=\sqrt{\frac{\eta }{2\nu }}$$
(21)   
which are used for the theoretical predictions of the previous section.

The simplest check of our predictions concerns the growth rates of the instability which, in the initial linear regime, are given by the dispersion relation (20). In physical (not rescaled) variables (20) becomes

$$\gamma =\frac{27}{8}\eta _ck^4+2\eta _c\epsilon ^2k^2$$
(22)   
which shows that the largest unstable wavenumber is $`k_{max}0.77\epsilon `$. The smallest transverse wavenumber is $`k_1=r`$ thus, in order to numerically observe the instability, it must be $`r0.77\epsilon `$.

In Figure 2 we report the growth rates of the first modes for a simulation with $`r=1/64`$, $`\nu =0.49`$ and $`\eta =1`$. We have $`\epsilon 0.1`$ and thus only the first $`4`$ modes are unstable. The initial perturbation is small, random and on all the first $`20`$ modes, thus we are able to observe also negative $`\gamma `$’s (stable modes). The comparison with the linear prediction (22) is very good even for $`\epsilon `$ not very small. The numerical data of Figure 2 are obtained by a linear fit of the logarithm of the mode amplitude versus time in the early stages of the simulation.

Let us now consider the nonlinear stage of the perturbation growth. We describe here a different simulation with $`r=1/16`$ and $`\epsilon 0.32`$ which was advanced for a very long lapse of time. The nonlinear evolution will ultimately lead to a fixed point by a succession of long lasting quasi equilibrium states of decreasing wavenumber. The evolution of the amplitudes of the first $`5`$ transverse modes computed from the direct numerical simulation is plotted in Figure 3. Observe that in this case the fifth mode $`k=5/16`$ is linearly stable but it is non-linearly driven by smaller wavenumber.

The typical linear time is now rather short, $`1/\gamma O(1)`$, and the final stationary state, dominated by the largest available mode $`k_1`$, is reached at very long times, $`t>1000`$. At intermediate times, almost stationary metastable states, characterized by decreasing leading mode, are punctuated by fast coalescence processes. Most of the energy dissipation takes place during this fast reconnection processes. In Figure 4 we display the period-two metastable state at time $`t=200`$ and the final, period-one state at $`t=20000`$. The dynamical picture arising from Figures 3 and 4 qualitatively agrees with the dynamics described by the Cahn-Hillard equation .

To check quantitatively the validity of non linear multiscale analysis we have numerically integrated the Cahn-Hilliard equation for the large scale magnetic flux (15) with the same parameters of the DNS. As shown in Figure 5 we find an impressive agreement even for very long times. The final relative amplitudes of the most energetic transverse modes is recovered within a $`10\%`$ accuracy.

As a further test of the multiscale predictions, we checked the relations (13-14) during the evolution. At leading order in $`\epsilon `$, $`\mathrm{\Psi }^{(0)}(y,t)`$ is obtained by subtracting the basic flow $`\mathrm{cos}x`$ from the magnetic flux $`\psi (x,y,t)`$. The resulting field, which reveals to be indeed $`x`$-independent, is then used to reconstruct the stream function by means of (14). The results for the configuration of Figure 4 is shown in Figure 6.

## VI Conclusions

We have investigated the issue of stability of highly anisotropic, magnetic-energy dominated equilibrium states of the MHD fluid model equations. These configurations are known to be unstable for small values of resistivity leading to the formation of thin boundary layers in the nearby of the neutral line of the magnetic field. At variance with the above case, we focused our attention on the range of moderate Lundquist/Reynolds numbers, where the boundary layer approximation is not fruitfully applicable. In this situation we have shown analytically, by means of multiple-scale analysis, that large scale instabilities can arise above a definite threshold, and that for a generic perturbation the maximum growth is achieved by modes transverse to the magnetic field lines of the basic state. On the basis of this result, we have performed the multiple scale analysis for the marginally unstable case and for purely transverse perturbations. The analytic procedure yields a couple of partial differential equations which describe the full nonlinear evolution of magnetic perturbations. It is possible to show that these equations possess a Lyapunov functional and thus their solutions asymptotically approach a fixed point which represents a nonlinear equilibrium different from the basic one. The kinetic perturbations are always linearly stable, and their growth is uniquely due to the nonlinear coupling to the magnetic field. Numerical simulations of two-dimensional MHD performed with a pseudospectral code reveal an excellent quantitative agreement with the first-order analytical prediction in a wide range of values of the perturbative parameter.

We therefore conclude that the loss of stability of parallel magnetic field configurations at moderate Reynolds number is due to the growth of large scale perturbations, and that the features of this instability can be captured by a multiple scale analysis. The transverse large-scale instability is likely to be the generic mechanism of instability of sheared magnetic fields even for large Lundquist/Reynolds numbers whenever the basic state admits a large number of neutral lines. When there is a single neutral line at large enough Lundquist/Reynolds numbers this mechanism is overcome by the formation of resistive boundary layers.

## VII Acknowledgments

The authors wish to acknowledge support and hospitality by the Istituto di Cosmogeofisica CNR, Torino. The calculations were partially performed with computer facilities of INFN, Sezione di Torino.