Nonlinear Explosive Magnetic Reconnection in a Collisionless System

Recently, magnetic reconnection has received great attention in various plasma environments of the Earth’s magnetosphere, laboratory plasmas, and astrophysical objects such as pulsar magnetospheres and magnetars [e.g., 1, 2, 3, 4, 5]. In these hot and dilute magneto-active plasmas, the magnetic energy stored in a plasma sheet with an anti-parallel magnetic field component can be efficiently converted into not only the bulk flow energy but also the energies of plasma heating and nonthermal particle acceleration. The reconnection paradigm has been established as a powerful magnetic energy dissipation process in the plasma universe [e.g., 6, 7, 8, 9, 10]. In addition, in the last couple of decades, our understanding of reconnection dynamics has gradually accumulated for a “thin” plasma sheet with a thickness of the order of the particle gyro-radius, and it is recognized that fast energy dissipation can occur for such a thin plasma sheet. However, the mechanism of fast energy dissipation of reconnection for a “thick” plasma sheet has remained an intriguing enigma for more than half a century. According to the linear theory of collisionless tearing modes, the growth rate $\gamma_{\rm l}\tau_{\rm A}$ is known to be proportional to $(r_{\rm g}/\lambda)^{3/2}$, where $r_{\rm g}$, $\lambda$, and $\tau_{A}$ are the gyro-radius of the thermal particle, the thickness of the current sheet, and the Alfvén transit time defined by $\lambda/v_{A}$, respectively. The linear theory suggests that the energy dissipation for a thick plasma sheet with $r_{\rm g}\ll\lambda$ is extremely slow.

To overcome the difficulty of the fast energy dissipation for a thick plasma sheet, many people have challenged the mechanisms of the explosive nonlinear growth of the reconnecting magnetic flux in the fluid approximation so far [e.g., 13, 14, 15, 16, 17, 18]. Although the standard resistive magnetohydrodynamics (MHD) model was known not to initiate a fast reconnection unless an anomalous resistivity is locally enhanced around the X-type neutral point [e.g. 11, 12], two-fluid/collisionless effects such as electron inertia, Hall current, ion gyro-radius have been argued to drive a fast reconnection. During the early time evolution, a rapid shrinkage of the current sheet can happen, and the thinning of the current sheet into such kinetic scales can lead to the explosive reconnection. It has been argued that the topological conservation of the Lagrangian invariant leads to the formation of the thin current sheet in collisionless reconnection [14, 17].

In addition to the shrinkage of the current sheet in collisionless reconnection, Galeev et al. [19] proposed another idea of the explosive reconnection mechanism by focusing on the time evolution of the collisionless inertia resistivity in association with the growth of the reconnecting magnetic field. They considered the enhancement of the inertia resistivity due to magnetized particles by the reconnecting magnetic field around the X-type neutral line, and estimated that the nonlinear growth rate $\gamma_{\rm G}\tau_{A}$ was of the order of $(r_{\rm g}/\lambda)^{1/2}b_{1}$, where $b_{1}=B_{1}/B_{0}$ is the ratio of the reconnection magnetic field $B_{1}$ to the magnetic field $B_{0}$ outside of the plasma sheet. This behavior has been demonstrated by particle-in-cell (PIC) simulation [20]. This nonlinear reconnection suggests not only much faster energy dissipation than in the linear stage, but also an explosive nature of growth, namely, $b_{1}(t)=b_{1}(t_{0})/(1-(t-t_{0})\gamma_{\rm l})$, where $t_{0}$ is the start time of the nonlinear reconnection. The nonlinear explosive growth, however, still depends on the initial plasma thickness of $r_{\rm g}/\lambda$, and the time scale is of the order of $\gamma_{\rm l}^{-1}$. Therefore, the fast reconnection issue for a thick current sheet remains unresolved.

