Fast Magnetic Reconnection: Secondary Tearing Instability and Role of the Hall Term

TL;DR

This paper investigates the ideal tearing instability in Hall-MHD simulations, revealing that the Hall term significantly accelerates magnetic reconnection and energy release in thin current sheets, explaining explosive astrophysical phenomena.

Contribution

It introduces detailed simulations of the Hall-MHD tearing instability, highlighting the role of the Hall term in enhancing reconnection rates and explosive energy release.

Findings

Hall term increases reconnection rate up to ten times compared to linear phase.

Secondary current sheets naturally evolve to ideal aspect ratios at high Lundquist numbers.

Reconnection occurs on super-Alfvénic timescales, explaining explosive astrophysical events.

Abstract

Magnetic reconnection provides the primary source for explosive energy release, plasma heating and particle acceleration in many astrophysical environments. The last years witnessed a revival of interest in the MHD tearing instability as a driver for efficient reconnection. It has been established that, provided the current sheet aspect ratio becomes small enough ($a/L \sim S^{-1/3}$ for a given Lundquist number $S\gg 1$), reconnection occurs on ideal Alfv\'en timescales and becomes independent on $S$. Here we investigate, by means of two-dimensional simulations, the \emph{ideal} tearing instability in the Hall-MHD regime, which is appropriate when the width of the resistive layer $\delta$ becomes comparable to the ion inertial length $d_i$. Moreover, we study in detail the spontaneous development and reconnection of secondary current sheets, which for high $S$ naturally adjust to the ideal aspect ratio and hence their evolution proceeds very rapidly. For moderate low $S$, the aspect ratio tends to the Sweet-Parker scaling ($a/L \sim S^{-1/2}$), in order to fulfill the condition $\delta \ll a$ necessary for the onset of a tearing instability. When the Hall term is included, the reconnection rate of this secondary nonlinear phase is enhanced and, depending on the ratio $d_i/\delta$, can be twice with respect to the pure MHD case, and up to ten times larger than the linear phase. Therefore, the evolution of the tearing instability in thin current sheets in the Hall-MHD regime naturally leads to an explosive disruption of the reconnecting site and to energy release on super-Alfv\'enic timescales, as required to explain astrophysical observations.

Tables

Table 1. Table 1: Physical and numerical parameters of all the simulations performed in this work for linear (L) and nonlinear (N) runs. From left to right: Lundquist number S 𝑆 S , aspect ratio a / L 𝑎 𝐿 a/L , Hall coefficient η H = d i / L subscript 𝜂 𝐻 subscript 𝑑 𝑖 𝐿 \eta_{H}=d_{i}/L , ratio of the ion inertial length d i subscript 𝑑 𝑖 d_{i} with respect to the inner resistive layer δ 𝛿 \delta , domain size across ( L x subscript 𝐿 𝑥 L_{x} ) and along ( L y subscript 𝐿 𝑦 L_{y} ) the current sheet, and grid size ( N x × N y subscript 𝑁 𝑥 subscript 𝑁 𝑦 N_{x}\times N_{y} ). Run 0L starts from a pressure equilibrium (PE) configuration, all the other runs use a force-free (FF) equilibrium.

Run	$S\left({\eta }^{-1}\right)$	$a/L$	${\eta }_H$	$d_i/\delta$	$L_x$	$L_y$	$N_x\times N_y$
0L	${10}^5$	$\sim 0.022$	0.01	$\sim 3.2$	$\sim 0.43$	$\sim 2.2\pi$	$1024\times 128$
1L	${10}^5$	$\sim 0.022$	0	0	$\sim 0.43$	$\sim 2.2\pi$	$1024\times 128$
2L	${10}^5$	$\sim 0.022$	0.002	$\sim 0.6$	$\sim 0.43$	$\sim 2.2\pi$	$1024\times 128$
3L	${10}^5$	$\sim 0.022$	0.005	$\sim 1.6$	$\sim 0.43$	$\sim 2.2\pi$	$1024\times 128$
4L	${10}^5$	$\sim 0.022$	0.01	$\sim 3.2$	$\sim 0.43$	$\sim 2.2\pi$	$1024\times 128$
5N	${10}^5$	$\sim 0.022$	0	0	$\sim 0.43$	$\sim 2.2\pi$	$4096\times 512$
6N	$6.7\cdot {10}^5$	$\sim 0.011$	0	0	$\sim 0.23$	$\sim 1.14\pi$	$4096\times 512$
7N	$8\cdot {10}^5$	$\sim 0.011$	0	0	$\sim 0.22$	$\sim 1.08\pi$	$4096\times 512$
8N	${10}^6$	0.01	0	0	0.2	$\pi$	$4096\times 512$
9N	${10}^6$	0.01	0.0002	$\sim 0.2$	$0.2$	$\pi$	$4096\times 512$
10N	${10}^6$	0.01	0.0014	$\sim 1.4$	$0.2$	$\pi$	$4096\times 512$
11N	${10}^6$	0.01	0.003	$\sim 3.0$	$0.2$	$\pi$	$4096\times 512$
12N	$2\cdot {10}^6$	$\sim 0.008$	0	0	$\sim 0.16$	$\sim 0.79\pi$	$4096\times 512$
13N	$5\cdot {10}^6$	$\sim 0.006$	0	0	$\sim 0.12$	$\sim 0.58\pi$	$4096\times 512$
14N	${10}^7$	$\sim 0.0046$	0	0	$\sim 0.093$	$\sim 0.46\pi$	$4096\times 512$
	${10}^8$	$\sim 0.0022$	0	0	$\sim 0.043$	$\sim 0.22\pi$	$4096\times 512$

Equations

$${\tau }_{SP}={\tau }_A(c_A/v_{\mathrm{in}})={\tau }_A(L/a)={\tau }_A{S}^{1/2}=\sqrt{{\tau }_A{\tau }_D},$$

$${\tau }_t\sim {{\tau }_A}_a{S}_a^{1/2}=(a/L{)}^{3/2}{\tau }_A{S}^{1/2},$$

