The ideal tearing mode: theory and resistive MHD simulations
Classical MHD reconnection theories, both the stationary Sweet-Parker model and the tearing instability, are known to provide rates which are too slow to explain the observations. However, a recent analysis has shown that there exists a critical threshold on current sheet’s thickness, namely , beyond which the tearing modes evolve on fast macroscopic Alfvénic timescales, provided the Lunquist number S is high enough, as invariably found in solar and astrophysical plasmas. Therefore, the classical Sweet-Parker scenario, for which the diffusive region scales as and thus can be up to times thinner than the critical value, is likely to be never realized in nature, as the current sheet itself disrupts in the elongation process. We present here two-dimensional, compressible, resistive MHD simulations, with S ranging from to , that fully confirm the linear analysis. Moreover, we show that a secondary plasmoid instability always occurs when the same critical scaling is reached on the local, smaller scale, leading to a cascading explosive process, reminiscent of the flaring activity.
1 Introduction: from slow to ideal reconnection
Magnetic reconnection is thought to be the primary mechanism providing fast energy release, readily channeled into heat and particle acceleration, in astrophysical and laboratory magnetically dominated plasmas. Within the macroscopic regime of resistive magnetohydrodynamics (MHD), however, classical reconnection models predict timescales, in highly conducting plasmas, which are too slow to explain bursty phenomena such as solar flares in the corona or tokamak disruptions. In the following we summarize the main results of classical theories and the latest developments in this subject.
Given a plasma flow pattern , the magnetic field evolution is governed within resistive MHD by the induction equation
where is the magnetic diffusivity, supposed to be a constant. A simple dimensional analysis allows one to separate the characteristic timescales of fluid advection and Ohmic diffusion
where is the Lundquist number, that is the magnetic Reynolds number defined via the macroscopic quantities and , the Alfvén speed. Astrophysical tenuous plasmas are characterized by very large values of , say up to , so that diffusion occurs on very long timescales and only where currents are strong, that is in current sheets.
The Sweet-Parker model (hereafter SP) of two-dimensional, steady, incompressible reconnection (Sweet:1958, ; Parker:1957, ) predicts, for an aspect ratio ( is the current sheet length or breadth, identified with the macroscopic scale, and its width), a rate of reconnected flux
leading to time scales , far too slow to explain solar flares with s.
The linear stability of a current sheet of infinite length was investigated by Furth:1963 (), who discovered that the equilibrium is unstable to the tearing mode, leading to the formation of X-points and plasmoids during the reconnection process. If measured on top of the only available scale, the sheet half width (from now on starred quantities refer to rather than , so they indicate local properties rather than global), the instability growth rate is, again, far too slow
where is the wave number at which the instability peaks. The scaling for is exactly the same found for the stationary rate in the SP model, so once more we would need extremely small scales to explain the observations.
In order to make a step forward, one could also investigate the behaviour of the extremely thin current sheets predicted by the SP stationary theory. As first demonstrated by means of 2D MHD simulations Biskamp:1986 (), reconnecting SP-like sites become unstable once the Lundquist number exceeds a critical value of order , and are subject to fast tearing modes and plasmoid formation when their aspect ratio becomes large enough, also increasing the local reconnection rate. Recent detailed linear analyses and simulations have confirmed these findings (Loureiro:2007, ; Lapenta:2008, ; Samtaney:2009, ; Bhattacharjee:2009, ; Cassak:2009a, ; Huang:2010a, ; Uzdensky:2010, ; Cassak:2012, ; Loureiro:2015, ). In particular, the SP current sheet, of inverse aspect ratio , in the presence of the typical inflow/outflow pattern characterizing steady reconnection, was shown to be tearing unstable with growth rate and wave number of maximal instability
However, the existence of instabilities with growth rates scaling as a positive power of poses severe conceptual problems, since the ideal limit, corresponding to would lead to infinitely fast instabilities, while it is well known that in ideal MHD reconnection is impossible. Moreover, the positive power of also for indicates an explosive behaviour of the tearing mode, with the formation of a large number of plasmoids. From these paradoxical findings it should be clear that, since SP current sheets are violently unstable in the ideal MHD limit , this could only imply that such extremely elongated equilibrium structures with ( in corona) simply cannot form in nature.
This issue was resolved by Pucci:2014 () (PV hereafter), who studied the stability of current sheets with generic inverse aspect ratios . The authors showed that a critical exponent separates current sheets subject to slow instabilities, with growth rates scaling as a negative power of , from the unphysical fast instabilities scaling as a positive power of . For this generic dependence the growth rate is
thus, there is a critical value for what they called the ideal tearing mode, given by
for which the growth rate of the fastest reconnecting mode becomes independent on the Lundquist number. Notice that for the threshold is 100 times larger than the SP one. Thus, in this limit reconnection starts to occur on ideal timescales, and the SP configuration is never established in the thinning process. This scenario should also survive in the presence of moderate viscosity (Tenerani:2015, ).
Going in deeper details, consider a current sheet defined by the equilibrium magnetic field (for example the Harris sheet), the usual tearing mode linear analysis for fluctuations in an incompressible medium leads to the system (Furth:1963, )
where and are the linear fluctuations across the sheet (the other components are retrieved by ). When letting and with , the dispersion relation for varying clearly shows curves with an ideal asymptotic limit for , as can be seen in figure 1, taken from PV.
The wave number corresponding to the maximum growth rates is seen to decrease with as
as expected. However, once the macroscopic normalization against the current sheet length is recovered, we do expect a number of islands actually increasing with
thus even this ideal tearing mode leads to the disruption of the current sheet in a large number of islands, well before the SP sheet configuration (thinner than this one) can be reached.
2 Numerical simulations: setup
In order to prove the analytical theory by PV, Landi:2015 () performed 2D resistive MHD simulations of the same scenario, that is they studied numerically the stability of a PV current sheet with , subject to initially small velocity perturbations. Here we summarize and report the main results.
The compressible, resistive MHD equations in the form
are integrated numerically, where is the Lundquist number defined above, is the adiabatic index, and other quantities retain their obvious meaning. Physical quantities are normalized using Alfvénic units, namely a characteristic length scale , a characteristic density , and a characteristic magnetic field strength (the background values measured far from the current sheet). Velocities are then expressed in terms of the Alfvén speed , time in terms of , the fluid pressure in terms of . Note that we are using the energy equation (14) written for the normalized temperature , where we use as a reference value . With the given normalizations the Lundquist number is basically the inverse of the magnetic diffusivity, namely .
The initial condition at for our two-dimensional simulations of the tearing instability is a main magnetic field along the direction , with and everywhere, asymptotic pressure (and temperature) and field magnitude far from , the current sheet center where . The generic equilibrium can be conveniently described as
where for the Harris sheet case with and pressure equilibrium (PE hereafter), with varying to preserve constant, and for the force-free equilibrium (FFE hereafter) with magnetic field rotating across the sheet. Intermediate cases of mixed fluid/magnetic pressure equilibrium with are also possible.
The compressible, resistive MHD equations (12-15) are solved in a rectangular numerical box with resolution and respectively. In the -direction, in order to resolve the steep gradients inside the current sheet using a reasonable number of grid points, we limit our domain to a few times the current sheet thickness, i.e. we set , with : this is a good compromise between the high resolution required inside the current sheet and the need to have boundaries sufficiently far from the reconnecting region. Along the -direction the length is chosen in order to resolve for the fastest growing modes of the instability (see (Landi:2015, ) for details). Finally, incompressible velocity perturbations of amplitude are applied at in order to trigger the instability.
The numerical simulations are performed by integrating equations (12-15) with an MHD code developed by our group. Along the current sheet, where periodicity is assumed, spatial integration is performed by using pseudo-spectral methods, while in the -direction integration is performed by the use of a fourth-order scheme based on compact finite-differences (Lele:1992, ). The boundary conditions in the non-periodic direction are treated with the method of projected characteristics, here assuming non-reflecting boundary conditions. Time integration is performed using a third-order Runge-Kutta method. Details of the code are described in Landi:2005 (). The resolution is adapted to the Lundquist number we use: for and we choose and , while for the number of cells in the direction is increased up to . In the periodic direction we use for the single-mode runs (in order to reduce the computational costs as many simulations are required to reproduce the instability dispersion relation curves), while we take in the nonlinear reference simulation. We have verified that this relatively low resolution along the periodic direction is adequate, due the extreme accuracy of Fourier methods and the rather smooth gradients observed in the direction. In spite of the relatively high values of S, in addition to the instability evolution, the initial equilibrium diffuses on timescales which, although long compared to the instability one, are still sufficient to underestimate the growth rates of linear modes. To avoid this problem, only for the single mode linear analysis, the diffusion term of the initial equilibrium (not that related to the perturbations, obviously needed for the tearing instability evolution) is removed from the induction equation at all times, as explained in Landi:2008 ().
3 Numerical simulations: single mode analysis
A first set of simulations of the tearing instability in current sheets with is performed to confirm the expected linear behavior presented in Section 1. We choose to test the two limits of our initial equilibria for the current sheet, namely PE () and FFE (). Moreover, we investigate both cases with () and (), where is the asymptotic plasma beta in equation (16), even if no differences are expected at the linear level from the incompressible, analytical analysis. Finally, three different values of the Lundquist number are tested here, namely , , and , for a total of 12 sets of simulations of the linear phase of the tearing instability, with the aim of reproducing the expected dispersion relations numericall, as shown in figure 2. As anticipated in the previous section, for each value of we vary while always selecting a single mode . The growth rate of the instability is computed by measuring the -averaged amplitude of the component of the perturbed magnetic field ( at the initial time).
The first thing to notice by inspecting the computed dispersion relations is that, as predicted by the classical linear theory, for each value of the curves have a maximum at a given , the peak location decreasing in as increases. The growth rate of the instability (normalized to the inverse of the large-scale Alfvén time ) has peaks ranging from for to for . In general we find that the simulations with (bottom panel) are more precise in matching the analytical results than those with (top panel), since a large beta is a condition closer to incompressibility (formally corresponding to an infinite value for the sound speed). Moreover, we find that simulations of the FFE scenario (squares) yield higher values of the growth rates, closer to the results from the analytic theory, as compared to those employing the PE settings (crosses): this is probably due to the fact that the purely force-free equilibrium leads to intrinsically less compressible fluctuations. Finally, rather large discrepancies are observed for small scales (large values of ), especially in the PE case.
In figure 3 we plot the profiles of the perturbations ( is determined by ), and , all normalized to their respective maximum, across the current sheet in the direction. In the top panels we show the analytical results, that is the eigenmodes of the linear analysis (here the PV calculations have been recomputed by imposing for ), and in the lower panels we report the numerical solutions for a simulation in the FFE scenario with and , at a given time of the linear evolution of the instability. In order to recover the theoretical eigenmodes, velocity and magnetic field perturbations are shown with a shift in , as expected. Notice the steep gradients arising within the current sheet (), where a high resolution is needed to resolve the small scales developed during the instability evolution. The eigenmodes are very well reproduced: the magnetic field perturbations are identical to the analytical expected ones, while in the velocity perturbations the only major difference is due to the non-reflecting free-outflow boundary conditions imposed, that do not force at as in the analytical solutions and result in a slightly different profile even in the vicinity of the reconnecting region.
4 Numerical simulations: fully nonlinear case
Now we describe the nonlinear stages of the evolution of the tearing instability. Since we are interested in its late development, where interaction and merging of plasmoids is expected, we trigger the instability by selecting an initial spectrum of modes, rather than a single one as in the previous set of simulations, and we choose a maximum mode number . We also choose and , so the modes with are all excited. From the theoretical curves in the previous section we expect mainly a competition between modes and as the fast growing ones. As a reference run, we analyze the instability of a force-free equilibrium with constant temperature (FFE, ), and we select the case with . This combination was shown to provide a linear phase which is the closest to the analytical expectations (see figure 2). The resolution employed for this run is , which is very high if one consider that the code employs high-order methods (compact finite-differences along and Fourier transforms along , where periodical boundary conditions apply).
The complete evolution is shown in figure 4, where 2D snapshots for selected times are provided, of the region . Note that the figure has different scales in and , the axis is normalized against , the axis against , with ). Panels on the left show the current component, whereas panels on the right show the distorted field lines superimposed on a map of temperature . It is easy to recognize the end of the linear phase of the tearing instability (top row), with the dominant mode still clearly apparent. When the tearing instability growth is over, the nonlinear phase sets in leading to further reconnection events and eventually to island coalescence, departing from the dominant phase with . First, due to the attraction of current concentrations of the same sign, the X-points elongate and stretch along the direction. This process leads to a strong increase of the electric current concentrations and thus to the formation of new, elongated current sheets (see the second row, details of this phase will be discussed further on). Beyond (third row) the evolution has become fully nonlinear and we clearly observe the process leading to the creation of a single, large magnetic island as arising from coalescence. The situation is very dynamic, especially near the major X-point where explosive expulsion of smaller and smaller islands is observed, typical of the plasmoid instability. These islands then move towards the largest one, which is continually fed and thus further increases its size in a sort of inverse cascade eventually leading to the largest mode. Notice also the temperature enhancement at the reconnection sites. At time (bottom row) both the current concentration and the plasma temperature around the X-point are so high that we need to saturate their values in order to retain an appropriate dynamical range in the color bar. The initial macroscopic current sheet is basically disrupted in a series of highly dynamical features. The current and temperature enhancements are stronger at the X-point and at the boundaries of the magnetic islands. Moreover, from the major reconnecting site we clearly see the production of magnetosonic waves, which propagate and soon steepen into shocks.
Let us now investigate in more detail the situation right after the end of the linear phase of the tearing instability, when the islands coalescence is about to set in and the local current sheets have just formed and start to further evolve. In figure 5 (left panel) we show the zoom around the X-point (which is just about to develop) near at , the first row of figure 4. Here we display the electric current by using the same scaling in both and directions, normalized against the macroscopic current sheet’s width , thus the position of the reconnection zone is now expressed as . The local current sheet has just formed as the result of a stretching process in the direction (and shrinking across the other direction), as a typical output of the nonlinear phase of the tearing instability. We are now in the phase in which, in the center of the current sheet, reconnection is about to take place, leading to a topology change in the magnetic structure and to the disruption of the current sheet itself, initially into two smaller strips. It is very interesting to measure the aspect ratio of this reconnection site. From the white box in the figure, defined by the rectangular region where the component has values which are roughly half those of the central peak, we estimate , where we have used an asterisk to indicate the local values and to differentiate with respect to the macroscopic ones (we also recall that in the whole paper we identify as the half width of a current sheet). Additional simulations for and , reported in the second and third panels, lead to local aspect ratios of and , respectively.
If we now compare these numbers with , trying to find a value of that fits the data best, it is easy to see that the value is a very good guess. Here is the local Lundquist number, which is obviously smaller than the macroscopic one, due to the much smaller length of the local current sheet (to be measured in each case). Therefore, based on our very limited data set, we derive the scaling
where , that is of order unity, as expected. In figure 6 we report the numerical results, assuming a error on both and , for the sequence of increasing (macroscopic) Lundquist numbers , , and . The scaling in equation (17) is over-plotted, as a solid line, for the average value of .
In the present paper we have summarized the main results of the analytical analysis by Pucci:2014 () and of the numerical simulations by Landi:2015 (), who studied, by means of compressible, resistive MHD simulations, the linear and nonlinear stages of the tearing instability of a current sheet with inverse aspect ratio , where is the Lundquist number measured on the macroscopic scale (the current sheet length) and the asymptotic Alfvén speed. Our results confirm on the one hand the linear analysis of PV of the ideal tearing mode, leading to extremely fast growth rates with when is sufficiently large, and on the other hand the nonlinear simulations show that the evolution follows what appears to be a quasi-self-similar path, with subsequent collapse, current sheet thinning, elongation, and destabilization, starting from the X-points formed in the original sheet. As scales become smaller, and the local Lundquist numbers decrease, the dynamical time-scales decrease too, leading to explosive behavior.
These findings are very important, in our opinion. For the first time we clearly see in simulations that, even in the nonlinear stages of the tearing instability, the new current sheets that form locally become unstable when the inverse aspect ratio of these structures reaches the critical threshold of (where the star indicates the local, smaller scale), precisely the same limit found by PV for the fast reconnection of the initial, macroscopic current sheet. After that, a new ideal tearing instability starts, with time occurring on a faster timescale , since typically . Furthermore, when smaller and smaller scales are produced nonlinearly, as observed at the time proceeds, each time the newly formed local current sheets elongate and reach their own critical value, that corresponding to equation (17), faster and faster reconnection will arise producing a cascading, accelerating process: this, we believe, is the real nature of the plasmoid (or super-tearing) instability.
There are many space and astrophysical applications of the ideal tearing, including geomagnetic storms in our planet, coronal heating and coronal mass ejections from the Sun. As future work we plan to move to 3D simulations (Bettarini:2009, ; Landi:2012, ) and to investigate the tearing instability of thin current sheets within relativistic MHD (Del-Zanna:2007, ; Bucciantini:2013, ; Del-Zanna:2014, ). Fast reconnection in magnetically dominated plasmas is of paramount importance in high-energy astrophysics too (Kagan:2015, ), invariably invoked in models for magnetar flares, acceleration of Poynting-dominated jets, and dissipation in pulsar winds (Tavani:2013, ; Porth:2013, ; Olmi:2015, ).
M.V. was supported by the NASA Solar Probe Plus Observatory Scientist grant. The research leading to these results has received funding from the European Commissions Seventh Framework Programme (FP7/2007-2013) under the grant agreement SHOCK (project number 284515).
- (1) Sweet P A 1958 Electromagnetic Phenomena in Cosmical Physics (IAU Symposium vol 6) ed Lehnert B p 123
- (2) Parker E N 1957 J. Geophys. Res. 62 509–520
- (3) Furth H P, Killeen J and Rosenbluth M N 1963 Physics of Fluids 6 459–484
- (4) Biskamp D 1986 Physics of Fluids 29 1520–1531
- (5) Loureiro N F, Schekochihin A A and Cowley S C 2007 Physics of Plasmas 14 100703 (Preprint astro-ph/0703631)
- (6) Lapenta G 2008 Physical Review Letters 100 235001 (Preprint 0805.0426)
- (7) Samtaney R, Loureiro N F, Uzdensky D A, Schekochihin A A and Cowley S C 2009 Physical Review Letters 103 105004 (Preprint 0903.0542)
- (8) Bhattacharjee A, Huang Y M, Yang H and Rogers B 2009 Physics of Plasmas 16 112102 (Preprint 0906.5599)
- (9) Cassak P A, Shay M A and Drake J F 2009 Physics of Plasmas 16 120702
- (10) Huang Y M and Bhattacharjee A 2010 Physics of Plasmas 17 062104 (Preprint 1003.5951)
- (11) Uzdensky D A, Loureiro N F and Schekochihin A A 2010 Physical Review Letters 105 235002 (Preprint 1008.3330)
- (12) Cassak P A and Shay M A 2012 Space Sci. Rev. 172 283–302
- (13) Loureiro N F and Uzdensky D A 2015 ArXiv e-prints (Preprint 1507.07756)
- (14) Pucci F and Velli M 2014 ApJLett 780 L19
- (15) Tenerani A, Rappazzo A F, Velli M and Pucci F 2015 ApJ 801 145 (Preprint 1412.0047)
- (16) Landi S, Del Zanna L, Papini E, Pucci F and Velli M 2015 ApJ 806 131 (Preprint 1504.07036)
- (17) Lele S K 1992 Journal of Computational Physics 103 16–42
- (18) Landi S, Velli M and Einaudi G 2005 ApJ 624 392–401
- (19) Landi S, Londrillo P, Velli M and Bettarini L 2008 Physics of Plasmas 15 012302
- (20) Bettarini L, Landi S, Velli M and Londrillo P 2009 Physics of Plasmas 16 062302 (Preprint 0906.5383)
- (21) Landi S and Bettarini L 2012 Space Sci. Rev. 172 253–269
- (22) Del Zanna L, Zanotti O, Bucciantini N and Londrillo P 2007 A&A 473 11–30 (Preprint 0704.3206)
- (23) Bucciantini N and Del Zanna L 2013 MNRAS 428 71–85 (Preprint 1205.2951)
- (24) Del Zanna L, Bugli M and Bucciantini N 2014 Astronomical Society of the Pacific Conference Series (Astronomical Society of the Pacific Conference Series vol 488) ed Pogorelov N V, Audit E and Zank G P p 217 (Preprint 1401.3223)
- (25) Kagan D, Sironi L, Cerutti B and Giannios D 2015 Space Sci. Rev. (Preprint 1412.2451)
- (26) Tavani M 2013 Nuclear Physics B Proceedings Supplements 243 131–140
- (27) Porth O, Komissarov S S and Keppens R 2013 MNRAS 431 L48–L52 (Preprint 1212.1382)
- (28) Olmi B, Del Zanna L, Amato E and Bucciantini N 2015 MNRAS 449 3149–3159 (Preprint 1502.06394)