In this paper, based on high-resolution two-dimensional PIC simulation, we propose a new mechanism of nonlinear reconnection in which the energy dissipation can occur in the time scale of the Alfvén transit time. The crucial factor is the enhancement of the inertia resistivity due to the effects of both thinning of the current sheet and the reconnecting magnetic field during reconnection. As the key process, we show that the current sheet shrinks quickly to a thickness of the gyro-radius in PIC simulations for a pair plasma. The theoretical model to explain the nonlinear explosive reconnection is also discussed.

Let us first study the time evolution of the tearing-mode instability using two-dimensional PIC simulation [21, 22, 23] with a periodic boundary in the $x$-direction and the conducting walls for the upper and lower boundaries at $y=\pm 8\lambda$, where $\lambda$ is the thickness of the initial plasma sheet. The total system size $L_{x}\times L_{y}$ is $108\lambda\times 16\lambda$ and the computational grid size is $5376\times 800$. The total number of particles is $1.2\times 10^{11}$ to calculate the high-$\beta$ plasma sheet as accurately as possible. This large number of particles is necessary for the linear growth of the tearing mode to overcome thermal fluctuations, the amplitudes of which are proportional to $1/\sqrt{n_{\rm p}}$, where $n_{\rm p}$ is the number of particles per grid cell.

For simplicity, we adopt the Harris solution [24] for the pair plasma with the same mass $m$ and temperature $T$, and we focus on a ubiquitous nonlinear reconnection mechanism in a collisionless plasma system. The plasma temperature $T=mc^{2}/10$ and the Alfvén speed $v_{\rm A}=B_{0}/\sqrt{8\pi mN_{0}}=0.447c$ are maintained for all simulation runs. The thermal velocity $v_{\rm th}$ is equal to $v_{\rm A}$. The magnetic field $B$ and plasma density $N$ are given by $B_{x}(y)=B_{0}\tanh(y/\lambda)$ and $N(y)=N_{0}\cosh^{-2}(y/\lambda)+N_{\rm b}$, respectively, where $N_{\rm b}$ is the background uniform plasma density. The background plasma density $N_{\rm b}$ is set to be $5\%$ of the maximum Harris plasma density $N_{0}$. In addition, we have the relationships of the pressure balance $B_{0}^{2}/8\pi=2N_{0}T$ and the force balance $2T/\lambda=|e|\beta B_{0}$, where $\beta c$ is the drift velocity. Three different cases of $\beta=0.05$, $0.1$, and $0.2$, which correspond to $r_{\rm g}/\lambda=0.11$, $0.22$, and $0.45$, are studied. Note that to date most simulation studies have investigated a thin plasma sheet such as $r_{\rm g}/\lambda\sim 0.5-2.0$. [e.g., 25, 26, 27, 28]

We add a small initial perturbation for the vector potential $\delta A_{z}$ to the Harris equilibrium, which is given by

where $\lambda_{x}=2\lambda$, $\lambda_{y}=\lambda/2$, $x_{0}=54\lambda$, and the initial amplitude of the reconnected magnetic field is set to be $\max(B_{y}/B_{0})=10^{-3}$ in the neutral sheet. The reason of the initial perturbation is to place the X-type point in the center of the simulation box and to save the computational CPU time. We have also studied the case without $\delta A_{z}(x,y)$, and observed the almost same time evolution as the case with $\delta A_{z}(x,y)$. Other remarks will be given in the section of Discussion and Summary. Although this perturbation mimics an X-type reconnection structure, the initial perturbation does not necessarily satisfy an eigenfunction for the linear tearing mode. Therefor, broadband electromagnetic waves are generated. To suppress these waves, a numerical filter outside the plasma sheet is used during the early phase with $t/\tau_{\rm A}<16$.