$$v_{\mathrm{in}}/c_A=a/L\le S^{-1/3}$$

$$\gamma {\tau }_A\simeq 0.63S^{(3\alpha -1)/2},$$

$$\gamma {\tau }_A\simeq 0.63,$$

$${\mbox{\boldmath${\mathrm{v}}$}}_{e}={\mbox{\boldmath${\mathrm{v}}$}}-\frac{{\mbox{\boldmath${\mathrm{J}}$}}}{en_{e}},\qquad{\mbox{\boldmath${\mathrm{J}}$}}=\frac{c}{4\pi}\nabla\times{\mbox{\boldmath${\mathrm{B}}$}}$$

$$\partial_{t}{\mbox{\boldmath${\mathrm{B}}$}}=\nabla\times\left({\mbox{\boldmath${\mathrm{v}}$}}_{e}\times{\mbox{\boldmath${\mathrm{B}}$}}\right)+\eta\nabla^{2}{\mbox{\boldmath${\mathrm{B}}$}}.$$

$$\displaystyle\partial_{t}\rho+\nabla\cdot{(\rho{\mbox{\boldmath${\mathrm{v}}$}})}=0$$

$$\displaystyle\rho\left(\partial_{t}+{\mbox{\boldmath${\mathrm{v}}$}}\cdot\nabla\right){\mbox{\boldmath${\mathrm{v}}$}}=-\nabla P+(\nabla\times{{\mbox{\boldmath${\mathrm{B}}$}}})\times{\mbox{\boldmath${\mathrm{B}}$}}$$

$$\displaystyle\left(\partial_{t}+{\mbox{\boldmath${\mathrm{v}}$}}\cdot\nabla\right)T=(\Gamma\!-\!1)\!\left[-(\nabla\cdot{{\mbox{\boldmath${\mathrm{v}}$}}})T\!+\!\frac{1}{S}\frac{|\nabla\times{{\mbox{\boldmath${\mathrm{B}}$}}}|^{2}}{\rho}\right]$$

$$\displaystyle\partial_{t}{\mbox{\boldmath${\mathrm{B}}$}}=\!\nabla\!\times\!\left({\mbox{\boldmath${\mathrm{v}}$}}\times{\mbox{\boldmath${\mathrm{B}}$}}\right)\!+\!\frac{1}{S}\nabla^{2}{\mbox{\boldmath${\mathrm{B}}$}}-\eta_{H}\nabla\!\times\!\frac{(\nabla\!\times\!{\mbox{\boldmath${\mathrm{B}}$}})\!\times\!{\mbox{\boldmath${\mathrm{B}}$}}}{\rho},$$

$$d_i=\frac{c}{{\omega }_{pi}}=c\sqrt{\frac{m_i}{4\pi n_i{e}^2}},$$

$$\delta /a\simeq S_a^{-3/10}{\Delta }^{\prime 1/5},$$

$${\Delta }^{\prime }=2\left(1/ka-ka\right).$$

$$d_i/\delta ={\eta }_H(L/a)S_a^{1/4}={\eta }_H{S}^{(3\alpha +1)/4}$$

$$d_i/\delta ={\eta }_H{S}^{1/2}$$

$$\displaystyle{\mbox{\boldmath${\mathrm{B}}$}}_{0}=\tanh(x/a){\mbox{\boldmath${\mathrm{\hat{y}}}$}}+\mathrm{sech}(x/a){\mbox{\boldmath${\mathrm{\hat{z}}}$}},$$

$$\varphi =\epsilon \mathrm{sech}(x/a)\sum _{n=1}^N\cos (k_{n}y+{\phi }_n),$$

$$\gamma (v_x^{\prime \prime }-k^2{v}_x)=$$

$$\gamma b_x=$$

$$+\frac{{\eta }_H}{{\rho }_0}\left(k^2{B}_{0y}{b}_z-ik{B}_{0z}^{\prime }{b}_x\right),$$

$$\gamma b_z=$$

$$+\frac{{\eta }_H}{{\rho }_0}\left[B_{0y}(b_x^{\prime \prime }-k^2{b}_x)-B_{0y}^{\prime \prime }{b}_x\right],$$

$$\gamma v_z=$$

$$1=\mathrm{P}(b_x^{\mathrm{R}})=-\mathrm{P}(b_z^{\mathrm{R}})=-\mathrm{P}(v_x^{\mathrm{I}})=\mathrm{P}(v_z^{\mathrm{I}}),$$

$$1=-\mathrm{P}(b_x^{\mathrm{I}})=\mathrm{P}(b_z^{\mathrm{I}})=\mathrm{P}(v_x^{\mathrm{R}})=-\mathrm{P}(v_z^{\mathrm{R}})$$

$${\gamma }_{\mathrm{rec}}=\frac{1}{N}\sum _{i=1}^N\frac{1}{{\Phi }_i}\frac{\mathrm{d}{\Phi }_i}{\mathrm{d}t},$$

$$\displaystyle{\mbox{\boldmath${\mathrm{B}}$}}^{*}(x)=\mathrm{B}_{0}^{*}\tanh(x/a^{*}){\mbox{\boldmath${\mathrm{\hat{y}}}$}},$$

Full text

plasmas — magnetohydrodynanics (MHD) — magnetic reconnection — instabilities

\turnoffedits

1 Introduction

The rapid conversion of magnetic energy into heat and particle acceleration is often

On macroscopic scales, magnetized astrophysical plasmas are invariably modeled by using the MHD approximation, with a finite conductivity to be employed in Ohm’s law. However, in astrophysical systems the magnetic diffusivity $\eta$ is so small that the diffusion time $\tau_{D}=L^{2}/\eta$ is incomparably longer than the (ideal) dynamical time scale $\tau_{A}=L/c_{A}$ required to explain such phenomena (here $L$ is the macroscopic length scale and $c_{A}$ the Alfvén speed). The presence of localized strong current sheets can speed up the magnetic annihilation by the mechanism of reconnection. Unfortunately, the classical MHD mechanisms for reconnection, namely the non-linear steady-state model by Sweet and Parker (Parker, 1957; Sweet, 1958), SP from now on, and the linear tearing instability (Furth et al., 1963) both predict a very inefficient reconnection rate, and the search for efficient reconnection has steadily moved from macroscopic MHD to kinetic regimes (e.g. Yamada et al., 2010).

