# Explosive Magnetic Reconnection in Double-current Sheet Systems: Ideal versus Resistive Tearing Mode

Hubert Baty

► To cite this version:

Hubert Baty. Explosive Magnetic Reconnection in Double-current Sheet Systems: Ideal versus Resistive Tearing Mode. The Astrophysical Journal, 2017, 837 (1), pp.74. 10.3847/1538-4357/aa60bd hal-03143024

HAL Id: hal-03143024

https://hal.science/hal-03143024v1

Submitted on 23 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Explosive Magnetic Reconnection in Double-current Sheet Systems: Ideal versus Resistive Tearing Mode

Hubert Baty

Observatoire Astronomique de Strasbourg, Université de Strasbourg, CNRS UMR 7550, F-67000 Strasbourg, France; hubert.baty@unistra.fr

Received 2016 December 12; revised 2017 February 13; accepted 2017 February 13; published 2017 March 6

## Abstract

Magnetic reconnection associated with the tearing instability occurring in double-current sheet systems is investigated within the framework of resistive magnetohydrodynamics (MHD) in a two-dimensional Cartesian geometry. A special emphasis on the existence of fast and explosive phases is taken. First, we extend the recent theory on the ideal tearing mode of a single-current sheet to a double-current layer configuration. A linear stability analysis shows that, in long and thin systems with (length to shear layer thickness) aspect ratios scaling as $S_L^{9/29}$ ($S_L$ being the Lundquist number based on the length scale $L$), tearing modes can develop on a fast Alfvénic timescale in the asymptotic limit $S_L \to \infty$. The linear results are confirmed by means of compressible resistive MHD simulations at relatively high $S_L$ values (up to $3 \times 10^6$) for different current sheet separations. Moreover, the nonlinear evolution of the ideal double tearing mode (IDTM) exhibits a richer dynamical behavior than its single-tearing counterpart, as a nonlinear explosive growth violently ends up with a disruption when the two current layers interact trough the merging of plasmons. The final outcome of the system is a relaxation toward a new state, free of magnetic field reversal. The IDTM dynamics is also compared to the resistive double tearing mode dynamics, which develops in similar systems with smaller aspect ratios, $\gtrsim 2\pi$, and exhibits an explosive secondary reconnection, following an initial slow resistive growth phase. Finally, our results are used to discuss the flaring activity in astrophysical magnetically dominated plasmas, with a particular emphasis on pulsar systems.

**Key words:** magnetic reconnection – magnetohydrodynamics (MHD) – plasmas – pulsars: general

## 1. Introduction

Magnetic reconnection is a fundamental process in laboratory and space plasmas that allows the conversion of magnetic energy into bulk flow and heating. Understanding magnetic reconnection is a major goal of theoretical plasma physics in order to explain explosive events like disruptions in fusion laboratory experiments and flares in astrophysical plasmas. For such collisional plasmas, the magnetohydrodynamic (MHD) model is commonly adopted as an excellent framework (Priest & Forbes 2000). The MHD approximation is a single-fluid description of plasma dynamics, where the electromagnetic properties are taken into account via coupling terms obeying Maxwell's equations. Typically, electrical currents can be generated via fluid motions producing Lorentz forces on the fluid. The plasma resistivity (inverse of the electrical conductivity) $\eta$ representing the magnetic field dissipation is very small in plasmas of interest, but it cannot be neglected, as it precisely drives the reconnection process. Magnetic reconnection timescales are traditionally expressed by using some nondimensional Lundquist number $S$, where $S = L_0 v_A / \eta$. $L_0$ and $v_A$ are the characteristic macroscopic length and Alfvèn velocity of the system, respectively.

Static current sheets with aspect ratios $L/a \gtrsim 2\pi$ ($L$ and $a$ being characteristic length and thickness, respectively) are prone to spontaneous reconnection, due to the development of the tearing instability. Tearing modes have been extensively studied, as they are believed to play a central role in the reconnection phenomenon. It is commonly admitted that these modes are slow instabilities growing on a characteristic time scaling as $S_a^{1/2}$, with $S_a = \alpha v_A / \eta$. This gives the well known Sweet-Parker (SP) time defined as $\tau_{SP} = (a/v_A) S_a^{1/2}$, leading to a timescale intermediate (geometric mean) between the fast Alfvèn time $\tau_A = a/v_A$ and the extremely slow diffusion time

$\tau_d = a^2/\eta$, as $\tau_{SP} = (\tau_d \tau_A)^{1/2}$. Consequently, as $S_a = 10^{8-12}$ in most plasmas of interest, tearing reconnection cannot account for fast timescales involved in flares or tokamak disruptions.

However, it was recently discovered that current sheets with relatively large aspect ratios $L/a$, scaling as $S_L^{1/3}$ (with $S_L = L v_A / \eta$), can disrupt on the rapid Alfvèn timescale $t_A = L/v_A$ provided $S_L \gtrsim 10^6$, thus independently of the resistivity (Landi et al. 2015; Tenerani et al. 2015). These studies have confirmed the existence of a new tearing mode developing with a linear growth rate $\gamma t_A \approx 0.62$, which has been discovered using a revisited stability analysis (Pucci & Velli 2014). The latter study is based on a brilliant formulation of renormalization in which $L$ defines the macroscopic length instead of $a$ (which is traditionally used in tearing mode theories), the aspect ratio $L/a$ thus becoming a free parameter. This theory of the so-called *ideal* tearing also solves the paradox of the linear growth rate $\gamma$ of super-Alfvénic tearing instabilities (also called plasmons) developing in SP current sheets (for which $L/a \sim S_L^{1/2}$), with a positive power of $S_L$, i.e., $\gamma t_A \sim S_L^{1/4}$, thus diverging in the ideal MHD limit (Samtaney et al. 2009; Huang & Bhattacharjee 2010; Loureiro et al. 2012).