Figure 1 shows the time evolution of plasma density for $r_{\rm g}/\lambda=0.11$. Part of the simulation domain is depicted. Three typical stages of the time evolution of reconnection are shown by the color scales indicated by the bars at the right-hand side of the panels. The white lines are magnetic field lines calculated from the contour of the vector potential $A_{z}(x,y)$. The top panel is the early linear stage at $t/\tau_{\rm A}=290$, when the amplitude of the reconnecting magnetic field is gradually growing; the second panel is the nonlinear stage at $t/\tau_{\rm A}=519$, when thinning of the plasma sheet is clearly observed around the X-type point; and the third panel at $t/\tau_{\rm A}=533$ is the late nonlinear stage, when the elongation of the current sheet and the formation of several small plasmoids are observed. The time evolution is basically same as many previous reconnection simulations [e.g., 20, 21, 25, 26, 27, 28], but shown here is a thick current sheet in which the gyro-radius is much smaller than the thickness of the plasma sheet.

Note that other simulation runs for $r_{g}/\lambda=0.22$ and $0.45$ show almost the same time evolution as that of Figure 1 (not shown here). Namely, the plasma sheet gradually shrinks during the early time evolution, and after the formation of X-type neutral point at the center of the simulation box, the elongated thin current sheet is formed, where several small scale plasmoids are embedded. We find the small scale plasmoids in the elongated current sheet has a tendency to be easily formed for the thick current sheet for $r_{g}/\lambda=0.11$.

The top three panels in Figure 2 show the time evolution of the power spectra for the Fourier components of $|A_{z}(k_{x})|$ around $y=0$ as the function of wave number $k_{x}$. $A_{z}(k_{x})$ is obtained by

From the left to the right-hand panel, the simulation results for $r_{g}/\lambda=0.11$, $0.22$, and $0.45$ are depicted, respectively. The dark blue line represents the initial stage, the blue to green colors are the linear to the early nonlinear phase of the tearing mode/magnetic reconnection, and the red one is the nonlinear stage after the formation of the elongated thin current sheet. The color bar on the top indicates the corresponding time stage. We find that from the early time evolution with the dark blue to the green colors the vector potential $|A_{z}(k_{x})|$ can grow only for $0<k_{x}\lambda<1$, suggesting that the tearing mode is unstable for the modes with long wavelengths, and it has the maximum growth rate around $k_{x}\lambda=0.2-0.5$.

During the early nonlinear phase with the green to the yellow colors, the wave number with the maximum wave power of $|A_{z}(k_{x})|$ is gradually shifted into the short wave number region. The reddish lines correspond to the stage after the formation of the elongated thin current sheet around the X-type neutral point, and we find a spiky structure due to the elongated current sheet by analogy with the Fourier transformation for Heaviside function. The wave number of the first spiky notch $k_{\rm notch}$ is related to the length of the elongated current sheet $l$ by $k_{\rm notch}=2\pi/l$, and its high harmonics appear at $k=2\pi n/l$ with $n=2,3,..$ The wave number $k_{\rm notch}$ is shifted to a small regime in association of the elongation of the current sheet.

The bottom three panels in Figure 2 show the time histories of each Fourier components of $|A_{z}(k_{x})|$ as the function of time $t/\tau_{A}$. The data sets are same as the above three panels. The color lines represent the wave number, which is indicated on the top color bar. The reddish colors are the small wave numbers of $k_{x}\lambda\sim 0$, and the bluish colors are the large wave numbers of $k_{x}\lambda\sim 1$. From the left to the right-hand panel, three simulation runs for $r_{g}/\lambda=0.11$, $0.22$, and $0.45$ are respectively shown. If we focus on the time evolution of the Fourier mode with the maximum wave power shown by the yellowish/reddish color line, we find that there are two time stages of the linear growth phase and the rapid growth phase followed after the linear growth phase. The rapid growth phases appear around $t/\tau_{A}\sim 520$ (left), $230$ (middle), and $110$ (right), respectively. In the following sections, we focus on the second phase with the rapid nonlinear growth. Before discussing the rapid nonlinear phase, however, we analyze the linear phase obtained in our PIC simulations in the next section.