The steady-state SP model for incompressible magnetic reconnection driven by a constant velocity inflow $v_{in}$ in a current sheet of length $L$ and width $a$, predicts a reconnection time which increases with the Lundquist number $S=\tau_{D}/\tau_{A}\gg 1$ as

[TABLE]

far too slow to explain the observed flare-like events, given that for astrophysical plasmas the usual estimate is $S\sim 10^{12}$ (note that the SP model also implies that the scaling for the aspect ratio of the current sheet is $L/a=S^{1/2}$). On the other hand, current sheets are known to be locally prone to the linear tearing instability, which leads to the formation of X-points and magnetic islands (also called plasmoids) during the reconnection process. The $e$-folding time $\tau_{t}$ of the fastest growing mode, calculated by using the current sheet half thickness $a$ as characteristic length, reads

[TABLE]

where ${\tau_{A}}_{a}=a/{c_{A}}$ and $S_{a}={ac_{A}}/{\eta}$. As far as $a$ is of the order of the macroscopic scale and $S\gg 1$, the timescale is again too large .

, it has been recognized that, for high Lunquist numbers $S~{}\gtrsim~{}10^{4}$, a SP like current sheet undergoes tearing instability which, once measured on the relevant scale $L$, is very fast. Indeed, Eq. (2) predicts a super-Alfvénic linear growth rate for the tearing instability $\gamma\sim 1/\tau_{t}>1/\tau_{A}$, and even increasing with $S$ as $\gamma\tau_{A}\sim S^{1/4}$. This implies a very efficient reconnection and an explosive nonlinear production of a chain of fast moving, merging plasmoids (the so-called plasmoid instability). This result is of course paradoxical, since the ideal MHD case (where reconnection is forbidden) cannot be retrieved for $S\to\infty$, and the only possibility to resolve this puzzle is that the SP current sheet cannot form in any dynamical thinning process (Pucci & Velli, 2014; Tenerani et al., 2016; Uzdensky & Loureiro, 2016; Landi et al., 2017).

Indeed, by analyzing the characteristic timescales involved in the dynamic evolution of a whenever the condition

[TABLE]

holds, the tearing instability evolves on super-ideal timescales and the sheet is disrupted. Therefore the SP configuration, even thinner than this critical threshold, can never form. We note that a criterion similar to Eq. (3) was found to hold in the formation of secondary tearing instabilities in a plasmoid-induced reconnection model (Shibata & Tanuma, 2001). If one assumes that the inverse aspect ratio scales with $S$ as $a/L=S^{-\alpha}$, the growth rate of the most unstable mode is

[TABLE]

where the numerical factor arises from the detailed analysis of the tearing instability . Notice that the SP case is correctly retrieved for $\alpha=1/2$. The linear phase of the tearing instability for the critical case $a/L=S^{-1/3}$, named ideal tearing instability, was first examined analytically by Pucci & Velli (2014), who calculated the eigenmodes and found that the growth rate of the fastest reconnecting mode indeed tends asymptotically (for $S\to\infty$) to

[TABLE]

that means reconnection on the macroscopic Alfvénic timescales. This result has been also retrieved and extended to the nonlinear regime using numerical simulations (Landi et al., 2015; Del Zanna et al., 2016a; Landi et al., 2017).

Related works have analyzed the evolution during the collapse of a current sheet (Tenerani et al., 2015b), the dependence on viscosity (Tenerani et al., 2015a) and on the equilibrium profile (Pucci et al., 2018), the inclusion of electron inertia (Del Sarto et al., 2016; Del Sarto & Ottaviani, 2017) and the extension to the relativistic regime, in which the linear and nonlinear cases have been analyzed for the first time (Del Zanna et al., 2016b).

Since the critical inverse aspect ratio can be very small, it is important to study how the ideal tearing instability is affected when approaches the ion inertial length $d_{i}$. In general, the growth rates of the tearing mode instability are known to be larger in the appropriate Hall regime (Terasawa, 1983; Shay et al., 2001; Shaikhislamov, 2004). A linear analysis for the thin current sheets required for fast reconnection has been performed in the Hall-MHD regime by Pucci et al. (2017). The scaling for the growth rates now depends on $d_{i}$ as well and the growth is confirmed to be faster with respect to the MHD case. Preliminary nonlinear simulations can be found in Papini et al. (2018), where it is shown that secondary instabilities are also more rapidly evolving when the Hall effect is included.

In the present paper we extend these works and study in detail, through two-dimensional MHD and Hall-MHD simulations, the development and the nonlinear stage of the tearing instability in critical current sheets with $a/L\sim S^{-1/3}$. In particular, we concentrate on the physical conditions holding at the time of the onset of secondary tearing instabilities occurring inside the main reconnecting sheet.

2 Equations and numerical setup

2.1 Hall-MHD model of the tearing instability

While the macroscopic MHD approximation is a one-fluid model, when spatial ion scales are reached the electron and the ion velocities decouple. When that happens, it is the electron velocity, defined by

[TABLE]

($n_{e}$ is the numerical density of electrons, $e$ is the unsigned fundamental electrical charge, and $c$ is the speed of light), that drives magnetic evolution by entering the induction equation

[TABLE]

The full system of compressible, nonlinear Hall-MHD equations then becomes

[TABLE]

where $\Gamma$ is the adiabatic index and the other variables retain their usual meaning. All quantities have been here normalized against the Alfvénic ones $L$, $B_{0}$, $\rho_{0}=m_{i}n_{0}$, $c_{A}=B_{0}/\sqrt{4\pi\rho_{0}}$, $P_{0}=\rho_{0}c_{A}^{2}$, $T_{0}=(k_{B}/m_{i})P_{0}/\rho_{0}$, with $m_{i}$ being the mass of the ions constituting the plasma. The Hall coefficient is defined as $\eta_{H}=d_{i}/L$, where $d_{i}$ is the reference value of the ion inertial length