The aim of the present study is twofold. First, we extend the linear stability analysis of the ideal tearing instability taking place in single-current sheets to double-current layer systems, where two identical layers of oppositely directed current density are separated by a distance $2l$ (see Figure 1). In this work, we focus on a range of $l/a$ values varying between approximately unity and 10. Linear and nonlinear phases of *ideal* tearing developing in relatively long and thin configurations are thus investigated by means of resistive MHD simulations. Second, similar systems with smaller aspect ratios, $L/a \gtrsim 2\pi$, have been previously shown to lead to explosive
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Schematic view of the initial double-current sheet configuration. The magnetic field is unidirectional, with a component $B_y(x)$ experiencing two reversals and forming two current layers of negative and positive current density of half-thickness *a* (blue and red regions). The two current layers of length $2L$ are separated by a half-distance *l*.</figcaption></figure>

and fast nonlinear reconnection phases (roughly independent of the resistivity), which follows a relatively long initial (linear and early nonlinear) phase characteristic of the double resistive tearing mode (i.e., on a slow resistive timescale), as long as an initial critical threshold $2L/l \gtrsim 5$ is taken (Janvier et al. 2011). Our aim is therefore to compare the dynamics (fastness and explosivity) of the two mechanisms, namely ideal versus resistive tearing in double-current sheet systems.

The paper is organized as follows. Section 2 presents a linear stability study for an idealized double-current sheet system, using the previously mentioned renormalization technique (Pucci & Velli 2014). Numerical simulations are presented in Section 3. Finally, conclusions and discussions are reported in Section 4.

## 2. The Double Tearing Mode: A Linear Analysis

In order to carry out a simple linear analysis, we start from an idealized double-current sheet system where the equilibrium magnetic field varies linearly inside the two current layers and is constant elsewhere (Otto & Birk 1992). The corresponding current density is thus constant (step function profile) inside the two anti-parallel sheets (see Figure 1). The system is characterized by different spatial scales with *L* as the half-length, *a* as the half-thickness, and *l* as the half-separation of the sheets. Such a simple configuration allows us to obtain analytical dispersion relations that can be easily solved by some numerical root-finding technique (Birk et al. 1997). One must note that, according to the symmetry properties of the perturbation, two different types of unstable normal mode can be distinguished. Indeed, a varicose-like deformation of the inner region between the two layers can be associated with a symmetric mode, and a kink-like deformation corresponds to an anti-symmetric mode. In the present work, we focus only on the anti-symmetric mode, as its symmetric counterpart is clearly dominated for the range of current sheet separation considered here (and it is even stabilized for $l \approx a$) (Birk et al. 1997; Shen & Liu 1998).

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Linear growth rate $\gamma\tau_A$ as a function of $ka$, obtained by solving a dispersion relation for a simplified double-current layer equilibrium (see the text) for a Lundquist number $S_a = 10^6$. The different curves correspond to different relative distances between the two layers, ranging between $l/a = 1.25$ and $l/a = 250$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Maximum linear growth rate $\gamma_m\tau_A$ and corresponding wavenumber $k_m a$ as a function of the Lundquist number $S_a$, obtained for a high $l/a$ value ($l/a = 250$). The calculated values (square points and asterisks) follow the expected scalings for the single-tearing mode in the $S_a \to \infty$ limit (see Equations (1) and (2) in the text).</figcaption></figure>

### 2.1. The Resistive Double Tearing Mode (RDTM)

We numerically solve the dispersion relation derived for the anti-symmetric mode in a double-current sheet system and given by Equation (5) in Birk et al. (1997).

The resulting linear growth rate $\gamma\tau_A$ is plotted in Figure 2 as a function of $ka$ for five values of the relative sheet separation, namely $l/a = 1.25, 2.5, 5, 10,$ and $250$. The curves are obtained for a relatively high fixed Lundquist number value $S_a = 10^6$. Note that only modes with a wavenumber $ka \lesssim 0.75$ are unstable, differing from the $ka \lesssim 1$ expected from classical resistive tearing theory. This difference is simply due to the linear magnetic reversal profile assumed in the sheet, which is used instead of the traditional hyperbolic-like profile.

We have checked that the resistive single-mode theory is recovered for $l/a = 250$, as required in the limit $l/a \to \infty$. Indeed, two single-tearing instabilities are expected to develop independently on the two layers in this limit. More precisely, the measured maximum growth rate obtained for $S_a = 10^6$ (and
<figure><img src="image_6.png" /><figcaption>Figure 4. Same as the previous figure obtained for four smaller $l/a$ values. The results in the high $S_a$ limit can be approximated by scaling laws that differ from the single-tearing mode, as $\gamma_m \tau_A \approx \alpha \times S_a^{-9/20}$ and $k_m a \approx \beta \times S_a^{-3/20}$ for the linearly fastest growing mode and corresponding wavenumber, respectively. The values of the two constants $\alpha$ and $\beta$ depend on $l/a$ (see the text).</figcaption></figure>

deduced from Figure 2) is $\gamma_m \tau_A = 6.06 \times 10^{-2}$ for the dominant wavenumber $k_m a = 3.95 \times 10^{-2}$. These values agree well with the values of $6.2 \times 10^{-2}$ and $4.2 \times 10^{-2}$ deduced from the theoretical formula (Priest & Forbes 2000),

$$ \gamma_m \tau_A \approx 0.62 \times S_a^{-1/2}, \qquad (1) $$

and

$$ k_m a \approx 1.4 \times S_a^{-1/4}, \qquad (2) $$

for the maximum growth rate and corresponding linearly fastest growing wavenumber, respectively.

The closeness effect of the two current layers is clearly to enhance the maximum growth rate by a factor of the order of two when comparing our two extreme $l/a$ cases, in agreement with previously reported results using hyperbolic tangent profiles for the equilibrium (Shen & Liu 1998). A similar effect is also observed for the corresponding wavenumber $k_m$, which is shifted toward higher values as $l/a$ decreases.