We shortly study the linear stage for the tearing instability/magnetic reconnection seen in Figure 2. Shown in Figure 3 is the linear growth rates calculated from the slopes of the time history of $|A_{z}(k_{x})|$ in the bottom panels in Figure 2 for three cases for $r_{g}/\lambda=0.11$ (red), $0.22$ (green), and $0.45$ (blue). In Figure 2, the slopes of the wave modes with $k\lambda\sim 0$ and $k\lambda\sim 1$ are more or less fluctuating compared with other wave numbers around $k\lambda=0.2-0.5$, because the particle-in-cell simulation has intrinsic “thermal” noises due to a finite number of particles in the grid cell. We made a least-square fitting for the slopes of $r_{g}/\lambda=0.11$, $0.22$, and $0.45$ by a function of $\exp(\gamma_{L}t)$ for the time interval of $t/\tau_{\rm A}=[10,150]$, $[10,70]$, and $[10,55]$, respectively. The closed circles are the linear growth rate $\gamma_{L}$ obtained in the simulation runs.

The dashed curves superposed on the simulation data points are the theoretical linear growth rates of Equation (14) in Appendix, which is given by

As the inertia conductivity may be an order of estimation due to the straight orbit of the meandering particle (see Appendix), we have multiplied a factor of $0.73$ on the right-hand side of the above equation. In addition, we assumed the background uniform plasma density $N_{b}$ in the initial condition, which may reduce the growth rate. From Figure 3, not only the dependence of the wave number $k\lambda$ but also the dependence of the the ratio of thickness $(r_{g}/\lambda)^{3/2}$ is consistent between the simulation results and the theory. Note that the factor of $k\lambda/(k\lambda+\sqrt{r_{g}/\lambda})$ comes from the correction of the non-adiabatic energy (see Appendix). The theoretical growth rate is applicable only for $\sqrt{r_{g}/\lambda}<k\lambda<1$, nevertheless the consistency in the small $k\lambda<\sqrt{r_{g}/\lambda}$ is good.

Although we found the linear growth rates is basically consistent with the theoretical prediction, the profile for the thick current sheet with $r_{g}/\lambda=0.11$ slightly differs from other two cases. In addition, we also see some discrepancy between the theory and the simulation around $k\lambda\sim 1$. As we have already mentioned, PIC simulations has a thermal fluctuation, which is proportional to $1/\sqrt{n_{p}}$, and it is not easy to demonstrate the linear growth phase for the wave modes with the relatively small growth rates.

The debate regarding fast energy dissipation has remained a fundamental issue in magnetic reconnection for a thick plasma sheet, because the linear growth rate is known to decrease as the thickness of the current sheet is increased. In this paper, we have argued that the nonlinear process of reconnection in association with the local enhancement of the inertia resistivity leads to fast energy dissipation, the rate of which is independent of the initial thickness of the plasma sheet. The local enhancement of the inertia resistivity around the X-type neutral point emerges by two effects of the shrinkage of the current sheet into the order of the particle gyro-radius and the magnetization of the non-adiabatic meandering particles by the reconnecting magnetic field $B_{y}$. In the master equation of the energy principle in Equation (1), the former effect appears in $d_{y}$, and the latter one is $d_{x}$ due to the kinetic effect.

The explosive nature of reconnection has been already studied mainly in the framework of fluid approximation so far, and our finding of the shrinkage of initially thick current sheets in PIC simulations is basically same, but the enhancement of inertia resistivity proposed by Galeev et al. [19] has not been included in those numerical studies in the fluid approximation. The constant nonlinear growth rate, which is independent of the initial thickness of current sheet, is realized by the help of the local enhancement of the inertia resistivity.

By integrating in time for both sides of Equation (4) by noting $\gamma_{\rm nl}=d(\log b_{1})/dt$, the amplitude increases explosively in time as