[TABLE]

which depends on the plasma frequency $\omega_{pi}$ of ions and in turn on $n_{i}\equiv n_{e}$ ($=n_{0}$ for $d_{i}$).

As anticipated, the Hall term is not negligible when the ion inertial length $d_{i}$ becomes comparable to the width $\delta$ of the inner resistive layer of the tearing instability, which is smaller than the sheet’s half thickness $a$. For the fastest growing mode, the inner width $\delta$ is described by the equation (Biskamp, 1993)

[TABLE]

where $S_{a}$ is the Lundquist number employed in Eq. (2) and $\Delta^{\prime}$ is an instability parameter which depends on the configuration considered for the equilibrium magnetic field. $\Delta^{\prime}$ may depend on $ka$, which scales as $ka\sim S_{a}^{1/4}$ for the fastest growing mode. For the Harris sheet configuration commonly employed in numerical works, including the present one (see the initial conditions further on)

[TABLE]

Thus, at high Lundquist numbers $S_{a}$ and for the fastest growing mode we find $\Delta^{\prime}\sim S_{a}^{1/4}$ and $\delta/a\sim S_{a}^{-1/4}$, hence the ratio

[TABLE]

is the quantity that determines whether Hall effects are important in the dynamics of reconnection (Terasawa, 1983). Here the second expression is referred to the generic aspect ratio $a/L=S^{-\alpha}$ considered above, recalling that $S_{a}=(a/L)S$. Notice that in the critical case $\alpha=1/3$ we find

[TABLE]

and we can identify three distinct regimes: an MHD regime ($\eta_{H}\ll S^{-1/2}$), where the Hall term does not play a relevant role, a mild Hall regime ($\eta_{H}\lesssim S^{-1/2}$), where the ion inertial length is comparable to the thickness of the inner layer, and a strong Hall regime ($\eta_{H}>S^{-1/2}$), where reconnection is dominated by the Hall effect and the classic theory of the tearing instability is no longer valid (see also Shaikhislamov, 2004). Recently, Pucci et al. (2017) extended the ideal tearing instability to include the Hall term in the case of a Harris current sheet in pressure equilibrium, that is an unperturbed magnetic field which is unidirectional. They found the existence of the regimes discussed above, and showed that the linear growth rate starts to increase roughly for values $d_{i}/\delta\sim 3$ . However, this threshold is likely to be actually even smaller, since Hall currents may affect the subsequent nonlinear evolution, where thinner current sheets that formed between the ejected plasmoids may eventually host secondary reconnection events, as we will show later in Section 5.

2.2 Numerical setup

In the present work we consider nonlinear simulations with initial condition for the magnetic field given by a two-dimensional force-free (FF) current sheet configuration, centered at $x=0$ and asymptotically aligned in the $y$-direction with |{\mbox{\boldmath{\mathrm{B}}}}|=B_{0}=1 according to the profile

[TABLE]

which reproduces a Harris profile for the in-plane component, but keeps a constant magnetic (and thus total) pressure by rotating the magnetic field around the $x$-axis. Moreover, we assume homogeneous density $\rho=1$, pressure and temperature $P=T=\beta/2$, and we do not impose initial velocities ({\mbox{\boldmath{\mathrm{v}}}}_{0}=0). The plasma beta is chosen as $\beta=0.5$, the adiabatic index is $\Gamma=5/3$, and we investigate the case appropriate for the ideal tearing, thus we choose the half thickness of the current sheet $a=S^{-1/3}=\eta^{1/3}$, for $L=1$. We follow the evolution of the plasma by integrating the system of Hall-MHD equations (8 -11) in a $[-L_{x},L_{x}]\times[0,L_{y}]$ domain in the $xy$-plane, but retaining all components of the 3D vectors. The in-plane magnetic field is evolved through a scalar potential $\phi$ (the $z$ component of the vector potential), so that $B_{x}=\partial_{y}\phi$ and $B_{y}=-\partial_{x}\phi$, in order to preserve the solenoidal constraint for the magnetic field.

At the beginning of the simulation, the tearing instability is triggered by in-plane magnetic perturbations localized inside the current sheet. In terms of the scalar potential these perturbations take the form

[TABLE]

where $k_{n}=2\pi n/L_{y}$ and $\varphi_{n}$ is a random phase different for each value of $n$. We choose $N=10$ and $\varepsilon=10^{-6}$ (the overall amplitude of the perturbed magnetic field is $\sim 10^{-4}$). The value of $L_{y}$ is chosen such that the lowest wavenumber resolved for the tearing instability is $ka=2\pi a/L_{y}=0.02$. This value is more than sufficient to capture the fastest growing mode of the tearing instability for the values of $S$ considered in this work. In the $x$-direction we set $L_{x}=20a$ to have boundaries sufficiently far from the reconnection region while retaining the high resolution required inside the current sheet.

The Hall-MHD equations are numerically solved by means of the same MHD code we used in Landi et al. (2015), modified to include the Hall term. Spatial derivatives are calculated using Fourier decomposition along the periodic direction and a fourth-order compact scheme (Lele, 1992) across the current sheet. Time integration is performed with a third-order Runge-Kutta scheme, taking into consideration the effect of the Hall term in the definition of the timestep. Boundary conditions are periodic along $y$ (only integer numbers of wavelengths are thus allowed) and of free outflow at $x=\pm L_{x}$, using the method of projected characteristics (Poinsot & Lele, 1992; Del Zanna et al., 2001; Landi et al., 2005). Unless differently specified, the employed grid consists of $N_{x}\times N_{y}=4096\times 512$ points, which allows us to resolve secondary reconnection events in both the $x$- and the $y$-directions. Table 1 report the full set of simulations used in this work.

3 Linear phase