We have investigated the dependence of $\gamma_m$ and $k_m$ on the Lundquist number for $S_a$ ranging from $10^2$ to $10^8$, in order to deduce some scaling law as $S_a \to \infty$. The results for $l/a = 250$ are plotted in Figure 3, while Figures 4(a)–(d) displays the other $l/a$ cases. First, note that for $l/a = 250$, scaling laws of $\gamma_m \tau_A \approx 0.62 \times S_a^{-1/2}$ and $k_m a \approx 1.4 \times S_a^{-1/4}$ are obtained in the limit $S_a \to \infty$, which compares very well with the expected theoretical expressions (given by Equations (1) and (2)), in spite

<figure><img src="image_7.png" /></figure>

of the use of our simplified magnetic equilibrium profile. The effect of the two layers’ closeness, with $l/a$ varying in the range between 1.5 and 10, can now be explored as discussed below. First, the asymptotic scaling law for the maximum linear growth rate is only slightly modified compared to the highest $l/a$ case. Indeed, when $S_a \to \infty$, our results follow a scaling formula,

$$ \gamma_m \tau_A \approx \alpha S_a^{-9/20}, \qquad (3) $$

where $\alpha$ is a parameter of order unity. The exact value of $\alpha$ is varying between 0.5 and 0.8 when $l/a$ decreases from 10 to 1.25. The $S_a$ exponent of $-9/20 = -0.45$ for our different $l/a$ cases remains thus very close to the $-0.5$ value representative of the single-tearing. However, this is different for the corresponding wavenumber, which is shown to follow a $S_a^{-0.15}$ dependency instead of a $S_a^{-1/4}$ one. More precisely, we obtain,

$$ k_m a \approx \beta S_a^{-3/20}, \qquad (4) $$

with $\beta$ being a parameter depending on the $l/a$ ratio value. Typically, $\beta$ is of order unity and increases from 0.35 to 0.77 when $l/a$ decreases from 10 to 1.25. These scaling laws roughly agree with previous results for small enough separation of the current sheets, where scaling of $S_a^{-3/7}$ and $S_a^{-1/7}$ were reported for the maximum growth rate and corresponding wavenumber, respectively (Otto & Birk 1992).
<table><caption>**Table 1**
Parameters Used for the Simulations</caption><thead><tr><th>Run n</th><th>l/a</th><th>S<sub>a</sub></th><th>S<sub>L</sub></th><th>L/a</th><th>S<sub>L</sub><sup>9/29</sup></th><th>S<sub>L</sub><sup>1/2</sup></th><th>γ<sub>th</sub>t<sub>*</sub></th><th>γ<sub>nu</sub>t<sub>*</sub></th><th>k<sub>th</sub>a</th><th>k<sub>nu</sub>a</th></tr></thead><tbody><tr><td>Run 1</td><td>5</td><td>1429</td><td>1 × 10<sup>5</sup></td><td>70</td><td>35</td><td>317</td><td>0.418</td><td>0.36</td><td>0.151</td><td>0.11</td></tr><tr><td>Run 2</td><td>5</td><td>3000</td><td>3 × 10<sup>5</sup></td><td>100</td><td>50</td><td>550</td><td>0.302</td><td>0.28</td><td>0.135</td><td>0.10</td></tr><tr><td>Run 3</td><td>5</td><td>6667</td><td>1 × 10<sup>6</sup></td><td>150</td><td>73</td><td>1000</td><td>0.209</td><td>0.22</td><td>0.120</td><td>0.12</td></tr><tr><td>Run 4</td><td>5</td><td>13636</td><td>3 × 10<sup>6</sup></td><td>220</td><td>102</td><td>1730</td><td>0.152</td><td>0.17</td><td>0.108</td><td>0.10</td></tr><tr><td>Run 5</td><td>2.5</td><td>13636</td><td>3 × 10<sup>6</sup></td><td>220</td><td>102</td><td>1730</td><td>0.184</td><td>0.21</td><td>0.139</td><td>0.12</td></tr><tr><td>Run 6</td><td>10</td><td>13636</td><td>3 × 10<sup>6</sup></td><td>220</td><td>102</td><td>1730</td><td>0.139</td><td>0.16</td><td>0.084</td><td>0.09</td></tr><tr><td>Run 7</td><td>2.5</td><td>1.4 × 10<sup>4</sup></td><td>1.4 × 10<sup>5</sup></td><td>20</td><td>...</td><td>...</td><td>0.085</td><td>0.08</td><td>...</td><td>...</td></tr><tr><td>Run 8</td><td>2.5</td><td>5.0 × 10<sup>4</sup></td><td>5.0 × 10<sup>5</sup></td><td>20</td><td>...</td><td>...</td><td>0.042</td><td>0.04</td><td>...</td><td>...</td></tr></tbody></table>

Note. Note that $a = 0.05$ is fixed for all the cases. In the first set (runs 1–4), IDTM is explored with different $S_L$ values for a fixed current sheet separation $l/a = 5$. In a second set (runs 5 and 6), IDTM is explored with different $l/a$ values for a fixed Lundquist number $S_L = 3 \times 10^6$. The third set of simulations (run 7 and 8) explores the RDTM for two Lundquist numbers. The aspect ratios chosen for the runs are reported in the fourth column and can be compared to the critical and Sweet-Parker ratios in the fifth and sixth columns. Estimated maximum linear growth rates (from simulations) $\gamma_{nu}$ and corresponding wavenumbers $k_{nu}$ are compared to the expected theoretical values ($\gamma_{th}$ and $k_{th}$) using Equations (3) and (4) in the last columns.

## 2.2. The Ideal Double Tearing Mode (IDTM)