where $b_{1}(0)$ is the amplitude at the start time $t_{0}$ of the nonlinear reconnection. This equation has a singularity at $\Delta t/\tau_{A}=(t-t_{0})/\tau_{A}=(4k\lambda)/(b_{1}(0)\sqrt{\pi})$. As the coefficient of $(4k\lambda)/(b_{1}(0)\sqrt{\pi})$ is of the order of unity, the explosive nonlinear reconnection can occur in the Alfvén transit time within a factor of the order of unity. In fact, we find that the explosive growth in our simulation occurs in tens of the Alfvén transit time. During this explosive phase, however, the free energy on the left-hand side of Equation (1) is quickly exhausted, and thus the explosive reconnection may then cease. The nonlinear saturation is beyond the scope of this paper.

We would like to mention the initial perturbation of $\delta A_{z}(x,y)$. We discussed that the explosive reconnection occurs when $b_{1}>r_{g}k$, and the critical amplitudes of $b_{1}=B_{y}/B_{0}=(r_{g}/\lambda)(k_{\rm min}\lambda)$ in our simulation can be given as $6.4\times 10^{-3}$, $1.3\times 10^{-2}$, and $2.6\times 10^{-2}$ for $r_{g}/\lambda=0.11$, $0.22$ and $0.45$, respectively. On the other hand, the initial perturbation in the simulation was set to be ${\rm max}(B_{y}/B_{0})=10^{-3}$, which is smaller than the above critical amplitudes of the nonlinear reconnection. This parameter setting makes possible to demonstrate the time evolution from the linear stage. In fact, the initial X-type perturbation of $\delta A_{z}(x,y)$ has the maximum amplitude at $k_{\rm min}\lambda=0.0558$ in the Fourier space, which can be seen in the bottom purple/crimson lines for the top three panels in Figure 2, while the peak amplitudes in the linearly growing stage (blueish lines) are situated around $k\lambda=0.2\sim 0.5$. And then, in the nonlinear stage, the most unstable mode is shifted toward the longest wavelength. This confirms that the initial perturbation was set to be small so as to discuss the separation of two time stages of the linear and the nonlinear evolution.

We investigated the nonlinear process to drive a fast magnetic reconnection for a pair plasma, but we think that this nonlinear process is applicable to the case for an ion-electron plasma as well. During the early time evolution, the meandering electron around the X-type neutral point first starts to be magnetized, and if the electron scale current sheet embedded in the ion scale current is formed, it would be possible to observe the explosive reconnection in the electron scale. During the later time evolution, the major player of the energy dissipation will be switched to the ion inertia resistivity, and we may observed it in the ion scale as well. However, if the current sheet thickness becomes of the order of ion inertia length, the Hall current and the ambipolar electrostatic field become important [21, 30, 36]. The electron stabilization effect due to the magnetized electrons is another open question [37, 38, 39].

It would be important to mention the guide field $B_{z}$ effect. We discussed the collisionless reconnection for an anti-parallel magnetic field topology without $B_{z}$, in which the inertia resistivity can be provided by unmagnetized particles. For the guide field case, however, the inertia resistivity can be supplied by magnetized particles moving along the magnetic field line [40]. Therefore, our claim of the enhancement of the inertia resistivity due to the magnetization effect of the reconnecting magnetic field will not occur for a strong guide field case.

Our final comment is the nonlinear stage. In the late nonlinear stage after the thinning of the plasma sheet, we observed the formation of the elongated current sheet [e.g., 13, 41, 42, 43] in association with several small-scale plasmoids in Figure 1. Although the stage appears after the saturation of the reconnecting flux, it can lead to dynamic evolution of reconnection by launching magnetic islands and turbulent waves, and this stage is another important phase of reconnection discussed by many papers [e.g. 43, 36, 44, 45, 46, 47, 48, 49]. Our claim of the explosive reconnection can bridge the gap between the late nonlinear phase and the linear phase.

The data support the findings of this study are available from the corresponding author upon reasonable request.