We now describe the evolution of the linear tearing instability in the case of the ideal limit $a/L=S^{-1/3}$. The results of this section confirm the findings of Pucci et al. (2017), where pressure equilibrium was imposed, and extend them to the force-free case, more appropriate for magnetically dominated systems, employed here and in many other works. The initial magnetic equilibrium is here considered in the general form {\mbox{\boldmath{\mathrm{B}}}}_{0}=(0,B_{0y}(x),B_{0z}(x)). The governing equations of the linear tearing instability read

[TABLE]

and for a given set of parameters $k,S,$ and $\eta_{H}$, the above system of equations constitute a 12th-order two-points eigenvalue problem. Here $v_{x},v_{z},b_{x},$ and $b_{z}$ are the (complex) eigenfunctions of the $x$ and the $z$ components of velocity and magnetic field respectively, and the apex denotes differentiation with respect to $x$. Each eigenmode perturbation, e.g. $b_{x}$, has the form $b_{x}(x,y,t)=b_{x}(x)e^{\gamma t+iky}$, where $k$ is the wavenumber associated to the perturbation in the $y$-direction and the eigenvalue $\gamma$ is the corresponding linear growth rate. In the following we will assume the same settings employed in our numerical simulations, that is, a FF equilibrium (see Eq. (17)).

The above set of equations holds also for the case of pressure equilibrium (PE), in which all terms with $B_{0z}^{\prime}$ vanish so that the equations simplify to the analogous ones employed in Pucci et al. (2017): the eigenvalue problem reduces to 6th-order and, in analogy with the MHD case, the eigenfunctions for the magnetic field and for the velocity are purely real and imaginary, respectively. Moreover, as for the MHD classical tearing instability, we see that in a PE configuration the presence of a constant guide field is ineffective. In the FF case, however, three additional terms involving $B_{0z}^{\prime}$ appear in the equations, and these will lead to different results, especially on the parity of some of the eigenfunctions.

In all simulations, a linear tearing instability develops at the beginning and with the same qualitative behavior. There are, however, some quantitative differences due to the presence of the Hall term. To highlight them, we calculated the linear dispersion relation for four simulations with the same value of $S=10^{5}$ but different values of $\eta_{H}=0$ (MHD case), $0.002,0.005,\text{ and }0.01$, corresponding to $d_{i}/\delta=0,0.6,1.6,\text{and }3.2$, respectively (see Run 1L-4L of Table 1). The linear growth rate $\gamma$ has been calculated by taking, at each time of the linear phase, the modulus of the Fourier transform along the $y$-direction of the average along the $x$-direction of $B_{x}$, since its eigenfunctions are even with respect to $x$ and since $B_{x}$ has no equilibrium component. An exponential fit has been then performed, separately for each Fourier component, to obtain the linear growth rate $\gamma$. These dispersion relations are reported in Fig. 1 and have a similar shape in all cases. In general, for larger $\eta_{H}$ values the corresponding curve yields larger values of $\gamma$. More precisely, the linear growth rate of each mode increases when $d_{i}$ exceeds the thickness of the inner resistive layer $\delta$, up to about $25\%$ more than the MHD case for $ka=0.14$ and $\eta_{H}=0.01$. This is in qualitative and quantitative agreement with Pucci et al. (2017), even though the initial equilibrium here is different and therefore the linear evolution may also be different, due to the additional terms present in Eqs. (19-22). Notice that the results are lower than expected theoretically. For instance, in the MHD case we observe a maximum rate roughly $20\%$ lower than the value predicted by theory. This difference was also encountered in Del Zanna et al. (2016b) and it is due to the diffusion of the initial equilibrium during the evolution. More accurate results were obtained in Landi et al. (2015), where this effect was properly treated.

Differences between the FF and the PE equilibrium arise in the eigenfunctions, shown in Fig. 2. The eigenfunctions have been obtained by using a linearized version of our Hall-MHD code (Landi et al., 2005) and are quantitatively and qualitatively similar to the ones observed in the linear phase of the fully nonlinear simulations. Indeed, the eigenfunctions $b_{x}$, $v_{x}$, $b_{z}$ and $v_{z}$ extracted by a numerical simulation with a PE configuration are even, odd, odd, and even, respectively, as in Pucci et al. (2017). Moreover, $v_{x}$ and $v_{z}$ are purely imaginary, while $b_{x}$ and $b_{z}$ are real. The parity relations can be written as

[TABLE]

where $\mathrm{P}$ denotes the parity operator, whereas the ’$\mathrm{R}$’ and ’$\mathrm{I}$’ superscripts indicate the real and the imaginary part, respectively. In the FF configuration, the eigenfunctions are complex. The parity relations (23) hold also in the FF case, complemented by the relations

[TABLE]

for the imaginary part of $b_{x}$ and $b_{z}$ and for the real part of $v_{x}$ and $v_{z}$.

4 Nonlinear phase: General Properties

We now focus on the nonlinear phase of the instability, and consider simulations with a higher Lundquist number, $S=10^{6}$, so that the settings for the ideal tearing lead to a half thickness $a=S^{-1/3}L=0.01\,L$. We firstly illustrate the general properties by discussing the results of the purely MHD case (Run 8N of Table 1), while differences due to the Hall effects will be discussed in .

In all simulations, as the linear phase evolves, the amplitudes of the perturbations increase exponentially, until the tearing instability saturates and the nonlinear phase begins, as shown in Fig. 3. There, two snapshots of the MHD simulation are taken at the beginning of the nonlinear phase, and a colored contour plot of $J_{z}$ is shown. At time $t=9~{}\tau_{A}$ (bottom left panel) the plasmoids have a size comparable to the thickness of the current sheet, and some of them have already merged. Among the plasmoids we also observe that secondary current sheets have formed, with a thickness of roughly one tenth of the original thickness. One of them is shown in the top left panel of the same figure, by zooming in the region centered at $y\simeq 2.3\,L$. The subsequent evolution is characterized by the coalescence and nonlinear growth of the plasmoids, but the most important feature is the onset of secondary reconnection events in the newly formed current sheets, which then drive the dynamics and eventually lead to the disruption of the whole system. These secondary tearing instabilities are indeed very fast, since already at time $t=9.875~{}\tau_{A}$, in less than one macroscopic Alfvén time, they are fully developed (see the right panels). A more detailed analysis of the evolution of these secondary current sheets is performed in section 5.