As explained in the introduction, it has recently been realized that single-tearing modes can evolve on a fast Alfvénic timescale, provided the aspect ratio of the current layer $L/a$ reaches a critical value of $S_L^{1/3}$ and $S_L \to \infty$ (Pucci & Velli 2014). More precisely, for configurations with the critical aspect ratio value, the growth rate of the tearing instability becomes independent of the resistivity already at $S_L \gtrsim 10^6$. In this high Lundquist number regime, the resulting linear growth rate is $\gamma_{tA} \approx 0.62$ (thus of order unity) with $t_A = L/v_A$, leading the authors to introduce the concept of an *ideal* tearing mode. The latter result has been obtained by a renormalization technique ($L$ is used instead of $a$ in the definition of the Alfvèn time). The ideal linear growth rate corresponds to the linear growth rate of the dominant tearing mode in classical resistive theory (Equation (1)), which becomes independent of $S_L$ for current layer aspect ratios equal to the critical value of $S_L^{1/3}$. The IDTM theory has been successfully checked in a series of papers (Landi et al. 2015; Del Sarto et al. 2016).

Following the same renormalization procedure, we can now investigate the existence of an IDTM. Indeed, the maximum growth rate scaling law (for $S_L \to \infty$) deduced in the previous subsection can be written as $\gamma_m t_A \approx \alpha S_L^{-9/20} (L/a)^{29/20}$, leading to $\gamma_m t_A \approx \alpha$ that is independent of $S_L$ for a current sheet with an aspect ratio $L/a = S_L^{9/29}$. This is only slightly different from the critical value obtained for the single ideal tearing aspect ratio of $L/a = S^{1/3}$. Note also that the maximum growth rate value varies slightly with $l/a$ (as $\alpha = 0.5 - 0.8$), and consequently remains close to the single ideal tearing value of 0.62. Concerning the corresponding wavenumber, our result of $k_m a \approx \beta S_a^{-3/20}$ is transformed into $k_m a \approx \beta S_L^{-3/20} (L/a)^{3/20}$, leading to

$$k_m L \approx \beta S_L^{6/29}, \quad (5)$$

when the critical aspect ratio $L/a = S_L^{9/29}$ is taken. This latter result differs substantially from the maximum wavenumber for the single ideal tearing mode for which $k_m L \approx 1.4 \times S_L^{1/6}$ (Del Sarto et al. 2016).

## 3. Numerical Simulations: IDTM versus RDTM

In order to deduce a characteristic timescale representative of the full time evolution of the configuration, it is also necessary

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Maximum absolute value of the magnetic field component $B_x$ measured over the computational domain as a function of time for run 4. An exponential dependence in $\exp(0.17t/t_*)$ is plotted in plain line in order to fit the linear phase.</figcaption></figure>

to explore the nonlinear development of the double tearing instability. In previous studies focusing on the single ideal tearing mode, resistive MHD simulations confirmed that current sheets with a critical aspect ratio $L/a = S_L^{1/3}$ are able to develop modes on the ideal Alfvénic time, with an initial growth rate as fast as $\gamma \approx 0.62 t_A^{-1}$, in agreement with linear stability analysis (Landi et al. 2015; Tenerani et al. 2015). The unstable perturbations take the form of the eigenmodes expected from the linear theory and manifest themselves as small magnetic islands (plasmoids) growing and later moving along the layer.

Beyond the linear evolution dominated by the linearly fastest mode, a double cascade process is observed. First, new secondary current sheets (of half-length $L^*$ and half-width $a^*$) can locally form on a smaller scale when local critical critical aspect ratios of $L^*/a^* = S_L^{1/3}$ are reached. This leads to new secondary plasmoids developing on a super fast Alfvénic timescale $t_A^* \approx L^*/V_A$. An explosive reconnection thus ensues as faster and faster reconnection produces a cascading accelerating process toward small scales. Second, an inverse process toward large scales takes place at the same time, as the plasmoids are able to merge, leading to the formation of larger plasmoids.
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Left panel: contour map of $B_x$ for run 1 at $t = 6t_*$ Right panel: corresponding $B_x(x)$ dependence obtained with a transverse cut at $y = 3.2$ through the left current layer.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Left panel: same as Figure 5 but obtained for runs 1–4 showing the $S_L$ effect for configurations using $L/a = 5$. Right panel: same as Figure 5 but for runs 4–6 showing the $l/a$ effect for configurations using $S_L = 3 \times 10^6$.</figcaption></figure>

In this work, we carry out nonlinear MHD simulations in order to explore such nonlinear development in the context of double-current sheet systems with aspect ratios of the order of the critical values (obtained in the previous section). As the critical aspect ratio scaling as $S_L^{9/29}$ is rather large for typical high Lundquist number plasmas, the configurations represent long and thin layers but not as much as SP ones, which scale as $S_L^{1/2}$. For example, this lead to a critical aspect ratio value of 102, which is one order of magnitude smaller compared to the SP ratio of 1732, when one takes the best Lundquist value employed in our simulations, i.e., $S_L = 3 \times 10^6$.

On the other hand, in much smaller aspect ratio systems with $L/a \gtrsim 2\pi$ RDTM have been shown to nonlinearly develop

with an explosive and fast phase (Wang et al. 2007; Janvier et al. 2011; Zhang & Ma 2011). Typically, the long resistivity-dependent linear and subsequent early nonlinear phases are followed by a third explosive evolution phase that scales quasi-independently with the resistivity, namely as $\eta^{0.06}$ (Zhang & Ma 2011). A secondary ideal nonlinear instability that becomes structurally unstable, provided $2L/l \gtrsim 5$, has been proposed (Janvier et al. 2011), and was later confirmed by other studies (Nemati et al. 2016). We thus propose to compare the explosiveness and timescale involved in RDTM for small aspect ratio configurations with an IDTM characterizing the evolution of much longer systems, by means of nonlinear MHD simulations.
<figure><img src="image_11.png" /><figcaption>Figure 8. Contour map of the thermal pressure for run 3 at different times. Only a part with x in the range [-1:1] and y in the range [-2:12.5] of the whole simulation domain is shown.</figcaption></figure>

### 3.1. Model and Numerical Setup

We solve the standard set of compressible resistive MHD equations (viscosity and thermal conductivity are ignored) written in the following dimensionless form:

$$ \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{v}) = 0, \qquad (6) $$

$$ \frac{\partial \rho \mathbf{v}}{\partial t} + \nabla \cdot (\rho \mathbf{v} \mathbf{v}^T + p_{\text{tot}} \mathbf{I} - \mathbf{B} \mathbf{B}^T) = 0, \qquad (7) $$

$$ \frac{\partial \mathbf{B}}{\partial t} + \nabla \cdot (\mathbf{v} \mathbf{B}^T - \mathbf{B} \mathbf{v}^T) = -\nabla \times (\eta \mathbf{J}), \qquad (8) $$

$$ \frac{\partial e}{\partial t} + \nabla \cdot (e \mathbf{v} + p_{\text{tot}} \mathbf{v} - \mathbf{v} \cdot \mathbf{B} \mathbf{B}^T) = \eta \mathbf{J}^2 - \mathbf{B} \cdot \nabla \times (\eta \mathbf{J}), \qquad (9) $$

$$ \mathbf{J} = \nabla \times \mathbf{B}, \qquad (10) $$

$$ \nabla \cdot \mathbf{B} = 0. \qquad (11) $$

Here, $\rho$ is the plasma density, $p$ is the thermal pressure, $\mathbf{v}$ is the fluid velocity, $\mathbf{B}$ is the magnetic fluid, and $\mathbf{J}$ is the electrical current density. The total pressure (thermal + magnetic) is defined as $p_{\text{tot}} = p + B^2/2$, and the total energy density is $e = \rho \mathbf{v}^2/2 + B^2/2 + p/(\Gamma - 1)$, where $\Gamma$ is the ratio of specific heats ($\Gamma = 5/3$). The resistivity coefficient is $\eta$. Note that the magnetic permeability is taken to be unity.

We assume a double Harris-like current sheet configuration, with a magnetic field parallel to the y-axis and varying with $x$,

$$ B_x = 0, \quad B_y = B_e \left[ 1 + \tanh\left[ \frac{(x-l)}{a} \right] - \tanh\left[ \frac{(x+l)}{a} \right] \right], \qquad (12) $$

where $B_e$ is the amplitude of the field and $a$ is the half-width of each current layer. We set $B_e = 1$, and $a = 0.05$, to define our normalization. Consequently, our normalized time is $t_* = 20a/V_A$, where $V_A$ is the Alfvén velocity taken to be unity in our units. The two current layers are localized at $x = \pm l$. A static equilibrium is considered, with force balance provided by the plasma density variation in an isothermal medium. Note that the plasma-$\beta$ taken at the outer x-boundary is 0.5.

In most of this work, a total of 1920 × 480 uniformly distributed spatial grid points is used in the spatial domain $-L_x \le x \le L_x$ and $0 \le y \le L_y$, of fixed dimension $L_x = 2$ and variable one $L_y$. In a few cases, we have also checked that doubling this spatial resolution does not change the overall scenario described below and our conclusions. Small differences during the very late evolution phase when relaxation occurs are unavoidable, due to turbulence inherent in using high Lundquist numbers (Keppens et al. 2013). This resolution is also similar to the one previously employed to numerically investigate the ideal single-tearing evolution (Landi et al. 2015). Fixed boundary conditions are chosen at $x = \pm L_x$ that are placed sufficiently far enough away from the double layer in order to have a negligible effect on the dynamics. Periodic boundary conditions are assumed in the y direction.

We use the general finite-volume based Versatile Advection Code MPI-AMRVAC (Porth et al. 2014). The MHD equations are solved by selecting a second-order shock-capturing scheme, namely the HLL solver combined with a Koren flux limiter.
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Maximum absolute value of the magnetic field component $B_x$ (left panel) and velocity field component $V_y$ (right panel) measured over the computational domain as a function of time for runs 7 and 8 (RDTM) vs. run 5 (IDTM).</figcaption></figure>

The numerically generated nonzero divergence of the magnetic field is handled by choosing the diffusive approach.

In order to trigger the instability, the initial equilibrium configuration is perturbed by adding divergence-free components of relatively small amplitude $\epsilon \approx 10^{-3}$,

$$ \delta B_x = \sum_{n=1}^{10} \epsilon B_e \operatorname{sech}[(x-l)/a] \cos(k_n^* y/a + \phi_n), \quad (13) $$

and

$$ \delta B_y = \sum_{n=1}^{10} \epsilon B_e / k_n^* \operatorname{sech}[(x-l)/a] \tanh[(x-l)/a] \times \sin(k_n^* y/a + \phi_n), \quad (14) $$

if $x > 0$. And if $x < 0$, the perturbation follows,

$$ \delta B_x = \sum_{n=1}^{10} \epsilon B_e \operatorname{sech}[(x+l)/a] \cos(k_n^* y/a + \phi_n), \quad (15) $$

and

$$ \delta B_y = \sum_{n=1}^{10} \epsilon B_e / k_n^* \operatorname{sech}[(x+l)/a] \tanh[(x+l)/a] \times \sin(k_n^* y/a + \phi_n), \quad (16) $$

where $k_n^* = 2\pi n a / L_y$ and $\phi_n$ are selected (normalized) wavenumber and random phase, respectively, for a given $n$ mode. This above form with sinusoidal shape in $y$ and concentrated shape around $x = \pm l$ is suitable, and is similar to the one used to investigate the relativistic single-tearing instability (Del Zanna et al. 2016). Note that we excite 10 modes along the configuration because the linearly fastest mode (according to the linear theory) can vary approximately from $n = 1$ to $n = 8$, depending on the run considered.