Let us now provide a more quantitative support to the above statements. We define the averaged reconnection rate as the quantity

[TABLE]

obtained by taking, at each time, the logarithmic time derivative of the reconnected flux $\Phi_{i}$ between the $i$-th pair of X- and O-points and then averaging over the number of pairs, $N$, in the main current sheet. Here, $\Phi_{i}$ is the difference between the scalar potential $\phi$ at the O- and the X-point in the $i$-th pair (for more details and applications to simulations of plasma turbulence see the Appendix of Papini et al., 2019). The top panel of Figure 4 shows $\gamma_{\mathrm{rec}}$ for our MHD reference run. As we can see, after the initial perturbations have rearranged to select the fastest tearing eigenmodes, $\gamma_{\mathrm{rec}}$ reaches a plateau with a value $\gamma_{\mathrm{rec}}\,\tau_{A}\simeq 0.4$. It is also possible to identify a second, more noisy, plateau between $t\simeq 9\,\tau_{A}$ and $t\simeq 10\,\tau_{A}$, roughly the temporal range selected in Fig. 3, that we interpret as the average reconnection rate of the secondary current sheets.

In order to support this conclusion, we estimated the reconnection rates by performing an exponential fit of the root-mean-square (rms) value of the $x$-component of the magnetic field, that we name $B_{x,\text{rms}}$, which is a good proxy of the reconnection rate. The bottom panel of Fig. 4 shows this quantity as a function of time. The black curve denotes the primary current sheet, while the blue curve has been calculated by restricting to the secondary current sheet (the values are lowered by a factor of 10 for ease of presentation). In the latter case we notice a steepening at $t\simeq 9.5~{}\tau_{A}$, a clear signature of the secondary tearing instability. The horizontal dashed and dot-dashed lines in the top panel, with values $\gamma_{\mathrm{rec}}\,\tau_{A}=0.48$ and $2.30$, respectively, correspond to the exponential fits indicated by the red dashed lines in the bottom panel and nicely match the two plateaux we identified. Indeed, the measured reconnection rate of the secondary current sheet is strongly super-Alfvénic.

In the final stage of the evolution the secondary reconnection events drive the dynamics, with new plasmoids being ejected by super-Alfvénic outflows and feeding the huge plasmoids generated by the first reconnection event. Eventually, the whole current sheet is disrupted.

5 Secondary ideal tearing instabilities

The study of the formation of secondary current sheets and their disruption by the onset of secondary tearing instabilities is obviously very important, since the observed secondary reconnection events have super-Alfvénic growth rates and drive the final evolution of the whole system. The aim of this section is to further characterize the spontaneously formed secondary current sheets before their evolution toward the final breakup. We will show that, in a dynamically evolving plasma environment, it is possible to form current sheets in local (and provisional) equilibrium which then evolve on small and local temporal and spatial scales in an explosive way. It is worth nothing that such substructures are within a global structure (the primary current sheet) which, on the contrary, is out of equilibrium and has already evolved in a highly turbulent state. Therefore, in a broader context results of this section have potential implications for what concerns the dynamics of turbulent systems.

Although morphologically different, the behavior of the evolution of all Hall-MHD runs is qualitatively similar to that of the MHD ones, the growth rates being larger and the final stage more violent for $\eta_{H}\neq 0$ , thus in the present section we focus only on purely MHD simulations.

We have already shown that secondary current sheets naturally form between consecutive plasmoids at the beginning of the nonlinear phase, with an approximate thickness which is $10\%$ to $15\%$ of the initial thickness $a$ (see Fig. 3). These current sheets further thin on a timescale of a couple of Alfvén times, until they reach a critical aspect ratio and become unstable to a secondary tearing instability. The formation of these secondary events is spontaneous, without any prior imposition on their equilibrium or their aspect ratio, therefore it is very interesting to characterize their evolution and the conditions under which the secondary tearing instabilities are triggered.

To that purpose, for a given simulation, we identified the region where a secondary current sheet had formed. Then we measured the position of its center, $y_{\mathrm{cs}}$, and we calculated its length $2L^{*}$ by measuring the full-width-half-maximum of the current density profile $J_{z}$ along $y$ at $x=0$, after the background current of the primary sheet had been subtracted. Moreover, by assuming a standard profile of the form

[TABLE]

and by performing a least square fit, we obtained the half thickness $a^{*}$ and the asymptotic magnetic field $\mathrm{B}_{0}^{*}$ of the secondary current sheet. Figure 5 shows an example of the fit, performed at the center of the secondary current sheet of the top-left panel of Fig. 3. The local Lundquist number is then found as $S_{a}^{*}=a^{*}\mathrm{B}_{0}^{*}/\eta$ (the density in the secondary sheets increases typically only by about $1\%$ of $\rho_{0}=1$, therefore we can safely identify $\mathrm{B}_{0}^{*}$ with the local Alfvén speed). Notice that, as it will be discussed later, these dynamically formed current sheets are in a state of almost perfect pressure equilibrium, hence we do not expect that a $z$ component of the magnetic field is needed to balance the magnetic pressure in a force-free state. We thus deem that Eq. (26) represents the best shape for the magnetic field to be used as a fit for the secondary current sheets.

In order to provide a statistically significant measure of $a^{*}$ and $\mathrm{B}_{0}^{*}$, a separate fit of the $y$-component of the magnetic field, $B_{y}(x,y_{i})$, is performed for each grid coordinate $y_{i}$ in the range $[y_{\mathrm{cs}}-L^{*}/2,y_{\mathrm{cs}}+L^{*}/2]$ We further note that the standard deviation $\Delta a^{*}$ is larger than $a^{*}$ when the secondary current sheet is either in its nonlinear reconnection phase or it has not formed yet. Therefore, values $\Delta a^{*}/a^{*}\ll 1$ indicate the phase in which a well defined secondary current sheet is present, while values $\Delta a^{*}/a^{*}>1$ give a rather precise indication of when secondary reconnection events are about to disrupt it.