## 3.2. Results for the IDTM

In the first set of simulations (runs 1–4), the effect of varying the Lundquist number $S_L$ is explored for a given half-distance between the two current layers, that is $l = 0.25$ (leading to $l/a = 5$). Typically, a range of $S_L$ values varying between $1 \times 10^5$ to $3 \times 10^6$ is employed (listed in Table 1). According to our stability analysis, this would correspond to critical aspect ratio values $L/a = S_L^{9/29}$ ranging between 25 and 102. However, in order to avoid the extra diffusion of the

equilibrium during the simulations at the lowest Lundquist numbers, we take aspect ratios that are significantly higher than the critical values (by typically a factor of two). We believe that this choice does not strongly influence our results as long as these aspect ratios remain much smaller compared to the corresponding SP ones, that is, by one order of magnitude for our different runs (see Table 1).

First, we measured the growth rate of the perturbations as a function of time. This is illustrated in Figure 5 for run 4, where we plot the maximum absolute value of the $B_x$ component, which is assumed to be a reasonable diagnostic for the instability magnitude. The linear phase can be clearly identified. Indeed, between $t \approx 2t_*$ and $t \approx 10t_*$, a growth with $B_x \propto \exp(0.17t/t_*)$ is obtained that compares well with the theoretical linear growth rate value of $\gamma_m t_* = 0.152$ expected from linear theory. As one can see in Table 1, the agreement between measured growth rate values and expected linear growth rate is also rather good for the three other runs (runs 1–3). The numerical values are slightly lower than the theoretical ones for runs 1 and 2 due to the resistive diffusion of the equilibrium for these configurations. Note that our numerically obtained growth rates can also be expressed using the Alfvén time $t_A$, based on the critical length value of each run. This leads to values of $\gamma_m t_A$ in the range [0.63 : 0.87], which slightly overestimates the maximum theoretical value expected for $l/a = 5$ of 0.55 because of our choice of supercritical aspect ratios. Nevertheless, our simulations clearly confirm that the IDTM can grow on a timescale of the order of $t_A$, as expected from linear theory. For the wavenumber of the linearly dominant mode, we have determined the fastest $n$ mode among the spectrum of 10 initially perturbed ones and deduced the corresponding value via $k_n = 2\pi n / L_y$. The results reported in Table 1 show rather good agreement. Moreover, we have investigated the profile of the perturbations during the linear growth phase. This is illustrated in Figure 6 using the $B_x$ magnetic field component, for run 1 at $t = 6t_*$. One can see a contour plot over the numerical domain on the left panel, and a corresponding transverse cut $B_x(x)$ over the left current layer on the right panel. The perturbation shape on each current layer is very similar to the eigenmode obtained from single ideal tearing mode theory (Landi et al. 2015), with an additional anti-symmetric deformation due to the presence of the second layer (the right current layer for the present case). The three other runs (runs 2–4) display similar structural features.

As seen in Figure 5, the linear phase is followed by a secondary rapid growth occurring again on a few Alfvénic
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Contour map of the thermal pressure for run 7 at different times of the RDTM development. The magnetic field vector is also overlaid.</figcaption></figure>

times (typically between $t \approx 20t_*$ and $t \approx 35t_*$). At the end of this second phase, the $B_x$ component is able to reach a maximum value of order unity, which is a magnitude comparable to the equilibrium component. A similar behavior is observed for the three other runs, as one can see in the left panel of Figure 7; however, there are some slight differences in the chronology, due to the different linear growth rates. In order to understand in more detail the temporal dynamics, in Figure 8 we have also drawn contour plots taken at different times for run 3 using the thermal pressure. One can clearly see the growing interaction of plasmons situated on the two layers, once they have reached a critical size. This is possible due to the inverse cascade effect (merging between plasmons on each layer) previously mentioned as occurring for a single layer. A cascade effect toward small scales (see the very small plasmons at $t = 32t_*$, for example) is also observed at the same time. The further evolution of the double-current layer system is different, however. Indeed, the interaction between the two layers is finally observed to violently disrupt the magnetic configuration situated in the central region at $t \approx 35t_*$. This strong interaction of the two layers via the plasmons seems similar to the process occurring during the nonlinear evolution of RDTM, where a secondary instability is structurally driven when a critical deformation of the magnetic islands is reached (Janvier et al. 2011). However, these are differences that should be addressed in a future study, as the process involves two layers with a complex irregular structure (i.e., a few plasmons of different sizes for IDTL versus one single island for RDTM). Then, a relaxation phase

(corresponding to $B_x$ decaying with time) occurs in which the inner region progressively retrieves a purely longitudinal magnetic field with the same direction as the outer regions. The final state (last panel of Figure 8) is finally reached when the central region is free of field-line reversal.

In a second set of simulations (runs 5 and 6), we have investigated the effect of varying the relative distance $l/a$ and compared it to run 4. The results using the maximum $B_x$ are reported in the right panel of Figure 7. First, the linear phase is not strongly dependent on the $l/a$ value and again our measured growth rate values follow the expected theoretical ones (see Table 1). However, this is different for the disruption time taken at the maximum obtained $B_x$. Indeed, the two layers interact sooner and later when $l/a$ is decreased and increased, respectively. This result is not surprising, as the plasmons growing on the two layers can interact once they are closer. The obtained disruption times are $t \approx 25t_*$ and $t \approx 45t_*$ for $l/a = 2.5$ and $l/a = 10$ cases, respectively. The final relaxed states are similar to the previously described runs.

### 3.3. A Comparison between IDTM and RDTM

In order to compare the results obtained for the IDTM in the previous section with RDTM counterparts developing in much smaller aspect ratio configurations, two additional runs (runs 7 and 8) are then carried out for two Lundquist number values. The numerical procedure is exactly the same as previously used. As the relative closeness between the tow layers is $l/a = 2.5$, the RDTM simulations can be closely compared to the previous (IDTM) run 5. An aspect ratio of $L/a = 20$ is also
taken as $L_y = 1$, leading to a linearly dominant wavenumber $k_m a = 0.31$ corresponding to a $n = 1$ mode, which is situated midway between the linearly fastest wavenumber $k_m a$ and the cutoff wavenumber $ka \approx 0.7$. One must also note that the value of $2L/l = 8$ exceeds the critical value of 5–6 necessary to trigger an explosive nonlinear phase (Janvier et al. 2011).

First, the time evolutions of the maximum $B_x$, as well as the maximum velocity $V_y$ component, are plotted in Figure 9 for runs 7 and 8 and compared to run 5. The results show that perturbations initially grow on a slow timescale, in agreement with the expected theoretical linear growth rates (see values in Table 1). Indeed, one must remember that classical RDTMs are known to follow an initial phase (including linear and early nonlinear growth) with an SP-like rate scaling as $\eta^{1/2}$ (Wang et al. 2007). A second, even slower phase is subsequently observed that was been previously identified as a Rutherford regime. During this Rutherford phase, the velocity perturbation in the right panel is slightly decaying while the perturbed magnetic field continues to grow on the extremely slow, purely diffusive timescale (i.e., a rate scaling as $\eta^1$). Contour maps of the thermal pressure (see Figure 10) clearly show the corresponding magnetic islands slowly growing on the two layers (see panels at $t = 50$ and $t = 80$) for run 7. Note that the perturbation is a $n = 1$ mode, as expected from linear theory, because $n > 1$ modes are stable for such large aspect ratio systems.

A third phase characterized by rapid growth starts at $t \approx 80t_*$ and $t \approx 100t_*$ for runs 7 and 8, respectively. Our results for this explosive evolution agree with all the characteristic features reported in detail elsewhere (Wang et al. 2007; Janvier et al. 2011; Zhang & Ma 2011). For example, the timescale of the phase is typically a few Alfvén times and is the same for both runs, thus independent of the resistivity. The maximum perturbation amplitude (in $V_y$ or $B_x$) that is reached at the end of the phase is also the same, independent of the run considered. The triggering is delayed for run 8 with respect to run 7, as the previous phases (linear, early nonlinear, and Rutherford) are strongly resistivity-dependent. The associated mechanism is due to a secondary instability that is structurally driven once a critical deformation of the two magnetic islands is reached (Janvier et al. 2011). As one can see in Figure 10, a magnetic reconnection is therefore observed to complete the process (at roughly the maximum peak perturbation), with an elongated current layer merging with the opposite island (panels at $t = 86t_*$ and $t = 88t_*$ in Figure 10). When compared to the dynamics of the IDTM described in the previous subsection, we thus conclude that the explosivity and fastness are comparable.

A final relaxed state (panel at $t = 100t_*$ in Figure 10) is finally reached, and is very similar to the previously relaxed states resulting from the IDTM disruptions described in the previous subsection. However, note that the relaxation time-scale is at least two times faster for the RDTM compared to the IDTM. This can be explained by the fact that the two layers are merging through bigger islands for RDTM cases compared to IDTM cases, thus leading to a more efficient interaction.

## 4. Conclusion

In this paper, we have studied the development of MHD tearing instabilities that occur in double-current sheet systems. We focus on configurations with relative current sheet

separation $l$ ranging between approximately 1 and 10 times the half-width of each current layer $a$.

First, IDTM is shown to linearly grow on an Alfvénic timescale in the asymptotic high Lundquist number limit ($S_L \to \infty$). This is already the case for the typical values used in this study of $S_L \approx 10^6$, as long as unstable configurations with aspect ratios $L/a$ scaling as $S_L^{9/29}$ are taken. Numerical simulations confirm the linear phase development expected from our linear analysis, extending previous studies that found that the single ideal tearing mode has a critical aspect ratio scaling as $S_L^{1/3}$ (Pucci & Velli 2014; Landi et al. 2015; Del Sarto et al. 2016). More explicitly, the linear growth rate value $\gamma_m t_A = 0.5 - 0.8$ varies slightly with $l/a$ and remains close to the single ideal tearing value of 0.62. The corresponding wavenumber is shown to be $k_m L \approx \beta S_L^{6/29}$ (with $\beta = 0.35 - 0.77$), and differs significantly from the $k_m L \approx 1.4 \times S_L^{1/6}$ obtained for the single-mode counterpart.

The nonlinear evolution of IDTM observed in our numerical simulations reveals some shared features with results obtained for single-current systems. Indeed, each current layer experiences a double cascade effect, as primary plasmons merge to form larger plasmons, while, simultaneously new secondary plasmons produce a cascading accelerating reconnection process toward small scales. The final issue is different, however, as a disruption violently terminates the evolution when the two current sheets interact, leading to a relaxation toward a state free of field-line reversals. In a single-current sheet layer, such numerical simulation occurs when a single plasmod reaches the largest size of the periodic domain (Landi et al. 2015). To our knowledge, the process of merging between the two layers developing plasmod chains has not yet been reported. We believe that it is the natural outcome of the inverse (large-scale) cascade effect on each current layer that drives the system in a structurally unstable state, in a way similar to the secondary instability process at work for RDTM (Janvier et al. 2011). A future study would be necessary to investigate the latter point at length.

Moreover, an interesting point concerns the comparison of IDTM dynamics (typical of long and thin systems) with the RDTM dynamics occurring in similar configurations with smaller aspect ratios (of order $2\pi$). Indeed, a fast explosive reconnection phase is triggered with a fast timescale (independent of the resistivity), which is close to the IDTM Alfvénic timescale. The disruption that characterizes the subsequent interaction of the two current layers is also very similar. And the same conclusion holds for the final state.