To track the evolution in time of the secondary current sheet, we performed the above fitting procedure for all the outputs of the simulation. Figure 6 shows the evolution of $a^{*}$ and $\mathrm{B}_{0}^{*}$ of the secondary current sheet already discussed in the MHD reference run. The value of $\Delta a^{*}/a^{*}$ is color coded, so to capture the formation of the secondary current sheet. In particular, the blue points indicate a well defined current sheet, with $\Delta a^{*}/a^{*}\ll 1$, whereas red points either denote the thinning of an X-point or the presence of a nonlinear secondary tearing instability .

At the beginning of the simulation, at $t=0$, the fit correctly gives the thickness $a$ and the amplitude $\mathrm{B}_{0}$ of the primary current sheet. As time proceeds, an X-point forms and then gets elongated due to the evacuation of nearby plasmoids. At about $t=7\,\tau_{A}$ a secondary current sheet has formed, since $\Delta a^{*}\ll a^{*}$ there. In the subsequent evolution, the current sheet further thins, but keeping a constant local Lundquist number, in this case $S_{a}^{*}\simeq 270$ (the dashed line in the figure). This happens because the magnetic field is kept frozen in the plasma inside the current sheet, since the diffusion time is much larger than the time scales of the thinning process. At $t=9~{}\tau_{A}$ a secondary linear tearing instability starts to develop inside the current sheet and the thinning stops concurrently.

The configuration of the secondary current sheet at this time clearly shows an almost perfect pressure equilibrium, as shown by the profiles of Fig. 7. Moreover, both an inflow perpendicular to the sheet and an outflow along its main direction are present, naturally formed because of the evacuation and merging of the plasmoids in the evolution of the primary reconnection process. As expected, the inflow is very weak , while the outflow peaks at roughly half of the local Alfvén speed. The secondary linear tearing instability that we have just described appears to be triggered by perturbations in the magnetic field with an amplitude of about $1\%$ with respect to $\mathrm{B}_{0}^{*}$, hence it is bound to develop very rapidly and we actually witness the disruption of the secondary current sheet in less than an Alfvénic time.

The dynamics is qualitatively similar in all the MHD and Hall-MHD simulations we performed. In order to assess the scaling of the these secondary tearing instabilities with the local Lundquist number $S^{*}=L^{*}\mathrm{B}_{0}^{*}/\eta$, here defined using the half length $L^{*}$ of the secondary current sheet as characteristic scale, we performed the same analysis on additional In Figure 8 (left panel) we report, for each simulation, the aspect ratio $a^{*}/L^{*}$ of the secondary current sheet, calculated at the time when the thinning stops (different for each simulation), against $S^{*}$. The plot clearly shows that the scaling is consistent with that characteristic of the ideal tearing ($a^{*}/L^{*}\propto S^{*-1/3}$), although the SP scaling seems to be more appropriate for the lowest values of $S^{*}$. This may suggest the existence of two different regimes, in agreement with the findings of Huang et al. (2017). This scenario is also confirmed by looking at the reconnection rate $\gamma$ of the secondary current sheets. Figure 8 ( panel) shows that, once rescaled to the local Alfvén time $\tau_{A}^{*}=L^{*}/c_{A}^{*}$, the growth rate is compatible with the value $\gamma\simeq 0.63$ of the ideal tearing instability, with the exception of few points, that seem to be more compatible with the SP scaling . Note, however, that here the local Lunquist number is close to the threshold minimum value of $10^{4}$ requested to allow super-tearing modes (see Shi et al., 2018, for an exploration of lower values). Moreover, the agreement with the critical scenario of the ideal tearing seems to be improving with increasing $S^{*}$, as expected, since we are moving toward the asymptotic regime ($S^{*}\rightarrow\infty$).

In this section we have demonstrated that, in the evolution of the nonlinear phase of MHD (and similarly for Hall-MHD) reconnection, secondary events occur inside the reconnecting sheet, which spontaneously adjust so to reach an ideal tearing regime: a local (inverse) aspect ratio of the secondary current sheet $a^{*}/L^{*}\sim S^{*-1/3}$ and a local growth rate of the linear tearing fastest mode $\gamma_{\text{rec}}\,\tau_{A}\simeq 0.63$, independent of the local Lundquist number when $S\to\infty$.

6 Role of the Hall term

Even though for values of $d_{i}/\delta\lesssim 1$ the Hall term does not affect the reconnection rates of the primary current sheet, as we have already shown in Section 3, the changes in the secondary reconnection events are substantial. Figure 9 shows a plot of $\gamma_{\mathrm{rec}}$ for the primary and the secondary reconnection events for simulations with $S=10^{6}$ and different values of the ratio $d_{i}/\delta$ (Runs 8N-11N of Table 1). As we can see, the reconnection rate of the secondary events (red lines) increases by $17\%$ for a value $d_{i}/\delta=0.2$, and almost doubles for $d_{i}/\delta=1.4$. As a consequence, the evolution of the overall reconnection process is faster, the growth rate being up to five times the one of the MHD primary instability, and leads to the disruption of the current sheet in a correspondingly shorter time.

Moving to larger ratios $d_{i}/\delta$ (the case of $d_{i}/\delta=3$ is shown in the plot) we enter the strong Hall regime: the primary reconnection events become more and more violent, and the formation of secondary current sheets seems to be inhibited. The cause for the suppression of the secondary reconnection events in the strong Hall regime is unknown, however we can identify three possible explanations. The first possibility is that, for $d_{i}\gg\delta$, the linear growth rate of the secondary Hall tearing instability is so fast that it disrupts any forming current sheet before it can sufficiently grow to become dynamically important. A second possibility is that the geometrical configuration (i.e., the quadrupolar structure) of a X-point in the Hall regime prevents the formation of a secondary current sheet. Alternatively, the numerical resolution employed here may not be sufficient to reproduce the dynamics of the secondary reconnection events.

7 Conclusions