Our results are relevant in the context of plasma environments where multiple current layers can form. This is the case for many astrophysical plasma systems, such as the solar photosphere, solar corona, solar wind, and pulsar wind nebula. For example, in the latter, the magnetic field structure around pulsars is believed to display stripes of alternating polarity separated by current sheets (Coroniti 1990). Recently, it has been shown that the variability of the Crab Nebula/pulsar system, seen as flares, can be explained by the explosive phase of the relativistic counterpart of the RDTM (Baty et al. 2013). The dynamics of the relativistic regime is similar to the Newtonian case when the light speed replaces the Alfvèn speed. The difficulty in the latter flare model is the somewhat long initial phase (resistivity-dependent) necessary for the development of magnetic islands before they can enter into the explosive growth phase, as described in this study. Indeed, this
would imply that RDTM instabilities never have enough time to grow, as they are developing in the wind frame that is itself radially propagating outward at a highly relativistic speed. Thus, some fast driven mechanism of reconnection is necessary to accelerate the initial phase, as suggested from propagating shock scenarios (Sironi & Spitkovsky 2014). We suggest that there is an alternate way based on the IDTM, as a similar explosive reconnection can spontaneously operate when the aspect ratio is large enough. In the pulsar system, the striped wind region is considered to contain infinitely long current sheets, as it is an extremely extended region where the magnetic field lines can be wound by a large number of turns (typically more than a few thousand times). Therefore, the triggering for the onset of IDTM dynamics can be easily reached even for a very Lundquist number. In this case, the optimal region for the plasmons' development would probably be situated further from the pulsar compared to the RDTM scenario. In order to account for the observations in a more complete way, future studies are required to investigate particle acceleration and the resulting radiation. Indeed, such explosive reconnection is a promising candidate for leading to significant thermal synchrotron radiation, as well as a nonthermal radiation (Takamoto et al. 2015; Lyutikov et al. 2016).

In the present work, we have followed the classical picture of magnetic reconnection, focusing on 2D dynamics representative of planar current sheets. Consequently, three-dimensional (3D) effects have not been considered. In most real systems, a guiding magnetic field perpendicular to the reconnection plane exists, so that the 3D field does not vanish in the sheets. Previous studies on 3D tearing modes have shown that new instabilities develop, as magnetic islands and plasmons appear as tube-like structures (called flux ropes) that are subject to secondary instabilities (Baalrud et al. 2012). The addition of a guiding field to our double-current sheet configuration could influence the results obtained in the present study, as the cascade of modes is expected to be different. An explosive phase driven by 3D MHD processes associated with a strongly turbulent reconnection could be also at work (Oishi et al. 2015). In a future study, it would thus be of interest to

examine how our present results and conclusions are affected in three dimensions.

We sincerely thank the anonymous referee for useful suggestions that helped improve the content of the paper. H.B. acknowledges support by the French National Research Agency (ANR) through Grant ANR-13-JS05-0003-01 (Project EMPERE). We also acknowledge computational facilities available at Equip@Meso of the Université de Strasbourg.

## References

* Baalrud, S. D., Bhattacharjee, A., & Huang, Y.-M. 2012, PhPl, 19, 022101

* Baty, H., Petri, J., & Zenitani, S. 2013, MNRAS, 436, L20

* Birk, G. T., Konz, C., & Otto, A. 1997, PhPl, 4, 4173

* Coroniti, F. V. 1990, ApJ, 349, 538

* Del Sarto, D., Pucci, F., Tenerani, A., & Velli, M. 2016, JGR, 121, 1857

* Del Zanna, L., landi, S., Papini, E., Pucci, F., & Velli, M. 2016, MNRAS, 460, 3753

* Huang, Y., & Bhattacharjee, A. 2010, PhPl, 17, 062104

* Janvier, M., Kishimoto, Y., & Li, J. Q. 2011, PhRvL, 107, 195001

* Keppens, R., Porth, O., Galsgaard, K., et al. 2013, PhPl, 20, 092109

* Landi, S., del Zanna, L., Papini, E., Pucci, F., & Velli, M. 2015, ApJ, 806, 131

* Loureiro, N. F., Samtaney, R., Schekochihin, A. A., & Uzdensky, D. A. 2012, PhPl, 19, 042303

* Lyutikov, M., Sironi, L., Komissarov, S., & Porth, O. 2016, arXiv:1603.05731

* Nemati, M. J. A., Wang, Z. X., & Wei, L. 2016, ApJ, 821, 128

* Oishi, J. S., mac Low, M.-M., Collins, D. C., & Tamura, M. 2015, ApJ, 806, L12

* Otto, A., & Birk, G. T. 1992, PhFLB, 4, 3811

* Porth, O., Xia, C., Hendrix, T., Moschou, S. P., & Keppens, R. 2014, ApJS, 214, 4

* Priest, E. R., & Forbes, T. G. 2000, Magnetic Reconnection (Cambridge: Cambridge Univ. Press)

* Pucci, F., & Velli, M. 2014, ApJL, 780, L19

* Samtaney, R., Loureiro, N. F., Uzdensky, D. A., Schekochihin, A. A., & Cowley, S. C. 2009, PhRvL, 103, 105004

* Shen, C., & Liu, Z. X. 1998, PhPl, 5, 2466

* Sironi, L., & Spitkovsky, A. 2014, ApJL, 783, L21

* Takamoto, M., Petri, J., & Baty, H. 2015, MNRAS, 454, 2972

* Tenerani, A., Velli, M., Rappazo, A. F., & Pucci, F. 2015, ApJ, 813, L32

* Wang, Z. X., Wang, X. G., Dong, J. Q., et al. 2007, PhRvL, 99, 185004

* Zhang, C. L., & Ma, Z. W. 2011, PhPl, 18, 052303