In this work we have presented a detailed study of the ideal tearing instability of thin current sheets in MHD and Hall-MHD plasmas, carried out by means of 2D compressible and fully nonlinear numerical simulations, along the same lines of Landi et al. (2015) and Del Zanna et al. (2016b). Our results confirm that magnetic reconnection via the ideal tearing instability is indeed an efficient mechanism of energy conversion, which is as fast as the ideal Alfvén timescales in MHD, and even faster in the Hall regime.

In the MHD regime, after the ideal tearing instability saturated and the nonlinear phase has begun, we observed the onset of secondary reconnection events in newly formed current sheets, thinner by one order of magnitude than the initial current sheet. These secondary tearing instabilities are strongly super-Alfvénic, with reconnection rates $\gamma_{\mathrm{rec}}\tau_{A}\simeq 2.3$, i.e., five times faster than the main instability one. The net result is a much more violent reconnection process and a speed up in the disruption of the current sheet. Moreover, numerical simulations in the Hall-MHD regime, performed with increasing values of $\eta_{H}=d_{i}/L$, have shown that even though the Hall effect is negligible in the linear phase , it considerably affects the secondary reconnection events in the nonlinear phase, by increasing the reconnection rates up to about $100\%$ (for $\eta_{H}=0.0014$, corresponding to $d_{i}/\delta=1.4$ for our reference simulation with $S=10^{6}$) with respect to the pure MHD case and about ten times the reconnection rate of the linear phase.

At higher values of $d_{i}/\delta$, the formation of secondary current sheets is not observed.

Particular attention has been devoted to the study of the conditions under which the secondary instability takes place. Previous studies already identified and highlighted the properties of ideal tearing instabilities triggered in secondary current sheets, that had formed either in presence of an artificially induced collapse (Tenerani et al., 2015b) or spontaneously (Landi et al., 2015, 2017) from the primary current sheet. Here we have quantitatively demonstrated for the first time that the new substructures, namely the thinning secondary current sheets formed among nearby X-points, spontaneously tend to the critical aspect ratio proper of the ideal tearing, $a^{*}/L^{*}\sim S^{*-1/3}$ for high $S^{*}$, this time calculated on the local spatial scales. In this phase the local Lundquist number remains constant and the sheet is in pressure balance with the external medium. Then the secondary instability fully develops, on timescales approaching the expected value $\gamma_{\text{rec}}\tau_{A}^{*}\simeq 0.63$, here using the shorter local value of $\tau_{A}^{*}\ll\tau_{A}^{\phantom{*}}$, thus on super-Alfvénic global timescales. This scenario has been investigated by performing several simulations varying the (global) Lundquist number in the range $S=10^{5}-10^{8}$: the ideal scaling for the locally formed secondary current sheets is retrieved for high $S^{*}(>3\cdot 10^{4})$, and the match with the asymptotic value for the instability growth rate improves with increasing values of $S$, as expected. Instead, for moderate low $S^{*}(\lesssim 3\cdot 10^{4})$ a regime compatible with a SP scaling is observed.

takes place at higher values of the Lundquist number (between $10^{5}$ and $10^{6}$). Moreover, they find a growth rate for the tearing instability which is much faster than 0.63, but still with the correct scalings in the two regimes. We believe that this discrepancy is only apparent, and it is due to the different normalization they used. In fact, in HU17 the current sheets form with a average half length $L^{*}=0.25$ (using our notation), which is correctly employed in the definition of $S^{*}=L^{*}c_{A}^{*}/\eta$. However, the Alfvén speed is set to the global one ($c_{A}^{*}=c_{A}^{\phantom{*}}=1$), without a strict check of its value, unlike done in this work. This may lead to overestimating $S^{*}$. Moreover, the measured growth rates (see Table 1 of HU17) are normalized with respect to a global Alfvén time $\tau_{A}=1$, and not with respect to the correct value $\tau_{A}^{*}=L^{*}/c_{A}^{*}$, that would be already four times smaller by using $L^{*}=0.25$.

Our detailed analysis on the reconnection dynamics of the secondary current sheets validates some predictions of Tenerani et al. (2015b). For instance, we retrieved a thickness of the secondary current sheet that roughly corresponds to the width of the inner resistive layer of the primary current sheet $a^{*}\simeq\delta\sim S^{1/2}$. Moreover, in the MHD case with $S=10^{6}$, our measured value of the reconnection rate $\gamma^{*}\tau_{A}\simeq 2.3$ matches their theoretical estimate for the linear growth rate of the secondary tearing instability.

Outcomes of this work have potential applications for explaining the explosive events in the strongly magnetized space and astrophysical plasmas mentioned in the introduction, and further extend results of previous works of recursive magnetic reconnection (Shibata & Tanuma, 2001; Tenerani et al., 2015b; Singh et al., 2019). Moreover, the interplay between fast reconnection and turbulence can be crucial, as predicted by reconnection-mediated turbulence models (Boldyrev & Perez, 2012; Loureiro & Boldyrev, 2017; Mallet et al., 2017) and by recent numerical simulations of the solar wind plasma retaining kinetic effects (Franci et al., 2017), where the role of reconnection in driving the turbulent cascade at sub-ion scales through the destabilization of current sheets of thickness $a\simeq d_{i}$ is established. However, a recent study by Papini et al. (2019) has demonstrated that the role of reconnection in shaping the spectrum of solar wind turbulence, including the change of slope at the ion inertial length $d_{i}$, may be captured even without resorting to (hybrid) particle-in-cell simulations, by just retaining the Hall effect within a macroscopic MHD description, as in the present study.

Acknowledgements

The authors wish to acknowledge valuable exchanges of ideas at T. Tullio. EP thanks Luca Franci and Daniele Del Sarto for fruitful discussion. SL and LDZ acknowledge support from the PRIN-MIUR project prot. 2015L5EE2Y Multi-scale simulations of high-energy astrophysical plasmas. This research was conducted with high performance computing (HPC) resources provided by the CINECA ISCRA initiative (grant HP10C2EARF and HP10B2DRR4).