2.1 Magnetic reconnection in tearing modes and formation of magnetic islands

A magnetic reconnection process is characterised by the formation of an $X$-point, where initially distinct magnetic lines have connected thanks to non-ideal effects that become important as the gradients of the magnetic field components, i.e. the components of the current density, are large enough. As a result of the magnetic reconnection event, the plasma system typically relaxes to a final state of lower ‘potential magnetic energy’, the diminished energy being converted into plasma kinetic and internal energy as well as into electron acceleration. During this process, the ‘magnetic topology’ in the plasma changes since the connection of fluid elements is globally modified – see figure 1.

Figure 1.

Cartoon schematising three steps of a simplified 2-D magnetic reconnection process that emphasises the notion of change of magnetic topology. (a) Two magnetic lines connecting different pairs of fluid elements, A $\leftrightarrow$ B and C $\leftrightarrow$ D, are twisted by the ideal fluid motion so as to locally generate strong spatial gradients, i.e. current densities. (b) Magnitude of the spatial gradients locally compensates the smallness of the non-ideal parameters in the generalised Ohm's law, by thus allowing the relaxation of the topological constraints that in the ideal limit, forbid the intersection of initially distinct magnetic lines; an $X$-point is formed. (c) Magnetic reconnection process undergone at the $X$-point, by which magnetic energy is dissipated and/or converted into other forms of plasma energy, which makes other magnetic configurations accessible to the plasma; the magnetic topology has globally changed as the magnetic lines now connect the pairs of fluid elements A $\leftrightarrow$ D and B $\leftrightarrow$ C.

Non-ideal effects allow a local violation of the topological conservations implied by the ideal Ohm's law (Cowling Reference Cowling1933; Alfvén Reference Alfvén1942; Batchelor Reference Batchelor1950; Elsasser Reference Elsasser1950a,Reference Elsasserb; Truesdell Reference Truesdell1950; Lundquist Reference Lundquist1951; Newcomb Reference Newcomb1958), which would otherwise prevent initially distinct magnetic lines embedded in the plasma to intersect. Mathematically speaking, this kind of conservation is a consequence of the fact that Faraday–Ohm's law in an ideal MHD plasma,

(2.6)\begin{equation} \frac{\partial {\boldsymbol B}}{\partial t}={\boldsymbol \nabla}\times({\boldsymbol U}\times{\boldsymbol B}), \end{equation}

is equivalent, thanks to a well-known vector identity (${\boldsymbol \nabla }\times ({\boldsymbol U}\times {\boldsymbol B})= {\boldsymbol U}({\boldsymbol \nabla }\boldsymbol {\cdot }{\boldsymbol B})- {\boldsymbol B}({\boldsymbol \nabla }\boldsymbol {\cdot }{\boldsymbol U})+ ({\boldsymbol B}\boldsymbol {\cdot }{\boldsymbol \nabla }){\boldsymbol U}-({\boldsymbol U}\boldsymbol {\cdot }{\boldsymbol \nabla }){\boldsymbol B}$) combined with the continuity equation (Truesdell Reference Truesdell1950), to a vector expression that corresponds to the null Lie-derivative of the vector field ${\boldsymbol B}/n$ dragged by the velocity field ${\boldsymbol U}$: this directly implies the Lagrangian invariance of magnetic lines from which a set of topological conservations follow, the most famous of which go under the names of ‘(Cowling–)Alfvén theorem’ (Cowling Reference Cowling1933; Alfvén Reference Alfvén1942), ‘Woltjer invariants’ Woltjer (Reference Woltjer1958) and ‘connection theorem’ (Newcomb Reference Newcomb1958) – see (Tur & Yanovsky Reference Tur and Yanovsky1993; Kuvshinov & Schep Reference Kuvshinov and Schep1997; Del Sarto et al. Reference Del Sarto, Califano and Pegoraro2006) and references therein for a more detailed discussion; see, e.g. Dubrovin, Nivikov & Fomenko (Reference Dubrovin, Nivikov and Fomenko1991, § III.23) for a definition of Lie derivative in tensor notation in the coordinate representation and Schouten (Reference Schouten1989, § II.8) for the Lie-derivative of tensor densities of arbitrary ‘weight’ (cf. also Lovelock & Rund (Reference Lovelock and Rund1989, p. 105) for the notion of a ‘relative tensor’). This is mathematically analogous to the set of topological conservations related to the fluid vorticity, which were well known to follow from the inviscid vorticity equation in a barotropic fluid (see Truesdell Reference Truesdell1954). In this context, the Alfvén theorem of magnetic flux conservation in an ideal MHD plasma is the mirror correspondent of the Helmholtz–Kelvin theorem of vorticity conservation in ideal hydrodynamics (Batchelor Reference Batchelor1950; Elsasser Reference Elsasser1950b; Truesdell Reference Truesdell1950; Axford Reference Axford1984; Greene Reference Greene1993). We incidentally note that some formal similarities can be also recognised in the eigenmode analysis of resistive tearing modes and of ideal instabilities in presence of viscosity in a Kolmogorov hydrodynamic flow (Fedele, Negulescu & Ottaviani Reference Fedele, Negulescu and Ottaviani2021). In terms of the RMHD equations above, the Lagrangian invariance of magnetic lines is expressed by the ideal limit ($d_e=\rho _s=S^{-1}=0$) of (2.1), although a finite $\rho _s$ alone allows preservation of the ideal MHD topological conservation, provided a redefinition of the stream function of the velocity field ${\boldsymbol U}$ according to $\varphi \to \varphi -\rho _s^2\nabla ^2\varphi$ (Pegoraro et al. Reference Pegoraro, Borgogno, Califano, Del Sarto, Echkina, Grasso, Liseikina and Porcelli2004). In particular, in the collisionless regime, even the nonlinear evolution of tearing-type instabilities can be shown to be ruled by the conservation of Lagrangian invariants whose existence is related to the condition $\rho _s\neq 0$ (Cafaro et al. Reference Cafaro, Grasso, Pegoraro, Porcelli and Saluzzi1998; Grasso et al. Reference Grasso, Califano, Pegoraro and Porcelli2001).

The early notion of magnetic reconnection has been formalised by Dungey (Reference Dungey1950, Reference Dungey1953) after the intuitions by Giovanelli and Hoyle. The former one first noted the occurrence of solar flares in correspondence with regions of local inversion of the magnetic field perpendicular to the Sun surface, and thus made the hypothesis that this could have been a signature of a mechanism of conversion from the magnetic energy to the plasma kinetic energy allowed by resistivity (Giovanelli Reference Giovanelli1946); the second one conjectured that the same may have occurred also in the terrestrial magnetotail (Hoyle Reference Hoyle1949). Resistivity is one of the non-ideal effects capable of violating alone the Lagrangian invariance of ${\boldsymbol B}$, so as the electron inertia, or the electron–electron viscosity, or a non-zero out-of-plane component of the rotational of the divergence of the pressure tensors, also are. All these non-ideal effects, together with other ones such as the Hall-term, can be synthetically expressed with the vector ${\boldsymbol \varPhi }_\varepsilon$ in generalised Ohm's law, $\varepsilon$ generally indicating the infinitesimally small parameter weighting the non-ideal contribution (we use the same normalisation of (2.1)–(2.2)):

(2.7)\begin{equation} {\boldsymbol E}+{\boldsymbol U}\times {\boldsymbol B}= {\boldsymbol \varPhi}_\varepsilon. \end{equation}

Magnetic reconnection at least requires ${\boldsymbol \nabla }\times {\boldsymbol \varPhi }_\varepsilon \neq 0$. Although ideal Ohm's law at the MHD scale takes the form of (2.7) with ${\boldsymbol \varPhi }_\varepsilon =0$, note that the dominant contribution to generalised Ohm's law, obtained by summing (A7)–(A8) multiplied by the respective charge over mass of the species, comes from electron momentum equation. That is, electrons fix the strongest ‘frozen-in’ condition between the plasma flow and the magnetic field, which, by neglecting all electron inertia, electron–electron viscosity and electron–ion viscosity (i.e. resistivity), and electron-pressure anisotropy, is expressed by the ‘ideal Hall-MHD Ohm's law’: ${\boldsymbol E}+{{\boldsymbol u}^e}\times {\boldsymbol B}= 0 \; \Leftrightarrow \; {\boldsymbol E}+{\boldsymbol u}^i\times {\boldsymbol B}= d_i ({\boldsymbol J}\times {\boldsymbol B})/{n}$. While at MHD scales, where ${\boldsymbol U}\simeq {\boldsymbol u}_i$, the frozen-in condition equally applies for both ions and electrons. At smaller spatial scales, the magnetic field can decouple from the ion fluid, but, as long as ${\boldsymbol E}+{{\boldsymbol u}^e}\times {\boldsymbol B}\simeq 0$ (or, better, as long as ${\boldsymbol E}+{{\boldsymbol u}^e}\times {\boldsymbol B}\simeq \nabla f$) holds, magnetic lines are dragged by the electron fluid flow. Note in this regard that (2.1) for the magnetic stream function $\psi$, corresponding to the $z$-component of Ohm's law, also expresses the variation of the $z$-component of the electron canonical momentum (cf. Appendix A).

The first analytical model of magnetic reconnection was the well-known Sweet–Parker model (Parker Reference Parker1957; Sweet Reference Sweet1958), which assumes a steady inflow condition in an $X$-point: in this case, the resistive reconnection steadily occurs in one point (the $X$-point) of a static current sheet of finite elongation $L$ and of large aspect ratio $L/a$, which asymptotically scales as $S$, when the Lundquist number is defined with respect to the current sheet thickness $a$, as in (2.1) – see figure 2(a) (when $S$ is instead defined – let us call it $S_L$ – with respect to the current sheet length, $L$, the scaling of the Sweet-Parker aspect ratio reads $L/a\sim S_L^{1/2}$). This steady reconnection scenario has been then extended to the inertia-driven regime by Wesson (Reference Wesson1990), for applications to the sawtooth crash in tokamaks, and by Bulanov, Pegoraro & Sakharov (Reference Bulanov, Pegoraro and Sakharov1992) for applications to reconnection in the electron-MHD regime. Variations to this model, in which a different choice of the boundary conditions around the reconnecting region has been made and different rates of magnetic reconnection have been obtained, have been done in subsequent works, starting from that of Petschek (Reference Petschek1964) – see also (Vasyliunas Reference Vasyliunas1975) and references therein.

Figure 2.

Comparison of four reconnection scenarios. (a) Sketch of the steady Sweet–Parker reconnection scenario. No particular hypotheses are done on the boundary conditions of the reconnecting current sheet, whereas steadiness and a specific scaling of the aspect ratio with non-ideal parameters like resistivity (i.e. $L/a\sim S^{-1}$, Parker Reference Parker1957; Sweet Reference Sweet1958) or electron inertia (Wesson Reference Wesson1990) are required. The configuration does not represent an instability: no distinction is made between equilibrium and perturbed quantities. (b) Typical ‘magnetic island’ configuration in the tearing-type reconnection. Periodic boundary conditions along the neutral line are required and hyperbolic patterns of the velocity field near the $X$- and $O$-points develop together with the island formation. The wavenumber of oscillations of the mode fixes the number of magnetic islands which are generated. In the linear stage, the island width is $w\simeq \delta \ll a$. (c,d) Merger configuration (c) of two undisplaced line currents (Syrovatskii Reference Syrovatskii1966a; Biskamp & Welter Reference Biskamp and Welter1980; Kleva et al. Reference Kleva, Drake and Waelbroeck1995), which resembles (d) a non-periodic version of the ‘coalescence instability’ scenario. The latter concerns a periodic structure constituted by magnetic islands in a chain (Finn & Kaw Reference Finn and Kaw1977; Pritchett, Lee & Drake Reference Pritchett, Lee and Drake1979; Bhattacharjee, Brunel & Tajima Reference Bhattacharjee, Brunel and Tajima1983). In the displaced current merger case, the attraction of two parallel current filaments squeezes the magnetic field that piles up and typically dissipates at the $X$-point: in this configuration, one may rather speak of ‘magnetic annihilation’ (cf. e.g. Priest & Sonnerup Reference Priest and Sonnerup1975; Watson & Craig Reference Watson and Craig1998; Gu et al. Reference Gu, Pegoraro, Sasarov, Golovin, Yogo, Korn and Bulanov2019). In the coalescence instability scenario, a somewhat similar process occurs (an out-of-plane current density being associated with the $O$-point of each island) once a non-rigid displacement along the neutral line perturbs the periodic configuration by making the islands get closer, pairwise. Similarly, a periodic 2-D ‘paving’ of magnetic cell-like structures of different shape (not shown here) can be destabilised the same way by planar, non-rigid displacements (Longcope & Strauss Reference Longcope and Strauss1993).

Instead, the notion of a tearing mode, the first example of a spontaneous reconnecting instability, was formalised a few years later in the pioneering work by Furth et al. (Reference Furth, Killeen and Rosenbluth1963). In this case, the current sheet is assumed to have periodic boundary conditions along the ‘resonant line’ where the wave vector of the linear perturbation is locally orthogonal to the equilibrium magnetic field. That is, if ${\boldsymbol B}_0(x_s)\boldsymbol {\cdot }{\boldsymbol k}=0$, then the equation $x=x_s$ defines a resonant line in two dimensions and a resonant surface in three dimensions. After a translation of the reference amplitude of the sheared magnetic field component, it can always be assumed that the sheared equilibrium magnetic field is zero at $x=x_s$, and changes sign in its neighbourhood: therefore, in slab geometry, the resonant line is usually termed the ‘neutral line’. It can be seen from the energy principle that the ${\boldsymbol B}_0(x_s)\boldsymbol {\cdot }{\boldsymbol k}=0$ condition minimises the stabilising role played by Alfvénic perturbations: this is why rational surfaces $(m,n)$ in large aspect ratio tokamaks are candidate resonant surfaces for tearing-type instabilities. We recall that a rational surface in a tokamak is a magnetic surface characterised by the condition $q(r)=-m/n$. In the large aspect ratio limit, i.e. the ‘cylindrical tokamak’ approximation, the safety factor reads $q(r)=B_\varphi r/(B_\theta R)$, where $B_\varphi$ and $B_\theta$ are respectively the toroidal and poloidal component of the magnetic field, and $r$ and $R$ are the minor and major radius of the tokamak, respectively. It is easy to verify that the condition ${\boldsymbol k}(r)\boldsymbol {\cdot }{\boldsymbol B}(r)=0$ defines a rational surface for ${\boldsymbol k}=m/(2{\rm \pi} r){\boldsymbol e}_\theta +n/(2{\rm \pi} R){\boldsymbol e}_\varphi$. The harmonic perturbation along the neutral line direction (the wave vector $k$ of (2.1)–(2.2)) induces a sinusoidal modulation of the perturbed magnetic field lines that gives rise to the characteristic ‘magnetic island’ pattern: in this case, a pair of $X$-points (at the local minima of $\psi _1$) delimits each magnetic island, and an elliptic point for the magnetic field, called the $O$-point (at the local maximum of $\psi _1$), is at the centre of the island – see figure 2(b) (see also White Reference White1983, Reference White1986) for a tutorial discussion about further features of tearing mode analysis that go beyond the purpose of the present work). In figure 2c,d), other reconnection scenarios related to the merger of current filaments (Syrovatskii Reference Syrovatskii1966a; Biskamp & Welter Reference Biskamp and Welter1980) and to the ‘coalescence’ of a chain of magnetic islands (Finn & Kaw Reference Finn and Kaw1977; Longcope & Strauss Reference Longcope and Strauss1993) are shown for comparison.

As tearing-type modes are linear instabilities, they grow exponentially in time, which marks an important, fundamental difference with respect to the Sweet–Parker scenario: the structure of the reconnecting current sheet is not static but is altered by the current density associated with the tearing mode perturbation, which makes the magnetic island grow in thickness until nonlinear saturation of the instability (see Ottaviani et al. Reference Ottaviani, Arcis, Escande, Grasso, Maget, Militello, Porcelli and Zwingmann2004 and references therein for a survey of different saturation scenarios). This dynamics is associated with hyperbolic patterns of the velocity field around both the $X$- and $O$-points, whereas only the hyperbolic flow at the $X$-point is required in the Sweet–Parker scenario. The hyperbolic flow at the $X$-point is responsible for the well-known quadrupole pattern of the stream function of the velocity field, and the elliptic and hyperbolic structure of $\psi$, $\varphi$ and of their derived fields near the critical points ($X$ and $O$) fix the symmetries of the eigenmode solutions. We incidentally note that this also provides a quite simple explanation (cf. Del Sarto et al. Reference Del Sarto, Pucci, Tenerani and Velli2016, App. A) of the quadrupole pattern associated with the out-of-plane magnetic field component at the $X$-point, in the regime where the Hall term is dominant in the generalised Ohm's law (quadrupole structure which is often recognised as an experimental proxy of Hall-mediated reconnection – cf. e.g. Deng & Matsumoto Reference Deng and Matsumoto2001): in this regime, the Lagrangian invariance of the magnetic field is expressed with respect to the electron flow (see Fruchtman Reference Fruchtman1991; Pegoraro et al. Reference Pegoraro, Borgogno, Califano, Del Sarto, Echkina, Grasso, Liseikina and Porcelli2004; Del Sarto et al. Reference Del Sarto, Califano and Pegoraro2006), which is in turn associated with the in-plane current density, and therefore to the spatial modulation of the out-of-plane magnetic field. If the equilibrium magnetic profile is even with respect to the neutral line, as it is typically assumed in analytical models for tearing modes (as we will do here), the matrix operators of (2.3) commute with the parity operator inducing the transformation $x\leftrightarrow -x$, which makes tearing modes have a fixed parity in $x$. Taking into account the parity in $y$ associated with the harmonic perturbation, the point symmetries of the $\psi _1$ and $\varphi _1$ fields near the critical points ($X$ and $O$) for tearing modes developing in an even magnetic equilibrium can be summarised as

(2.8)\begin{equation} {\text{X-point}: \left\{\begin{array}{@{}ll} \psi\ \text{even in} \ x, & \text{even in}\ y ,\\ \varphi\ \text{odd in}\ x, & \text{odd in}\ y, \end{array}\right. \quad \text{O-point}: \left\{\begin{array}{@{}ll} \psi\ \text{even in}\ x, & \text{even in}\ y ,\\ \varphi\ \text{odd in}\ x, & \text{odd in}\ y .\end{array}\right.} \end{equation}

These symmetries nonlinearly mix via the Poisson bracket operators, of which the terms of (2.1)–(2.2) are a rewriting, after linearisation (cf. (A1)–(A2), Appendix A). A local quadratic expansion of the eigenmodes near the critical points, which accounts for these nonlinear couplings, allows interesting insight in the early nonlinear, local dynamics of tearing-type instabilities (Pegoraro et al. Reference Pegoraro, Kuvshinov, Romanelli and Schep1995; Del Sarto et al. Reference Del Sarto, Marchetto, Pegoraro and Califano2011).

We conclude by recalling that several effects, in addition to the inclusion of further non-ideal terms in Ohm's law, have been considered, over the years, in the linear theory of tearing modes in reduced MHD. These include the effect of flows parallel to the neutral line, which can have a stabilising effect on the tearing mode instability (Bulanov, Sakai & Syrovatskii Reference Bulanov, Sakai and Syrovatskii1979; Syrovatskii Reference Syrovatskii1981; Faganello et al. Reference Faganello, Pegoraro, Califano and Marradi2010), and which can also make the latter compete with the Kelvin–Helmholtz instability (Hofmann Reference Hofmann1974; Chen & Morrison Reference Chen and Morrison1990; Ofman et al. Reference Ofman, Chen, Morrison and Steinolfson1991; Chen, Otto & Lee Reference Chen, Otto and Lee1997) – cf. also (Einaudi & Rubini Reference Einaudi and Rubini1986; Wang, Lee & Wei Reference Wang, Lee and Wei1988; Bettarini et al. Reference Bettarini, Landi, Rappazzo, Velli and Opher2006; Li & Ma Reference Li and Ma2012) for a numerical linear study; the stabilising effect of an in-plane magnetic field component orthogonal to the neutral line (see Nishikawa Reference Nishikawa1982; Somov & Verneta Reference Somov and Verneta1988, Reference Somov and Verneta1989); the role of equilibrium density gradients, which via the diamagnetic drift frequency gives a time-resonant character to the so-called drift-tearing modes (Coppi Reference Coppi1965; Drake & Lee Reference Drake and Lee1977; Coppi et al. Reference Coppi, Mark, Sugiyama and Bertin1979). These diamagnetic effects are particularly important in tokamaks, where they are related to the rotation of tearing structures and contribute to the nonlinear or resonant coupling of tearing modes with different mode numbers (see, e.g. Hicks, Carreras & Holmes Reference Hicks, Carreras and Holmes1984; Cowley & Hastie Reference Cowley and Hastie1988; Fitzpatrick et al. Reference Fitzpatrick, Hastie, Martin and Roach1993), but also influence their linear evolution and stability threshold (see, e.g. Mahajan et al. Reference Mahajan, Hazeltine, Strauss and Ross1978, Reference Mahajan, Hazeltine, Strauss and Ross1979; Migliuolo, Pegoraro & Porcelli Reference Migliuolo, Pegoraro and Porcelli1991; Grasso, Ottaviani & Porcelli Reference Grasso, Ottaviani and Porcelli2001; Yu Reference Yu2010). In tokamaks, also temperature gradients may induce, in the framework of a gyrokinetic description, the onset of the so-called ‘micro-tearing’ modes (see, e.g. Hazeltine & Ross Reference Hazeltine and Ross1975; Drake & Lee Reference Drake and Lee1977; Gladd et al. Reference Gladd, Drake, Chang and Liu1980; Connor, Cowley & Hastie Reference Connor, Cowley and Hastie1990), which cause the formation of microscopic magnetic islands that can affect turbulent transport (see Doerk et al. Reference Doerk, Jenko, Pueschel and Hatch2011) and foster stochastic fluctuations of the magnetic field, which can in turn affect tearing modes (Carreras, Rosenbluth & Hicks Reference Carreras, Rosenbluth and Hicks1981) and the magnetic confinement (see, e.g. Firpo Reference Firpo2015). Also, we recall that the linear tearing mode theory has been extended and adapted so as to study modes simultaneously and interdependently occurring on two (or more) sufficiently close resonant surfaces, i.e. the double- (or multiple-) tearing modes (Furth, Rutherford & Selberg Reference Furth, Rutherford and Selberg1973; Pritchett, Lee & Drake Reference Pritchett, Lee and Drake1980; Wang et al. Reference Wang, Wei, Wang, Zheng and Liu2011); or to include the so-called ‘neo-classical effects’ related to pressure gradient and toroidal curvature in tokamak plasmas. Although the latter effects are intrinsically nonlinear, they can be self-consistently included in the study of the linear growth of the so-called ‘neo-classical tearing’ modes (NTMs, in specialised ‘jargon’): these are modes developing out of a seed magnetic island (e.g. produced by previous tearing-type instabilities which have saturated and have become stable) and which are driven by the so-called ‘boot-strap current’ (see, e.g. Hahm Reference Hahm1988; Lütjens & Luciani Reference Lütjens and Luciani2002; Wilson Reference Wilson2012). In addition to the large number of works developed over the years on these subjects, attention has also been recently drawn to the role that a background plasma turbulence can have in driving the growth of NTMs (see, e.g. Muraglia et al. Reference Muraglia, Agullo, Benkadda, Garbet, Beyer and Sen2009, Reference Muraglia, Agullo, Benkadda, Yagi, Garbet and Sen2011; Agullo et al. Reference Agullo, Muraglia, Benkadda, Poyé, Dubuit, Garbet and Sen2017a,Reference Agullo, Muraglia, Benkadda, Poyé, Dubuit, Garbet and Senb; Choi Reference Choi2021). These latter subjects are somewhat related (see, e.g. Brennan et al. Reference Brennan, Strait, Turnbull, Chu, La Haye, Luce, Taylor, Kruger and Pletzer2002) to the further topic of forced magnetic reconnection, which has also been considered in the framework of tearing mode theory, by looking at the way an external forcing affects the stability and growth of the linear modes (see, e.g. Hahm & Kulsrud Reference Hahm and Kulsrud1985; Wang & Bhattacharjee Reference Wang and Bhattacharjee1997; Fitzpatrick Reference Fitzpatrick2008). None of these subjects is however of further concern to us; here, in the following, we will focus on the linear problem related to (2.1)–(2.2) only, with respect to which of each of these ingredients would provide further elements of analytical complication that would lead us beyond the purpose of this work.

2.2 Eigenvalue of the linear problem and reconnection rate in tearing-type modes

In the literature, the notion of ‘reconnection rate’, let us name it as $\mathcal {R}_{rc}$, formally refers to the rate at which the magnetic flux changes over time close to the $X$-point during a magnetic reconnection event associated with it. In general, we can write

(2.9)\begin{equation} \mathcal{R}_{rc}\equiv \frac{1}{\varPhi_\varOmega({\boldsymbol B})}\frac{{\rm d} \varPhi_\varOmega({\boldsymbol B})}{{\rm d} t}, \quad \varPhi_\varOmega({\boldsymbol B})\equiv \int_\varOmega {\boldsymbol B}\boldsymbol{\cdot} {\rm d}{\boldsymbol S}, \end{equation}

where $\varOmega$ is the resonant surface, which, in a planar configuration, corresponds to the surface generated by translation of the neutral line along the invariant coordinate ($z$ in the slab coordinate system we have chosen). If we assume the planar current sheet (or the neutral line) to be static and along the $y$ direction, then $\varOmega =\ell _y\ell _z$, where the extension of the length $\ell _y$ depends on the specific reconnection process considered.

The way this rate is evaluated clearly depends on the reconnection scenario which is considered. A few words on this topic are therefore useful, since the contrasting terminology acquired on this subject by specialised literature dealing with different reconnection scenarios (e.g. steady, Sweet–Parker-like versus tearing-type reconnection) can be sometimes confusing (cf. also Biskamp (Reference Biskamp2000, p. 54) for a brief discussion as well as Appendix B). This is especially important, since the reconnection rate is related to the rate at which magnetic energy is converted into other forms of plasma energy and therefore to the power released in the form of thermal energy and particle acceleration during a reconnection event. This subject is at the basis of many open questions turning around magnetic reconnection in both laboratory and astrophysical plasmas. These subjects, still debated, concern for example: the model capable of accounting for the short time scale of the sawtooth crash and of other reconnection-related disruptions in tokamaks (see, e.g. Wesson Reference Wesson1986, Reference Wesson2004; Porcelli, Boucher & Rosenbluth Reference Porcelli, Boucher and Rosenbluth1996; Boozer Reference Boozer2012; Jardin, Krebs & Ferraro Reference Jardin, Krebs and Ferraro2020); a theoretical model of the fast release of energy during solar flares and coronal mass ejections (see, e.g. Shibata Reference Shibata1998; Cassak & Shay Reference Cassak and Shay2012; Aschwanden, Xu & Jing Reference Aschwanden, Xu and Jing2014; Aschwanden et al. Reference Aschwanden, Holman, O'Flannagain, Caspi, McTiernan and Kontar2016; Janvier Reference Janvier2017; Aschwanden et al. Reference Aschwanden, Caspi, Cohen, Holman, Jing, Kretzschmar, Kontar, McTiernan, Mewaldt and O'Flannagain2019; Aschwanden Reference Aschwanden2020); the mechanism of energy release at the basis of $X$-rays and $\gamma$-ray emissions (gamma-ray bursts) in pulsars, magnetars and nebulae (see, e.g. Lyutikov Reference Lyutikov2003; Tavani et al. Reference Tavani, Bulgarelli, Vittorini, Pellizzoni, Striani, Caraveo, Weisskopf, Tennant, Pucella and Trois2011; Uzdensky Reference Uzdensky2011); the processes by which the kinetic and magnetic energy of the photospheric plasma are likely to heat up the expanding solar wind by a factor $\sim 10^6$ over the distance of just a solar radius, supposedly via reconnection in the turbulent coronal plasma (see, e.g. Leamon et al. Reference Leamon, Matthaeus, Smith, Zank, Mullan and Oughton2000; Matthaeus & Velli Reference Matthaeus and Velli2011); the energy injected in the magnetosphere via solar wind–magnetosphere interactions and released during geomagnetic storms (see, e.g. Lakhina & Tsurutani Reference Lakhina and Tsurutani2016), and the impact that the ensuing space–weather perturbations may have on terrestrial biosphere and human activity (see National Research Council (2008) and further more specific studies like, e.g. Gopalswamy Reference Gopalswamy2016; Nelson Reference Nelson2016; Eastwood et al. Reference Eastwood, Hapgood, Biffis, Benedetti, Bisi, Green, Bentley and Burnett2018; Knipp et al. Reference Knipp, Fraser, Shea and Smart2018). Regardless of whether the tearing mode theory be actually relevant and capable to explain these phenomena or not, it should be noted that the comparison of theoretical estimates of the reconnection rate with the characteristic time scales directly measured or inferred from experiments or simulations is probably the main criterion presently used to assess the pertinence of tearing-type instabilities or of alternative magnetic reconnection models for these energy releasing processes.

The operational definition by which the reconnected flux can be computed according to (2.9) has been refined over time in different contexts, which also keep account of observational difficulties related to issues encountered in both direct experimental and numerical measurements (see Parker Reference Parker1957; Petschek Reference Petschek1964; Syrovatskii Reference Syrovatskii1966a,Reference Syrovatskiib; Bratenahl & Yeates Reference Bratenahl and Yeates1970; Baum, Bratenahl & White Reference Baum, Bratenahl and White1971; Parker Reference Parker1973; Schnack Reference Schnack1978; Park, Monticello & White Reference Park, Monticello and White1984; Chen et al. Reference Chen, Otto and Lee1997; Hesse, Forbes & Birn Reference Hesse, Forbes and Birn2005; Comisso & Bhattacharjee Reference Comisso and Bhattacharjee2016; Grasso et al. Reference Grasso, Borgogno, Tassi and Perona2020). For the purpose of the present discussion, we here rely on the formal definition: combination of (2.9) with (2.7) via Faraday's law and Kelvin–Stokes theorem of circulation allows us to write

(2.10)\begin{equation} \mathcal{R}_{rc} = \frac{1}{\varPhi_S({\boldsymbol B})} \int_\varOmega \left(\frac{\partial {\boldsymbol B}}{\partial t}-{\boldsymbol \nabla}\times({\boldsymbol U}\times{\boldsymbol B})\right)\boldsymbol{\cdot} {\rm d}{\boldsymbol S} ={-}\frac{c}{\varPhi_S({\boldsymbol B})} \int_\varOmega ({\boldsymbol \nabla}\times\mathcal{\boldsymbol F}_\varepsilon )\boldsymbol{\cdot} {\rm d}{\boldsymbol S}, \end{equation}

where the curl term in the first integral generally accounts for the evolution of the surface over which the magnetic flux is calculated, as it is dragged by the bulk plasma velocity ${\boldsymbol U}$. If the surface over which the flux is evaluated is static, like it can be chosen in the Sweet–Parker scenario or for non-propagating tearing-like modes, that contribution is null. This is the case in which we are interested, where we consider a non-evolving integration surface along the neutral line, whose extension in the $y$ direction depends on the reconnection scenario considered, as it will be specified below. If, under this assumption, we specify the surviving terms of the first integral of (2.10) for the slab geometry, RMHD variables ${\boldsymbol B}=\boldsymbol {\nabla }\psi \times \boldsymbol {e}_z+b\boldsymbol {e}_z$ of (2.1)–(2.2), we can write ${\rm d}{S}={{\rm d} y}\,{\rm d}z\,{\boldsymbol e}_x$ and therefore

(2.11)\begin{equation} \varPhi_\varOmega({\boldsymbol B})=\int_\varOmega {\boldsymbol B}\boldsymbol{\cdot} {\rm d}{\boldsymbol S} = \int_\varOmega {B}_{x,1} \,{{\rm d} y} \,{\rm d}z =\int_{l_z}\,{\rm d}z \int_{l_y} \frac{\partial \psi_1}{\partial y} \,{{\rm d} y}. \end{equation}

Thus, posing the $X$-point to lay on the line $x=0$, we obtain

(2.12)\begin{equation} \mathcal{R}_{rc} = \frac{\int_{l_y}{{\rm d} y}\;\partial^2 \psi_1(0,y,t)/\partial t \partial y}{\int_{l_y}{{\rm d}y}\;\partial \psi_1(0,y,t)/\partial y}. \end{equation}

In the case of tearing modes, it is easy to see from the symmetries of the eigenmode solutions in (2.8) that it is appropriate to consider the integration interval in $y$ to not have the $X$-point in the middle (otherwise the spatial integral would vanish). For a single tearing mode, the interval can be therefore taken to extend from the $X$-point to the $O$-point (although each $O$-point is delimited by two $X$-points, the periodicity of the configuration makes the numbering of $X$- and $O$-points be in a 1:1 correspondence). By construction, the eigenmode of the magnetic stream function has a separable dependence on time and on space of the form $\psi _1(x,y,t)=f(x,y) {\rm e}^{\gamma t}$. In reality, although we are here anticipating the formal results, as we will discuss in the next sections, it is easy to determine from the magnetic island shape that the solution is of the form $\psi _1(x,y,t)=g(x)\cos (ky){\rm e}^{\gamma t}$, if $y=0$ is taken at the $O$-point – cf. figure 2(b). Hence, one immediately obtains (the label $TM$ standing for ‘tearing mode’)

(2.13)\begin{equation} \mathcal{R}_{TM}(k) = \gamma(k), \end{equation}

that is, for a tearing-type instability of wavenumber $k$, the notions of reconnection rate and of growth rate are equivalent. The quadratic dependence of the magnetic energy on $\psi$ (cf. (A5)) trivially implies that for a single tearing mode, the rate at which magnetic energy $\mathcal {E}_{B}=\int |{\boldsymbol \nabla }\psi _1|^2\,{{\rm d} x}\,{{\rm d} y}$ is converted into other forms of energy during the linear stage of the instability is $2 \mathcal {R}_{TM} = 2\gamma$.

Of course, if several tearing modes are simultaneously unstable, and/or if one or more reconnection instabilities are in their nonlinear stage, the identification $\mathcal {R}_{rc}\leftrightarrow \gamma (k)$ is not correct any longer, even if a dominant mode can be assumed to dominate the reconnection process. In these cases, an evaluation of the global reconnection rate in a given volume, more accurate at least from a numerical point of view, can be made by taking

(2.14)\begin{equation} \mathcal{R}_{rc} \simeq \frac{1}{2 \mathcal{E}_c} \left|\frac{\partial\mathcal{E}_c}{\partial t}\right|, \end{equation}

where $\mathcal {E}_c$ is a component of the total energy (cf. (A5)) involved in the reconnecting process and the factor $1/2$ accounts for the quadratic relation between the fields and the energy components. Usually, $\mathcal {E}_c$ is the magnetic energy component $\mathcal {E}_{B}$, which in most reconnection regimes (see, e.g. White Reference White1983, Reference White1986) can be shown to provide the ‘reservoir’ of energy of the instability, which is ‘dissipated’ (or, more generally, converted into other forms) during the process. It should be however noted that in some regimes of tearing-like instabilities, numerically studied by using ‘large’ values of $d_e/L_0$ (i.e. $d_e/L_0\lesssim 1$ but not asymptotically smaller than unity), the role of $\mathcal {E}_{B}$ has been observed to be played, instead, by the electron kinetic energy $\mathcal {E}_{J}=\int \,{\rm d}_e^2|{\nabla }^2\psi _1|^2\,{{\rm d} x}\,{{\rm d}\kern 0.05em y}$ associated with the equilibrium current sheet (Del Sarto, Califano & Pegoraro Reference Del Sarto, Califano and Pegoraro2005; Del Sarto et al. Reference Del Sarto, Califano and Pegoraro2006). Also, the increase of some kinetic energy component could be used for the estimate of (2.14), if it assumed that the electromagnetic energy released during the reconnection event is dominantly transferred to particle heating and/or acceleration, or to radiation. This is the way by which, for example, the reconnection rate of processes at the basis of solar flares and coronal mass ejections, or of the sawtooth cycle in tokamaks, is indirectly inferred from observational measures.

Some alternative definitions of ‘reconnection rate’, provided in terms of local estimates of quantities, are also frequently used in the literature. These further definitions, briefly outlined and discussed in Appendix B for comparison, are not of concern for the purpose of this work, since they mostly deal with stationary reconnection or provide approximations to the more accurate (2.13), valid in the case of interest to us.

We conclude this section with a note on the notion of ‘fast’, often used in the literature for models devised in the attempt of modelling the rapid reconnection processes observed in Nature and in experiments, which display growth rates faster – in an asymptotic sense – than those predicted by the early reconnection models (i.e. resistive Sweet–Parker and resistive tearing modes). On the one hand, the notion of ‘fast’ has changed over the years, both depending on the specific reconnection context (space or laboratory) and on the physical process allowing a relative increase of the reconnection rate (e.g. electron inertia, the inclusion of the Hall-effect in Ohm's law, the accounting for other nonlinear effects, 3-D effects and secondary instabilities – see appendices of Del Sarto & Ottaviani Reference Del Sarto and Ottaviani2017) for a short historical review on these subjects). On the other hand, it is unlikely that a primary spontaneous tearing mode developing on a static current sheet overtake values of order unity, when measured with respect to the normalisation time $\tau _A$ defined as below ((2.4)–(2.5)) – see Pucci & Velli (Reference Pucci and Velli2014), Del Sarto et al. (Reference Del Sarto, Ottaviani, Pucci, Tenerani and Velli2018), and § IX of Betar et al. (Reference Betar, Del Sarto, Ottaviani and Ghizzo2020) for more detailed discussions. It is however worth noting that for linear tearing modes, the normalised growth rate evaluated from (2.1)–(2.2) typically becomes of the order of unity or, better, of some decimal fraction of it, when the shear length becomes comparable to the microscopic non-ideal scales of interest: that is, $\gamma \tau _A\sim O(10^{-1})-O(1)$, when some among $d_e$, $\rho _s$ or $(S^{-1}/\gamma )^{1/2}$ become of the order of $O(a/10)-O(a)$. However, this also generally means that we are out of the limits of applicability of the boundary layer theory (cf. § 3, next, and note that we leave aside, here, the open problem of discussing the validity of an extended fluid modelling at these kinetic scales): in this case, the analytical estimates of $\gamma (k)$, which we are going to develop below, are not applicable and a numerical computation is instead required. See figure 3 and table 2 of Del Sarto et al. (Reference Del Sarto, Marchetto, Pegoraro and Califano2011) for an example in the collisionless limit.

Figure 3.

Growth rates of collisionless tearing modes numerically computed from (2.1)–(2.2) with the solver of Betar et al. (Reference Betar, Del Sarto, Ottaviani and Ghizzo2020), for values of the non-ideal parameters at the limits of (or beyond the) applicability of boundary layer theory, and for a wavelength $ka\sim O(1)$ (likely value for modes on large aspect ratio current sheets having a microscopic thickness – a condition typically met, e.g. in MHD turbulence (see, e.g. Franci et al. Reference Franci, Landi, Matteini, Verdini and Hellinger2016; Del Sarto & Pegoraro Reference Del Sarto and Pegoraro2017; Franci et al. Reference Franci, Papini, Del Sarto, Hellinger, Burgess, Matteini, Landi and Montagud-Camps2022). Note that once $\gamma$ approaches unity, the effect of the non-ideal parameters on the magnetic equilibrium $\psi _0$ is not negligible anymore, and (2.1) and (2.2) in the resistive limit (not shown here) should be modified accordingly (cf. Appendix A). (a) Case of cold electrons. (b) Case of warm electrons, for different values of $d_e/a$ and for values of $\rho _s/a$ approaching one.

3.1 Orderings and instability parameter $\varDelta '$

In the present section, we discuss the different orderings of the operators in the eigenvalue problem, which will allow us to obtain the governing equations in different regions of the domain, and we introduce the notion of an ‘instability parameter’.

In the outer, ideal MHD region, the ordering of different operators in (2.3) is

(3.1)\begin{equation} \gamma\sim 0, \quad \mathcal{L}\sim 1, \quad \mathcal{F}\simeq 1, \quad \mathcal{A}\sim 1, \quad \mathcal{B}=\mathcal{C}={\rm i}k\psi_0'\sim 1. \end{equation}

Solving the equations resulting from this ordering allows us to obtain the solution in the outer region, $\psi _{\text {out}}$.

At the neutral line, the equilibrium magnetic field (shear field) vanishes and changes its sign when reconnection occurs. Therefore, this field should be an odd function in the vicinity of $x=0$. Hence, $\psi _0$, whose $x$-derivative represents the shear field $B_y^0(x)$, is an even function. Moreover, close to the neutral line, $\psi _0(x)$ is continuous at least up to its first derivative. We then expand it as a Taylor series which, for symmetric tearing modes in slab geometry, has non-zero coefficients only for even powers of $x$,

(3.2)\begin{equation} \psi_0(x)|_{x\to 0}\simeq C_0+C_2x^2+C_4x^4+\cdots. \end{equation}

Therefore, we can relate $\delta$ to the intensity of the equilibrium magnetic field and to the gradient of the associated current density inside of the inner layer, when looked at from the matching (or from the outer) region. So we write

(3.3)\begin{equation} \left.\psi_0'\right|_{x\to 0}\sim \left.\psi_0'''\right|_{x\to 0}\sim x|_{\text{inner}}\simeq {\delta}. \end{equation}

Before going further, we note that condition (3.2) is too restrictive, in general, for slab tearing modes on magnetic equilibria that are not symmetric with respect to the shear coordinate. An example is provided by tearing modes in a cylinder, for which $C_1=0$ but $C_{3} \neq 0$ (Bertin Reference Bertin1982; Militello et al. Reference Militello, Huysmans, Ottaviani and Porcelli2004, Reference Militello, Borgogno, Grasso, Marchetto and Ottaviani2011). It must be also noted that the inner layer width, $\delta$, has not yet been ‘defined’ here. As a quantity, it only explicitly appears in boundary layer calculations as a normalisation scale for the differential equations in the non-ideal region, whose boundaries are defined by the condition $|x|/ \delta \sim O(1)$. To date, indeed, no general criterion exists to quantitatively define $\delta$, whose estimation is made, when possible, thanks to further hypothesis and heuristic ansatz (e.g. comparison with the characteristic scales of the problem). In the present section, we therefore assume $\delta$ to be simply defined as the innermost layer width (in § 8, we will propose a quantitative general definition of $\delta$, whose appropriateness in the different regimes we will prove numerically by comparing it with theoretical estimates).

Unless a large aspect ratio current sheet is considered, for which an almost continuous spectrum of wavenumbers can be destabilised that also allows $k$ to be large, the further ordering $\partial _x\gg k$ is assumed for $x\ll 1$. This implies the following orderings inside the non-ideal region:

(3.4)\begin{align} \mathcal{L}\simeq \partial_x^2,\quad \mathcal{F}= 1-d_e^2\mathcal{L}, \quad \mathcal{A}={-} k\delta(1-\mathcal{L})\sim{-}k\delta\mathcal{L}, \quad \mathcal{B}= ik\delta\left[(1-d_e^2)-\rho_s^2\mathcal{L}\right]. \end{align}

The corresponding equations are solved for the ‘inner’ functions $\psi _1=\psi _{\text {in}}$ and $\varphi _{1}=\varphi _{\text {in}}$.

The matching with the outer solution is assumed as a boundary condition to be imposed across an intermediate matching layer, where one must compare the asymptotic series $\sum a_{n}^{(\text {in})}(x/ \delta )^n$ for $(x/ \delta )\gg 1$ and $\sum _n a_{n}^{(\text {out})}x^n$ for $x\ll 1$ (we recall that all lengths appear here as normalised to $L_0=a$), respectively representing the inner and outer solutions – see Bender & Orszag (Reference Bender and Orszag1978, § 9). This translates in the condition

(3.5)\begin{equation} \lim_{x/ \delta\to\pm\infty} a_{m}^{(\text{in})}\,(x/ \delta)^m\simeq \lim_{x\to \pm 0} a_{m}^{(\text{out})}\,x^m\quad \begin{array}{c} \text{(for }m\text{ that corresponds to}\\ \text{the leading term of the series).} \end{array} \end{equation}

Condition (3.5) is usually expressed in a looser notation as

(3.6)\begin{equation} \lim_{x/ \delta\to\pm\infty}\psi_{\text{in}}(x/ \delta)=\lim_{x\to \pm 0}\psi_{\text{out}}(x). \end{equation}

Although (3.6) formally compares the two numerical values of the limits of the eigenfunctions solving the differential equations, in the following, we will mean this expression as a shortcut writing of (3.5), as it is usually done in tearing mode analysis.

The dependence on the outer solution becomes then explicit through the relation

(3.7)\begin{equation} \int_{-\infty}^{+\infty}\psi_{\text{in}}''\,{\rm d}\left({x}/ {\delta}\right)= \psi_{\text{out}}(0)\varDelta', \end{equation}

which introduces the instability parameter (Furth et al. Reference Furth, Killeen and Rosenbluth1963) related to the discontinuity which is met in the first derivative of $\psi _{\text {out}}$ as $x\to \pm 0$,

(3.8)\begin{equation} \varDelta'(k;\psi_0)\equiv \lim_{\epsilon\to 0} \frac{\psi_{\text{out}}'(\epsilon)-\psi_{\text{out}}'(-\epsilon)}{\psi_{\text{out}}(0)}. \end{equation}

In the expression above, $\epsilon > 0$. Definition (3.8) enters in (3.6) as

(3.9)\begin{equation} \int_{-\infty}^{+\infty} \psi''_{\text{in}}\,{\rm d}(x/ \delta) = \psi'_{\text{in}}(+\infty) - \psi'_{\text{in}}(-\infty) = \psi'_{\text{out}}({+}0) - \psi'_{\text{out}}(- 0) = \psi_{\text{out}}(0)\varDelta', \end{equation}

where we have used the key fact that the outer solution $\psi _{\text {out}}$ is continuous as $x/L_0\to 0$ while this is generally not the case for its derivative $\psi _{\text {out}}'$ (as it is found a posteriori when solving the inner equation for tearing-type modes).

The linear problem in the inner region is closed by combining (3.7) with the relation that can be established between $\psi _{\text {in}}$ and $\varphi _{\text {in}}$. According to the inner region ordering, (2.2) becomes $\gamma \varphi ''_{\text {in}} = {\rm i}kx\psi ''_{\text {in}}$. Therefore,

(3.10)\begin{equation} \int_{-\infty}^{+\infty} \psi''_{\text{in}}\,{\rm d}\left({x}/ {\delta}\right) ={-} ({\rm i} \gamma/ k) \int_{-\infty}^{+\infty} \frac{\varphi''_{\text{in}}}{ x/ \delta} \,{\rm d}\left({x}/ {\delta}\right). \end{equation}

It must be however noted that when more than one non-ideal parameter is involved and/or when the microscopic scale of variation of $\psi$ and $\varphi$ differ (this happens, for example, in 2-D electron-MHD (Bulanov et al. Reference Bulanov, Pegoraro and Sakharov1992), where the fluid stream function is related to the fluctuation $b$ of the $B_z$ magnetic component and an equation for $b$ replaces that for $\varphi$), the matching procedure above requires more care, as we are going to see in § 7. In these cases, more than two boundary layers must be considered and the corresponding solutions must be matched – see Pegoraro & Schep (Reference Pegoraro and Schep1986), Porcelli (Reference Porcelli1991) and Bulanov et al. (Reference Bulanov, Pegoraro and Sakharov1992). These are the cases for which heuristic estimations are difficult, since a different width must be associated with each boundary layer (e.g. $\delta _2$ and $\delta _1$, with $\delta _1<\delta _2$). In this kind of analysis, the width of each layer appears as the characteristic normalisation scale, say $l_{\text {norm}}$, with respect to which to consider the limits $x/l_{\text {norm}}\ll 1$ and $x/l_{\text {norm}}\gg 1$ while performing the asymptotic matching in the intermediate region between the layers. In this sense, the limit $\varepsilon \to 0$ of (3.8) also must be reinterpreted, which we should read as $\varepsilon =|x|/l_{\text {norm}}$ for $l_{\text {norm}}=L_0=a$. At the same time, also the relations between $\varDelta '$ and the spatial gradients of $\varphi _{\text {in}}$ expressed by combinations of (3.7) and (3.10) become non-trivial, as we will discuss in § 10.

3.2 Small- and large-$\varDelta '$ limits in slab geometry and in tokamaks

For each eigenvalue problem in slab geometry, the numerical value of $\varDelta '$ depends both on the choice of the magnetic equilibrium profile and on the value of $k$. A classification of regimes of slab tearing modes can be generally done depending on the numerical comparison of $(\varDelta ')^{-1}$ and of the innermost layer width, $\delta _1$, regardless of the number of boundary layers involved. For simplicity of notation, we generically use for them the symbol $\delta$, although it is only later (§ 8) that we will provide some argument to identify the ‘reconnecting layer width’. In particular, one speaks of a small-$\varDelta '$ limit for $\varDelta '\delta \ll 1$ and of a large-$\varDelta '$ limit for $\varDelta '\delta \gg 1$.

Once the scaling of $\delta$ on the non-ideal parameters involved is known, the large- and small-$\varDelta '$ limits in slab geometry can be made to respectively correspond to the small and large wavelength limits of the tearing dispersion relation, which for each reconnection regime can be defined with respect to the comparison of $kL_0=ka$ with powers of the non-ideal parameters involved (as it has been discussed for example by Bulanov (Reference Bulanov2017) for the purely resistive case).

In particular, since $\varDelta '$ defined through (3.7)–(3.8) is, when analytically obtained (cf. example below – § 4.1), a continuous function of the variable $k$, the large wavelength limit can be postulated to correspond to a power-law dependence of the kind

(3.11)\begin{equation} \lim_{ka\to 0} \varDelta'(ka)=(ka)^{{-}p},\quad p>0, \end{equation}

where $p$ depends on the initial equilibrium profile (Del Sarto et al. Reference Del Sarto, Pucci, Tenerani and Velli2016; Pucci et al. Reference Pucci, Velli, Tenerani and Del Sarto2018; Betar et al. Reference Betar, Del Sarto, Ottaviani and Ghizzo2020). While the limit above corresponds to $\varDelta '\to \infty$, the marginal stability condition $a\varDelta '(ka)\to 1$ is approached as $ka\to 1$.

A fundamental difference between large- and small-$\varDelta '$ tearing-type modes in tokamak devices and their corresponding small- and large-$k$ limits in slab geometry must be however remembered, when the results of the two models are compared: in tokamak devices, the wave vector is fixed by the resonant surface on which reconnection occurs, and the large- or small-$\varDelta '$ condition is therefore determined by the specific shape of the unstable magnetic profile (sometimes in turn determined by some ideal instability, which has previously occurred and which has modified the otherwise stable magnetic configuration). In slab Cartesian geometry, instead, the transition between the small- and large-$\varDelta '$ limit occurs by moving along the $k$ interval for a fixed, tearing-unstable magnetic profile, and depends on the aspect ratio of the associated current sheet (see figure 5).

Figure 5.

Scheme of the correspondence between linear reconnecting instabilities in the RMHD slab geometry limit and in the RMHD ‘cylindrical tokamak’ approximation, which is compatible with the strong guide field assumption (Strauss Reference Strauss1976). At the bottom of the figure, we highlight the fact that the slab geometry easily allows modelling both of tearing modes of different wavenumber (centre frame) and of the double tearing mode, i.e. the tearing-like instability simultaneously occurring on two sufficiently close resonant surfaces (Furth et al. Reference Furth, Rutherford and Selberg1973), of which we have shown here only the symmetric case (rightmost frame). Note however that we have considered here the ‘cylindrical tokamak’, or ‘large aspect ratio tokamak’ approximation, in the strict limit of $\partial / \partial \varphi \to 0$, which maps into the $\partial / \partial z \to 0$ assumption that we consider in this article. In reality, the RMHD modelling also allows for inclusion of the neglected derivatives with respect to the axis-symmetric coordinate ordered as $\partial / \partial \varphi \sim \partial / \partial z \sim \varepsilon _B$, in terms of the ratio $\varepsilon _B$ between the in-plane and the guide magnetic field components at equilibrium (cf. (A11) and discussion in the last paragraph of Appendix A).

We remark indeed that the distinction between small- and large-$\varDelta '$ limits has been historically introduced to characterise, in terms of the instability parameter $\varDelta '$ defined by Furth et al. (Reference Furth, Killeen and Rosenbluth1963), two different types of unstable modes observed in tokamaks. In tokamak physics, for positive values of $\varDelta '$, a distinction is made between the tearing mode or constant-$\psi$ mode, first studied in the cylindrical geometry approximation by Furth et al. (Reference Furth, Killeen and Rosenbluth1963) (and which formally corresponds to $\varDelta '\delta \ll 1$) and the internal kink mode or $m=n=1$ non-ideal kink mode, first identified in cylindrical geometry by Coppi et al. (Reference Coppi, Galvao, Pellat, Rosenbluth and Rutherford1976) (and which formally has $\varDelta '=\infty$) (we recall indeed that, even if $\varDelta '<0$, ideal instabilities (i.e. the ‘ideal kink mode’) can also occur in a tokamak, depending on the value of the safety factor $q$ – cf. Wesson (Reference Wesson1990), § 6). In a cylindrical tokamak approximation, these two modes generally occur on different magnetic surfaces and for specific values of the wavenumbers, and also display different relations between the eigenfunctions $\varphi _1$ and $\psi _1$: differently from the tearing mode, for which the fluid displacement $\xi \equiv k^2 \varphi _1/\gamma$ is proportional to $\sim \psi _1/x$, the internal kink mode corresponds to a rigid displacement $\xi \sim const.$ of the plasma inside the inner region; more specifically, for the internal kink mode, $\xi \sim const.$ in the whole subdomain $r\leq r_s$, where $r_s$ is the radius of the resonant surface, whereas $\xi = 0$ for $r>r_s$ (see, e.g. Porcelli Reference Porcelli1987; Del Sarto & Ottaviani Reference Del Sarto and Ottaviani2017 for a more detailed discussion). It has been however shown by Ara et al. (Reference Ara, Basu, Coppi, Laval, Rosenbluth and Waddell1978) that, thanks to the formal analogy between the two corresponding eigenvalue problems in cylindrical and Cartesian geometry, the transition between the two types of modes can be modelled in a slab geometry configuration by varying the value of the wavenumber for a fixed magnetic equilibrium and for fixed values of the non-ideal parameters. This corresponds to a transition between the small- and the large-$\varDelta '$ limits. In this context, both the tearing (constant-$\psi$) mode and the internal-kink mode of tokamak physics can be considered as examples of tearing-type modes when modelled in slab, Cartesian geometry. In this case, a larger ‘free energy’ can be associated with tearing-type modes in the large-$\varDelta '$ limit. This is the approach we take here.

Another fundamental difference between tearing and internal kink modes in tokamaks on the one side, and tearing-type modes in slab geometry on the other side, is the fact that in slab geometry, it makes sense to identify a further wavenumber ‘regime’ that is characterised by the condition $\varDelta '\delta \sim 1$ (see Loureiro, Schekochihin & Cowley Reference Loureiro, Schekochihin and Cowley2007; Bhattacharjee et al. Reference Bhattacharjee, Huang, Yang and Rogers2009; Del Sarto et al. Reference Del Sarto, Pucci, Tenerani and Velli2016; Betar et al. Reference Betar, Del Sarto, Ottaviani and Ghizzo2020) at the varying of the non-ideal parameters involved. This wavelength limit, which is met thanks to the possibility to perform a ‘continuous’ transition from $\varDelta '\delta \ll 1$ to $\varDelta '\delta \gg 1$, corresponds to the fastest growing mode that can be destabilised in a periodic current sheet when a continuum spectrum of wavenumbers is admitted (see Furth et al. Reference Furth, Killeen and Rosenbluth1963, appendix D, and also Biskamp Reference Biskamp1982). Its scaling is exemplified in figure 6 for the case of purely resistive tearing modes (i.e. $d_e=0,\rho _s=\nu =0$ in (2.1)–(2.2)) destabilised on an equilibrium profile $\psi _0=\cos (x/a)$, whose $\varDelta '(k)$ formula has been discussed by Ottaviani & Porcelli (Reference Ottaviani and Porcelli1993, Reference Ottaviani and Porcelli1995) and for which the scaling $\gamma _M\sim S^{-1/2}$ can be deduced (dashed line in the figure).

Figure 6.

Dispersion relation of the purely resistive tearing as a function of $ka$ (dots correspond to numerically computed values), evaluated for the equilibrium $\psi _0=\cos (x/a)$. The smaller box is a ‘zoom’ of the plot for low values of $\gamma$, which correspond to the range of smaller values of $S^{-1}$ indicated in the caption. The dashed line corresponds to the fastest mode scaling, which for $\psi _0=\cos (x/a)$ is $\gamma _M\sim S^{-1/2}$. This can be obtained by using in (24) of Del Sarto et al. (Reference Del Sarto, Pucci, Tenerani and Velli2016) the value $p=2$, which can be deduced from the $ka\ll 1$ limit of the corresponding $\varDelta '(ka)$ formula (see Ottaviani & Porcelli Reference Ottaviani and Porcelli1993 or (8) of Ottaviani & Porcelli Reference Ottaviani and Porcelli1995). An 8th-order accurate scheme has been used here for integration on a uniform grid by using the solver discussed by Betar et al. (Reference Betar, Del Sarto, Ottaviani and Ghizzo2020).

When tearing modes are destabilised in a large aspect ratio, periodic current sheet, so that also the $\varDelta '\delta \sim 1$ modes are excited, several magnetic islands emerge and mutually interact during the nonlinear evolution, e.g. via the coalescence instability (cf. § 2.1), by moving along the neutral line. This dynamics is due to the superposition of several unstable modes with different wavenumbers and growth rates, among which the fastest mode is obviously the dominant one. In general, an aspect ratio $L/a$ of the order of ${\sim }20$ is sufficient to grant the instability of a wavenumber sufficiently close to the fastest growing mode Velli & Hood (Reference Velli and Hood1989) – see also figure 1 of Betar et al. (Reference Betar, Del Sarto, Ottaviani and Ghizzo2020). In this scenario, the ‘soliton-like’ behaviour of magnetic islands, noted already by Biskamp (Reference Biskamp1982) in his early simulations, has earned them the name of ‘plasmoids’ (see, e.g. Biskamp & Welter Reference Biskamp and Welter1989; Biskamp Reference Biskamp1996; Shibata & Tanuma Reference Shibata and Tanuma2001; Loureiro et al. Reference Loureiro, Schekochihin and Cowley2007; Bhattacharjee et al. Reference Bhattacharjee, Huang, Yang and Rogers2009), by analogy with the similarly behaving ‘plasma blobs’, also termed ‘plasmoids’, that have been identified and referenced in astrophysical plasma literature since the end of the sixties (Cowling Reference Cowling1967). Plasmoid structures displaying ‘soliton-like’ features have been observed in space plasmas, in particular, in magnetotail reconnection processes (see, e.g. early works such as those by Hones (Reference Hones1979) and Birn & Hones (Reference Birn and Hones1981) or later works such as Hautz & Scholer (Reference Hautz and Scholer1987), Scholer (Reference Scholer1987) and also Lin, Cranmer & Farrugia (Reference Lin, Cranmer and Farrugia2008), for a review of the usage of the term plasmoid in this context) and – although in a different geometrical setting – they have been also observed in connection to coronal mass ejections, arguably induced by reconnection processes in the solar corona (see, e.g. early works such as Pneuman Reference Pneuman1983; Gosling & McComas Reference Gosling and McComas1987). In this sense, in the recent literature, the term ‘plasmoid regime’ of reconnection (Uzdensky, Loureiro & Schekochihin Reference Uzdensky, Loureiro and Schekochihin2010; Huang & Bhattacharjee Reference Huang and Bhattacharjee2010; Comisso et al. Reference Comisso, Huang, Lingam, Hirvijoki and Bhattacharjee2018; Zhou et al. Reference Zhou, Wu, Loureiro and Uzdensky2021) is typically used to identify both tearing reconnection in the $\varDelta '\delta \sim 1$ wavelength limit (cf. the early numerical study by Biskamp (Reference Biskamp1986) on the instability of a Sweet–Parker current sheet to tearing modes) and the generation of magnetic islands in current sheets developed as a consequence of a turbulent motion. Here, the formation of ‘plasmoids’ is observed, either numerically (see Biskamp & Welter Reference Biskamp and Welter1989; Biskamp & Bremer Reference Biskamp and Bremer1994; Wei et al. Reference Wei, Hu, Feng and Schwen2000; Servidio et al. Reference Servidio, Matthaeus, Shay, Cassak and Dmitruk2009, Reference Servidio, Dmitruk, Greco, Wan, Donato, Cassak, Shay, Carbone and Matthaeus2011; Wan et al. Reference Wan, Matthaeus, Servidio and Oughton2013; Franci et al. Reference Franci, Cerri, Califano, Landi, Papini, Verdini, Matteini, Jenko and Hellinger2017 and several more recent works) or by experimental measurements (see, e.g. Nishizuka et al. Reference Nishizuka, Karlicky, Janvier and Bárta2015; Kumar et al. Reference Kumar, Karpen, Antiochos, Wyper and DeVore2019; Yan et al. Reference Yan, Xue, Jiang, Priest, Kliem, Yang, Wang, Kong, Song and Feng2022). We recall in this regard that also the term ‘plasmoid instability’ has been used in recent literature, with reference to the tearing instability associated with the fastest growing mode, which is destabilised on a Sweet–Parker current sheet of length L, and for which a scaling $\gamma _M\sim S_{L}^{1/4}$ is obtained (Tajima & Shibata Reference Tajima and Shibata1997; Loureiro et al. Reference Loureiro, Schekochihin and Cowley2007; Bhattacharjee et al. Reference Bhattacharjee, Huang, Yang and Rogers2009). The diverging scaling with respect to a vanishing non-ideal parameter – otherwise impossible for a spontaneous reconnecting instability – is here due to the fact that the fastest growing mode develops in this case as a ‘secondary’ reconnection process to a primary steady reconnection scenario (cf. figure 2): for $S_{L}^{-1}=0$, the Sweet–Parker reconnection does not occur, by thus forbidding the onset of the tearing modes. Due to the violently unstable nature of the tearing mode developing in the Sweet–Parker scenario, the possibility to measure, in Nature, a steady, Sweet–Parker-like reconnecting current sheet has been first questioned and discussed by Pucci & Velli (Reference Pucci and Velli2014).

Interpreting the magnetic islands observed in current sheets developed by turbulence as due to tearing mode instabilities is an hypothesis which is grounded on the assumption that a standard Fourier representation can be used in a Wentzel–Kramers–Brillouin (WKB)-sense also for reconnecting modes in a non-periodic, large aspect ratio current sheet. In this case, the fastest growing mode can be assumed to be, in an asymptotic sense, the dominant reconnecting instability. Although a proof of the correctness of this assumption for finite length current sheets has not been provided yet, this kind of approach has attracted since the beginning much attention in the context of both space plasmas (see, e.g. Cross & van Hoven Reference Cross and van Hoven1971; Van Hoven & Cross Reference Van Hoven and Cross1971) and solar physics (see, e.g. Priest Reference Priest1976; Velli & Hood Reference Velli and Hood1989), and has been more recently been considered for applications to secondary instabilities to primary reconnection processes (see, e.g. Tajima & Shibata Reference Tajima and Shibata1997; Shibata & Tanuma Reference Shibata and Tanuma2001; Tanuma et al. Reference Tanuma, Yokoyama, Kudoh and Shibata2001; Loureiro et al. Reference Loureiro, Cowley, Dorland, Haines and Schekochihin2005; Daughton & Scudder Reference Daughton and Scudder2006; Drake et al. Reference Drake, Swisdak, Schoeffler, Rogers and Kobayashi2006; Loureiro et al. Reference Loureiro, Schekochihin and Cowley2007; Bhattacharjee et al. Reference Bhattacharjee, Huang, Yang and Rogers2009; Landi et al. Reference Landi, Del Zanna, Papini, Pucci and Velli2015; Tenerani et al. Reference Tenerani, Velli, Rappazzo and Pucci2015; Del Sarto et al. Reference Del Sarto, Pucci, Tenerani and Velli2016; Del Sarto & Ottaviani Reference Del Sarto and Ottaviani2017; Papini, Landi & Del Zanna Reference Papini, Landi and Del Zanna2019b; Singh et al. Reference Singh, Pucci, Tenerani, Shibata, Hillier and Velli2019) and to reconnection in turbulence (see, e.g. Pucci & Velli Reference Pucci and Velli2014; Loureiro & Uzdensky Reference Loureiro and Uzdensky2016; Tenerani et al. Reference Tenerani, Velli, Pucci, Landi and Rappazzo2016; Comisso et al. Reference Comisso, Huang, Lingam, Hirvijoki and Bhattacharjee2018; Papini et al. Reference Papini, Franci, Landi, Hellinger, Verdini and Matteini2019a; Betar et al. Reference Betar, Del Sarto, Ottaviani and Ghizzo2020; Kowal et al. Reference Kowal, Falceta-Gonçalves, Lazarian and Vishniac2020; Pucci et al. Reference Pucci, Velli, Shi, Singh, Tenerani, Alladio, Ambrosino, Buratti, Fox and Jara-Almonte2020; Schekochihin Reference Schekochihin2020; Tenerani & Velli Reference Tenerani and Velli2020; Franci et al. Reference Franci, Papini, Del Sarto, Hellinger, Burgess, Matteini, Landi and Montagud-Camps2022) in the context of the so-called ‘turbulent (or turbulence-mediated or turbulence-driven) reconnection’ scenario (Matthaeus & Lamkin Reference Matthaeus and Lamkin1986; Strauss Reference Strauss1988; Loureiro et al. Reference Loureiro, Uzdensky, Schekochihin, Cowley and Yousef2009; Matthaeus & Velli Reference Matthaeus and Velli2011; Schekochihin Reference Schekochihin2020) – not to be confused with the study of ‘turbulent-driven magnetic island’ in tokamaks (cf. end of § 2.1). This topic is not of concern here, where we consider strictly periodic current sheets only. We therefore direct the interested reader to the aforementioned references. A survey of the scalings of the fastest growing mode in different reconnection regimes (among which the collisionless regimes we consider in the following) for periodic current sheets can be found in (Betar et al. Reference Betar, Del Sarto, Ottaviani and Ghizzo2020).

In the remainder of the article, we just focus on the formal solution of the eigenvalue problems in a slab current sheet in the $\varDelta '\delta \ll 1$ and $\varDelta '\delta \gg 1$ limits. In this regard, we note that, in a few cases, a closed form of the integral (3.7) has been provided: in resistive and collisionless regimes, formulae involving integration in the Fourier space have been provided by Pegoraro & Schep (Reference Pegoraro and Schep1986), Pegoraro et al. (Reference Pegoraro, Porcelli and Schep1989), Schep et al. (Reference Schep, Pegoraro and Kuvshinov1994), Basu & Coppi (Reference Basu and Coppi1981) and Porcelli (Reference Porcelli1991). These calculations provide the available analytical formulae for the growth rates in the ‘internal-kink’ and ‘constant-$\psi$’ RMHD, collisionless regimes.

Table 1 summarises the known results of the boundary layer analysis for the regimes of reconnection of interest in this article. The general dispersion relation for RMHD tearing modes, from which the large- and small-$\varDelta '$ limits can be obtained, has been written as in Del Sarto & Ottaviani (Reference Del Sarto and Ottaviani2017, appendix C) in terms of a characteristic scale length $\delta _L$ that coincides with the relevant layer width $\delta$ in the large- and small-$\varDelta '$ regimes. In §§ 5–7, we will detail the analytical calculations that allow the relevant limits of these dispersion relations to be recovered.

Table 1.

Known scalings available in the previous literature as obtained from a boundary layer analysis in the collisionless and resistive regimes. The general dispersion relation of Porcelli (Reference Porcelli1991) and Ottaviani & Porcelli (Reference Ottaviani and Porcelli1995) for $d_e^2\gg \rho _s^2$ and that of Porcelli (Reference Porcelli1991) and Comisso et al. (Reference Comisso, Grasso, Waelbroeck and Borgogno2013) for $\rho _s>d_e\neq 0$ are reported. In both cases, we have adopted the notation used in Del Sarto & Ottaviani (Reference Del Sarto and Ottaviani2017, appendix C), from which the small and large wavelength limits can be obtained once $\delta _L\to \delta$ is assumed and some estimation is made for $\delta =\delta _{LD}$ or $\delta =\delta _{SD}$. The general dispersion relation valid in all wavelength limits of the cold-electron resistive regimes (Ara et al. Reference Ara, Basu, Coppi, Laval, Rosenbluth and Waddell1978) and of the warm-electron resistive regimes (Pegoraro & Schep Reference Pegoraro and Schep1986) is not reported explicitly but it can be obtained from the collisionless cases previous substitution $d_e\to (S\gamma )^{-1/2}$. Note that, although available since the early solutions of the boundary layer calculations, the scalings of $\delta _{LD}$ and of $\delta _{SD}$ we have reported here in the different regimes have not been always pointed out in the related, reference works. These scalings correspond indeed to the characteristic width of the innermost boundary layer (i.e. $\delta _1$) in each reconnection regime, but the identification of the latter as the reconnecting layer width requires indeed further hypotheses which are discussed in § 8.

Finally, we draw attention to the fact that the definition (3.8) in a slab geometry changes meaning as soon as the condition $\varDelta '\delta \gtrsim 1$ is achieved, since application of (3.11) for $ka\ll 1$ and $\varDelta '\delta \sim 1$ formally implies a dependence of $k$ on $\delta$. In general, indeed, $(\varDelta ')^{-1}$ becomes a microscopic scale as soon as $\varDelta 'L_0\gg 1$, and for $\varDelta '\delta \gtrsim 1$, the condition $(\varDelta '(k))^{-1}\ll \delta$ is satisfied by a range of wavenumbers that depend on the non-ideal parameters which define $\delta$.

3.3 Boundary layer integration of collisionless tearing in the previous literature

In the following, we will discuss the method that can be used to find the analytical solution of the boundary layer problem associated with the scalings of table 1.

Differently from Pegoraro & Schep (Reference Pegoraro and Schep1986) and Porcelli (Reference Porcelli1991), who tackled the boundary layer calculations in the Fourier space, we are going to address the integration in the coordinate space: although a little more cumbersome, if one's interest is limited to obtaining the eigenvalue scaling alone, the analysis performed in the coordinate space allows a more intuitive understanding of the physics of the problem. Moreover, it spares one from the sometimes non-trivial task of reversing the Fourier transform so as to obtain the eigensolutions in the coordinate space. To the best of our knowledge, these kinds of calculations in the coordinate space for the warm-electron regime ($\rho _s\gtrsim d_e$) have never been reported in the literature, before.

Only in the work by Zocco & Schekochihin (Reference Zocco and Schekochihin2011) have some details of the analysis for the double boundary layer integration in the presence of two matching layers been discussed in the coordinate space (see appendix B therein). However, in that work, a different analytical approach has been taken for the integration of a different analytical reduced-MHD model including ion-FLR effects in the warm-collisionless regime. In particular, starting from gyrokinetic equations, these authors derived a reduced model consisting of three equations corresponding: to the scalar potential related to the ${\boldsymbol E} \times {\boldsymbol B}$ drift velocity ($\varphi$); to the parallel component of the vector potential ($A_{\|} = -\psi$); and to the ‘reduced’ electron distribution function ($g_e$, therein), which accounts for the moments of the electron distribution function, except for the density and the mean parallel velocity. The equation describing the evolution of $g_e$ is there coupled to the other two equations via the perturbation of the parallel electron temperature ($\delta T_{e,\|}$). In the linear limit, the latter can be written in terms of $\varphi$: this allowed the authors to obtain the customary system of two slab-geometry tearing equations for $\varphi$ and $\psi$, although modified with respect to (2.1)–(2.2) we consider here, because of the contribution of the gyrokinetic ion operator that appears in the vorticity equation after integration of the ion gyrokinetic equation (Pegoraro & Schep Reference Pegoraro and Schep1981); this contribution clearly vanishes in the limit of null ion temperature in which we are interested. In this way, Zocco and Schekochihin looked for asymptotic analytical solutions of the boundary layer problem in the coordinate space by identifying two distinct layers: the electron and the ion regions whose widths are defined by the reconnection layer width $\delta$ and by the ion-sound Larmor radius ($\rho _s$), respectively. However, differently from the approach we are going to develop below, in which we will look for a direct integration of the asymptotic solutions in the warm-collisionless regime, these authors solved the linear equations in the two asymptotic regions for the warm-collisionless regime by using perturbation methods (Zakharov & Rogers Reference Zakharov and Rogers1992). By matching the asymptotic solutions, they then recovered the scaling laws first evaluated in the Fourier space (Porcelli Reference Porcelli1991).

In the remainder of the article, we will treat both the collisionless and resistive limits in a unified way, by separately considering the warm ($\rho _s\gtrsim d_e$) and cold ($d_e\gg \rho _s$) regimes. Then, we will specify the results of the dispersion relation in each reconnection regime. We will thus reobtain the resistive scaling laws first computed by Furth et al. (Reference Furth, Killeen and Rosenbluth1963) in the cold-electron regime and by Pegoraro & Schep (Reference Pegoraro and Schep1986) in the warm-electron regime, and the collisionless scaling laws first evaluated by Coppi (Reference Coppi1964c,Reference Coppia) in the cold-electron limit and by Porcelli (Reference Porcelli1991) in the warm one.

In Appendix C, we make a brief review of previous works which have approached boundary layer calculations in the Fourier space (included are some more recent works which explicitly made a comparison between eigenmode solutions in the real and in the Fourier space – Connor et al. Reference Connor, Hastie and Zocco2012b), and we compare some key points of that procedure with the integration procedure that we follow in this work.

4.1 A specific example of evaluation of $\varDelta '(\psi _0;k)$

Having in mind the numerical results to be presented in the next sections, let us solve the outer equations for the equilibrium profile

(4.3)\begin{equation} \psi'_0(x) = \frac{\tanh(x)}{\cosh^2(x)}, \end{equation}

which is the one used by Betar et al. (Reference Betar, Del Sarto, Ottaviani and Ghizzo2020) to obtain the scalings to which we will make comparison. This profile was first proposed by Porcelli et al. (Reference Porcelli, Borgogno, Califano, Grasso, Ottaviani and Pegoraro2002). Substituting the previous relation into (4.2), one finds

(4.4)\begin{equation} \psi''_{\text{out}} - \left(\alpha^2 - \frac{12}{\cosh^2(x)}\right) \psi_{\text{out}} = 0, \end{equation}

where $\alpha ^2 = k^2 + 4$ and ${\psi ^{'''}_0}/{\psi ^{'}_0} = 4 - {12}/{\cosh ^2(x)}$. Grasso first evaluated the corresponding $\varDelta '$, reported by Porcelli et al. (Reference Porcelli, Borgogno, Califano, Grasso, Ottaviani and Pegoraro2002), by looking for a polynomial solution expressed in powers of $\tanh (x)$ (D. Grasso, private communication 2017). Here, we formally revise the problem, and we reduce (4.4) to a Legendre-type equation that can be shown to be valid for different kinds of equilibria, which we will discuss elsewhere.

First, we notice that $\lim _{x\to \infty } \cosh ^{-2}(x)=0$. This represents a main requirement if one seeks to apply periodic boundary conditions on the problem by considering $\psi _0$ as a periodic function in $x$ with spatial period of infinite extension. Therefore, far enough from the neutral line where the reconnection event takes place, the hyperbolic term in (4.4) can be neglected and (4.4) becomes

(4.5)\begin{equation} \psi''_{f} - \alpha^2 \psi_f = 0, \end{equation}

whose solution is

(4.6)\begin{equation} \psi_f(x) = A e^{-\alpha x}, \end{equation}

with $A$ an integration constant. One then expects the solution of (4.4) to take the form

(4.7)\begin{equation} \psi_{\text{out}}(x) = \psi_f(x) f(x) = A e^{-\alpha x} f(x). \end{equation}

Substituting the previous equation in (4.4), one obtains

(4.8)\begin{equation} f'' - 2 \alpha f' + \frac{12}{\cosh^2(x)} f = 0. \end{equation}

Motivated by the fact that $\tanh '(x) =\cosh ^{-2}(x)$, we perform the change of variables $z=\tanh (x)$. Then, (4.8) for $f(z)$ reads

(4.9)\begin{equation} (1 - z^2) f'' - 2 (z + \alpha) f' + 12 f = 0, \end{equation}

where the ‘$'$’ refers to the derivation with respect to $z$. The neutral line is now at $z = \tanh (x) = 0$. When $\alpha =0$, (4.9) becomes a Legendre equation of order three (we recall that a Legendre equation of order $p$ reads $(1-z^2)f''-2zf'+p(p+1)f=0$ – see Abramowitz & Stegun Reference Abramowitz and Stegun1964, § 8.1.1). This equation has two regular singular points occurring at $z=\pm 1 \ (x \to \pm \infty )$, and one of its two linearly independent solutions is a polynomial of third degree. This reads

(4.10)\begin{equation} f = a + b z + c z^2 + d z^3, \end{equation}

where the coefficients can be found by substituting the polynomial solution into (4.9) and equating to zero the coefficients with equal powers of $z$. After obtaining these coefficients, using $z = \tanh (x)$ and (4.7), the outer solution of (4.2) becomes

(4.11)\begin{align} \psi_{\text{out}} = A e^{-\alpha |x|} \left\{ 1 + \frac{6\alpha^2 - 9}{\alpha (\alpha^2 - 4)} \tanh(|x|) + \frac{15}{(\alpha^2 - 4)} \tanh^2(|x|) + \frac{15}{\alpha (\alpha^2 - 4)} \tanh^3(|x|) \right\}. \end{align}

The equation in the outer region accepts both an even and an odd solution. The even solution is the one required for tearing modes (the magnetic island shape is symmetric with respect to the neutral line – cf. also (2.8)). This solution displays a discontinuity in $\psi '_{\text {out}}$ at $x = 0$. Taking the limit of the first derivative of (4.11),

(4.12)\begin{equation} \lim_{x\to\pm 0} \psi'_{\text{out}} = A \left \{{\mp} \sqrt{k^2 + 4} \pm \frac{6k^2 + 15}{k^2 \sqrt{k^2+4}} \right \}, \end{equation}

we obtain, using (3.8),

(4.13)\begin{equation} \varDelta' = 2 \frac{-k^4 + 2k^2 + 15}{k^2 \sqrt{k^2 + 4}}. \end{equation}

It follows that $\varDelta ' > 0$ when $k \in [0, \sqrt {5}]$, which represents the spectrum of unstable wavenumbers to tearing-type instabilities for the equilibrium profile given by (4.3).

The outer solution for the case in which the magnetic equilibrium is given by Harris’ profile $\psi _0'(x) = \tanh (x)$ (Harris Reference Harris1962) can be obtained in the same manner, as already noted by Furth (Reference Furth1963) (§ 4, therein). In this case, the equivalent of (4.9) is a Legendre equation of degree $p=2$. Therefore, $f(z) = a + b z$, which leads to an even solution $\psi _{\text {out}}(x) = \psi _0 e^{-k|x|} \{ k + \tanh (|x|) \}$. Following the same arguments as before, one recovers the result $\varDelta ' = 2 ({1}/{k} - k)$, where the instability condition $\varDelta ' > 0$ is met when $k\in [0,1]$ (Furth et al. Reference Furth, Killeen and Rosenbluth1963). The point we underline here is that the solution of the outer equation and the evaluation of $\varDelta '$ by splitting the problem in (4.5)–(4.6) can be, in principle, used for a quite large class of magnetic equilibria.

As anticipated, the condition $\varDelta '>0$ determines the spectrum of unstable wavenumbers for a given equilibrium profile. This can be formally seen as related to the change of concavity that the solution $\psi _1$ must undergo while moving from the outer to the inner region, for the singular eigenfunction obtained in the ideal-MHD limit to become ‘regular’, thanks to the presence of non-ideal effects that allow reconnection (Furth et al. Reference Furth, Killeen and Rosenbluth1963). In general, exponentially growing solutions of the linear problem are obtained when $\varDelta '>0$ (cf. figure 7).

Figure 7.

Profile of the eigenfunction $\psi _1$ in the purely resistive regime, computed via the numerical solver of Betar et al. (Reference Betar, Del Sarto, Ottaviani and Ghizzo2020) for $S^{-1}=10^{-5}$ and $k=10$ (green curve), $k=1.7$ (orange curve) and $k=0.1$ (blue curve). These three values are respectively exemplificative of: the stable regime ($\varDelta '=-19.19<0$, green curve); of the small-$\varDelta '$ wavelength limit ($\varDelta '=3.2766$, red curve); and of the large-$\varDelta '$ wavelength limit ($\varDelta '=48.9$, blue curve). In the stable case ($\varDelta '<0$), the eigenfunction does not change concavity as $x$ approaches the neutral line, whereas for the unstable tearing-type solutions obtained for $\varDelta '>0$, the value $c_0\equiv \lim _{x\to 0}\psi _1(x)$ decreases approaching zero as the value of $\varDelta '(k)$ increases approaching $+\infty$.

5.1 Differential equations in the non-ideal region

The fact non-ideal effects are important in a region which is microscopic with respect to the equilibrium shear scale $a$, allows us to use (3.2) so as to write $\psi _0'(x)\approx 2x C_2= x\, \psi _0''|_{x=0}\equiv x J_0$, where we have named $J_0 \equiv \psi _0''(x)|_{x=0}$ the value (normalised to $B_0$ and $L_0=a$) of the second derivative of the equilibrium flux function at the neutral line. Using the ordering given by (3.4), (5.2)–(2.1) become

(5.4)\begin{gather} \gamma \psi_1 - \gamma \bar{d_e}^2 \psi_1'' = {\rm i} k J_0 x \varphi_1 - {\rm i} k J_0 \rho_s^2 x \varphi_1'', \end{gather}

(5.5)\begin{gather} \gamma \varphi_1'' = {\rm i} k J_0 x \psi_1''. \end{gather}

It can be noted that a direct comparison between the relative weight of the parameters $\bar {d_e}^2$ and $\rho _s^2$ in (5.4) can be done after eliminating $\varphi _{1}''$ in (5.4) via (5.5). Doing so, (5.4) takes the form

(5.6)\begin{equation} \gamma \psi_1 = {\rm i} k J_0 x \varphi_1 + \gamma\left(\bar{d_e}^2 + \rho_s^2\frac{ k^2 J_0^2}{\gamma^2} x^2 \right)\psi_1''. \end{equation}

The coefficient multiplying $\psi _1''$ at the right-hand side of (5.6) can be read as a generalised, space-dependent, ‘Lundquist number’, i.e. it can be assimilated to a space-dependent resistivity/conductivity. (In more complex reconnection models based on a gyrokinetic modelling and including also ion-FLR effects, it can be shown (Cowley et al. Reference Cowley, Kulsrud and Hahm1986; Zocco & Schekochihin Reference Zocco and Schekochihin2011) that the generalised conductivity, which rules the collisionless reconnection process inside the innermost layer dominated by electron dynamics, takes a rational-polynomial form of the kind $\sigma _e= ({a+b (x/ \delta _1)})/({c+d (x/ \delta _1)^2 + e(x/ \delta _1)^4})$, where $\delta _1$ is the innermost reconnecting layer width and $a, b, c, d$ and $e$ are coefficients that depend on the plasma parameters. This form has been used by Connor et al. (Reference Connor, Hastie and Zocco2012b) to obtain, via Fourier-space integration, a general dispersion relation encompassing drift-tearing modes, kinetic Alfvén modes and the internal kink mode at low values of the plasma $\beta$ parameter.) To have magnetic reconnection, i.e. to allow for the existence of tearing-type unstable modes, the coefficient $\bar {d_e}^2$ must be non-null since $\rho _s^2$ alone does not allow the relaxation of the topological constraints forbidding the intersection and breaking of initially distinct magnetic lines (cf. § 2.1). Once the importance of the $\bar {d_e}^2$ term in the non-ideal region is established, one sees that the condition for which the last term in parenthesis of (5.6) is negligible in a subdomain of the integration domain, even for a non-vanishing value of $\rho _s$, is

(5.7)\begin{equation} \frac{\rho_s^2}{\bar{d_e}^2}\ll \left(\frac{\gamma}{ k x J_0}\right)^2. \end{equation}

This assumption will be later done in some regimes and will be heuristically verified.

5.2 Normalisation and distinction of the boundary layers

More in general, to solve the system (5.4)–(5.5) with the boundary layer approach, further hypotheses are done about the relative magnitude of the different contributions inside of the inner layer, and some auxiliary function (traditionally noted as $\chi$) relating $\varphi _1$ to $\psi _1$ via (5.5) is introduced, depending on the $\varDelta '$-regime, so as to bring the system to a single ordinary differential equation (ODE). In each regime, appropriate normalisation of the spatial scales can be therefore chosen so as to better identify the extension of the sub-intervals of the inner region, where some term dominates over or is negligible with respect to the others. This is why, in the following, in the reconnection regimes and $\varDelta '$ limits which will be of interest, at each time, we will perform several changes of variable based on different normalisation choices.

This is an important step in the boundary layer procedure, since it is at this level that one introduces the layer widths $\delta _1$ or $\delta _2$ (cf. figure 4), although only implicitly: in practice, a ‘stretched’ coordinate is introduced by referring the coordinate $x$ to a characteristic length, say $\ell$, which replaces $x$ in different subdomains of the non-ideal region. It is when the limits $|x|/\ell \gg 1$ or $|x|/\ell \ll 1$ will be taken for the purpose of finding solutions of the differential equations within the boundary layer approach that the length $\ell$ will be identified as the layer width, e.g. $\delta _2$ or $\delta _1$. Also note that the length $\ell$ may initially depend on some unknown parameter like the eigenvalue $\gamma$: its scaling, and therefore that of $\ell$ and subsequently that of the layer width to which it corresponds, will be therefore determined a posteriori, from the conditions on the solutions obtained by integrating and matching the equations for small and large values of these stretched variables. Note that since the initial choice of each length scale $\ell$ is essentially arbitrary, it must be a posteriori verified for consistency that, asymptotically, $\delta _1\ll \delta _2 \ll 1$ (when both $\delta _2$ and $\delta _1$ are expressed in units of $L_0=a$).

It is also worth stressing that the identification of each boundary layer is subordered to the identification of the negligibility of some term of the non-ideal equations. This means that we must have knowledge of the relative ordering of the non-ideal parameters at play (for example, had we wanted to include viscosity, a third scale length – let us say $\lambda _{\nu }$ – would have appeared. If the hypothesis $\lambda _{\nu }\gg \rho _s \gtrsim \bar {d}_e$ had been fulfilled, one could have a priori expected three subdomains in the non-ideal region: one where ${\nu }$ alone dominates, one where both $\nu$ and $\rho _s$ in principle dominate, and one where all three non-ideal parameters are in principle important. Then, one should evaluate if, in the two innermost regions, specific conditions can be satisfied so that the dominant contribution of one or two non-ideal terms alone can be isolated: this possibly refines further the subdomains of interest for the purpose of the integration. It is at this level that we finally define the boundary layers, with respect to which the matching conditions will be fixed), or, if this is not the case, specific assumptions must be done on them by separately considering the different possible combinations, until in each case, the ‘thinnest’ sub-domain is this way singled out: this is the innermost region where the eigenmode solution must be first integrated by imposing the boundary conditions from outside.

As an example of general interest for the calculations that we are going to develop next, we can consider a normalisation of (5.4)–(5.5) to an arbitrary length $\ell$, which we will specify later, in each case examined. This will also let us get rid of some ‘superfluous’ parameters in the equations. We thus define

(5.8)\begin{equation} \zeta = \frac{x}{\ell}, \quad \mathcal{G}=\frac{\gamma}{k \ell J_0},\quad \tilde{\varphi}_1 ={-}{\rm i}\frac{\gamma }{k \ell J_0} \varphi_1(\zeta), \quad \tilde{\psi}_1 = \psi_1(\zeta), \end{equation}

which allow us to re-write (5.4)–(5.5) as

(5.9)\begin{gather} \tilde{\psi}_1 - \left(\frac{\bar{d}_e^2}{\ell^2} + \frac{\rho_s^2}{\ell^2} \frac{\zeta^2}{\mathcal{G}^2}\right)\tilde{\psi}_1'' ={-}\frac{\zeta}{\mathcal{G}^2}\, \tilde{\varphi}_1, \end{gather}

(5.10)\begin{gather} \tilde{\varphi}_1'' = \zeta \tilde{\psi}_1'', \end{gather}

where both functions $\tilde {\varphi }_1$ and $\tilde {\psi }_1$ depend now on $\zeta$.

We note that $\mathcal {G}$ depends on $\gamma$, which generally makes both quantities complex numbers: they are real for usual tearing-type instabilities, which, in the absence of diamagnetic effects related to equilibrium density gradients, do not propagate; instead, they are purely imaginary if the magnetic profile is stable. Because of this feature, in the following, it will be useful to perform some integration in the complex plane.

Also note that if $\ell =1$ (in units of $L_0=a$), then, in dimensional units, we can write $\mathcal {G} =\gamma \tau _0/(ka)$, where $\gamma \tau _0$ expresses the growth rate of the unstable mode measured with respect to the ‘natural timescale’ $\tau _0\sim \sqrt {mn_0/(4{\rm \pi} )}(c/J_0)$, which equals here the transition time of a shear Alfvén wave across the shear length $a$ (cf. Betar et al. Reference Betar, Del Sarto, Ottaviani and Ghizzo2020). Therefore, $\mathcal {G}|_{\ell =1} \ll 1$ in an asymptotic sense, as long as the wavenumber is fixed and then $ka$ is unordered with respect to $\gamma \tau _0$. If, instead, the scale $\ell$ is microscopic, the ordering of $\mathcal {G}$ with respect to unity must be a posteriori evaluated.

5.3 Auxiliary function ‘$\chi$’ and large- and small-$\varDelta '$ limits

The strategy by which the inner equations are integrated depends on mathematical features that we are going to discuss in detail in each specific regime and wavelength limit on which they depend (see also § 5.5). In general, some approximations and hypotheses are required, since they allow us to combine (5.9)–(5.10) into a single equation. An integration strategy which is particularly efficient, especially in the large-$\varDelta '$ limit, consists in casting the equations for $\tilde{\psi}_1$ and $\tilde {\varphi }_1$ into an equation for an auxiliary variable, even with respect to $\zeta$, which is typically noted $\chi (\zeta )$, after the notation first chosen for it by (Coppi et al. Reference Coppi, Galvao, Pellat, Rosenbluth and Rutherford1976; Ara et al. Reference Ara, Basu, Coppi, Laval, Rosenbluth and Waddell1978)

(5.11)\begin{equation} \chi(\zeta) \equiv \zeta \tilde{\psi}_1'(\zeta) - \tilde{\psi}_1(\zeta). \end{equation}

This is directly connected to (5.10), of which it is an integral. Indeed, we have

(5.12)\begin{equation} \tilde{\varphi}_1''(\zeta) = \chi'(\zeta) \quad \implies \quad \tilde{\varphi}_1'(\zeta) = \chi(\zeta) - \lim_{\zeta\to\infty}\chi(\zeta) = \chi(\zeta) - \chi_\infty. \end{equation}

In the last passage, we have named $\chi _\infty$ the value of $\chi (\zeta )$ as $\zeta \to \infty$. Note that in boundary layer calculations, when $\ell$ will be identified as $\delta _2$ or $\delta _1$, the limit $\zeta \to \infty$ for which $\chi _\infty$ is defined must be taken in the matching layer outside of the innermost layer, but still inside the non-ideal region, i.e. for $\delta _1\ll x\ll 1$ (in units of $L_0=a$): since $\tilde {\varphi }_1'$ represents the $y-$component of the velocity, it approaches zero as $\zeta \to \infty$ inside the ideal region. Determining the asymptotic behaviour of $\chi _\infty$ at both limits of the unstable wavelength spectrum is important since the scaling laws of the eigenvalue problem depend on it (Ara et al. Reference Ara, Basu, Coppi, Laval, Rosenbluth and Waddell1978). In practice, depending on whether $\chi _\infty$ is zero or not, the differential equation for $\chi$ is homogeneous or not (see (5.19) below). These topics are discussed in the remainder of this section by considering, for simplicity, the case of a single boundary layer (i.e. assuming $\ell =\delta _1$).

5.3.1 Asymptotic estimate of $\chi _\infty$ in the small-$\varDelta '$ limit

In the case of the small-$\varDelta '$ limit and when $x \to 0$, according to the definition of $\chi (\zeta )$ and to the matching condition $\lim _{\zeta \to \infty }\tilde {\psi }_1(\zeta )\sim \lim _{x\to 0}\psi _{\text {out}}(x)$, one finds

(5.13)\begin{equation} \chi_\infty=\lim_{\zeta\to\infty}\chi(\zeta) \sim \lim_{x\to 0}x \psi_{\text{out}}'(0) - \psi_{\text{out}}(0) \simeq{-} c_0, \quad \text{(constant-$\psi$ ordering)}. \end{equation}

In this case, we have used the ‘constant-$\psi$’ approximation (Furth et al. Reference Furth, Killeen and Rosenbluth1963), i.e. the fact that in the small-$\varDelta '$ regime, the inner solution tends to a constant value $\psi _1(0)$ when it approaches the innermost region. Such constant value is approximatively obtained already from the $x\to 0$ limit of the outer solution $\psi _{\text {out}}(x)$, so that in this wavelength limit, we can write $\lim _{x\to 0}\psi _{\text {out}}(x) = c_0$.

The constant-$\psi$ ordering can be assumed and heuristically verified after integration of the eigenmode problem for the single boundary layer limit (cf. Furth et al. Reference Furth, Killeen and Rosenbluth1963; White Reference White1983), and in the small-$\varDelta '$ regime, it can be shown to be equivalent to the condition $\delta _1 \varDelta '\ll 1$. In this case, we can also estimate $\lim _{x\to 0}x \psi _1'(x)\sim \lim _{x\to 0} x\varDelta ' \psi _{\text {out}}(x)\sim \delta _1 \varDelta 'c_0$, where in the last passage, we have (over)estimated $\lim _{x\to 0}x$ by the extent of the innermost layer width: owing to the $\varDelta ' \delta _1\ll 1$ condition, the last passage of (5.13) is thus justified, whence one deduces

(5.14)\begin{equation} \chi_\infty ={-}c_0, \quad \text{for } \delta_1 \varDelta' \ll 1. \end{equation}

5.3.2 Asymptotic estimate of $\chi _\infty$ in the large-$\varDelta '$ limit

In the large-$\varDelta '$ limit for $x\to 0$, instead, one finds that $\chi _\infty =0$. This has been first discussed by Coppi et al. (Reference Coppi, Galvao, Pellat, Rosenbluth and Rutherford1976) and Ara et al. (Reference Ara, Basu, Coppi, Laval, Rosenbluth and Waddell1978), where it has been shown that the vanishing of $\chi (\zeta )$ far out from the resonant region (equivalent to the innermost layer in our notation) is consistent with the rigid plasma displacement that characterises the internal kink mode in a cylindrical tokamak: in the notation and slab geometry assumption we use in this work, this corresponds to write $\tilde {\varphi }_1(\zeta )=-\gamma \zeta$ for $|\zeta |\leq 1$ and $\tilde {\varphi }_1(\zeta )=0$ for $|\zeta |> 1$ (cf. also Porcelli Reference Porcelli1987; Del Sarto & Ottaviani Reference Del Sarto and Ottaviani2017), whence the condition $\chi _\infty =0$ immediately follows by direct comparison with the definition in (5.12). An alternative way to look at the consistency of this result is to consider the Taylor expansion of the outer solution given by (4.7):

(5.15)\begin{equation} \psi_{\text{out}}(x) \sim c_0 + c_1 |x| \implies \psi_{\text{out}} \sim c_0 + c_1 |x| \implies \varDelta' \sim \frac{2c_1}{c_0}. \end{equation}

This solution formally holds in the outermost matching layer, where the instability parameter $\varDelta '$ is defined by neglecting non-ideal effects, i.e. by taking the limit $x\to 0$ of the eigenmode in the ideal region. Therefore, the use of $\lim _{x\to 0}\psi _{\text {out}}(x)$ to match the solution $\lim _{\zeta \to \infty }\psi _{1}(\zeta )$ is a priori not always justified. This however is not the case of the large-$\varDelta '$ limit, in which it can be a posteriori shown, from heuristic arguments (see Ottaviani & Porcelli (Reference Ottaviani and Porcelli1995) for the ‘cold’ reconnection regime) or by numerical integration, that the discontinuity in the derivative of $\psi _{\text {out}}$ occurs in a position which gets progressively closer to the neutral line, as long as the numerical value of $\varDelta '(k)$ increases (cf. also §§ 8–10). In this limit, using therefore $\lim _{\zeta \to \infty }\tilde {\psi }_1\simeq \lim _{x\to 0}\psi _{\text {out}}(x)$ combined with (5.15) in the definition of (5.11), one finds

(5.16)\begin{equation} \lim_{\zeta\to\infty}\chi(\zeta) \sim c_1 |x| - c_0 - c_1 |x| ={-}c_0 \sim \frac{c_1}{\varDelta'}. \end{equation}

The fact that $c_0\to 0$ as $\varDelta '$ increases is shown in figure 7, where the spatial profile of $\psi _1$, obtained after numerical integration of the purely resistive regime (for $S^{-1}=10^{-5}$), is shown for the values of $\varDelta '(k)= -19.19$ (green curve), $\varDelta '(k)= 48.9$ (blue curve) and $\varDelta '(k)= 3.2766$ (orange curve), evaluated using (4.13) for $k=10, 0.1$ and $1.7$, respectively. Thus, in the large-$\varDelta '$ limit where $\varDelta '\to \infty$, (5.16) implies that

(5.17)\begin{equation} \chi_\infty = 0 \quad \text{for } \varDelta' \to \infty. \end{equation}

The fact that condition $\varDelta ' \to \infty$ can be replaced by $\varDelta '$ larger than the inverse of some characteristic scale length will be discussed in § 10. As already anticipated (§ 3.2), this condition can be generally read as $\varDelta '\delta _1\gg 1$.

5.4 Auxiliary equation for the function $\chi (\zeta )$ in the non-ideal region

Definitions in (5.11)–(5.12) can be used to cast (5.9) in an equation for $\chi (\zeta )$. For analytical convenience, i.e. to facilitate the substitution of variables, it is preferable to evaluate the derivative of (5.9) after it has been divided by $\zeta$. After opportunely regrouping the terms of the equation so obtained, the result reads

(5.18)\begin{equation} (\zeta\tilde{\psi}_1' - \tilde{\psi}_1) - \frac{\bar{d}_e^2}{\ell^2} (\zeta\tilde{\psi}_1'''+\tilde{\psi}_1'' -2\tilde{\psi}_1'') -\frac{1}{\mathcal{G}^2}\frac{\rho_s^2}{\ell^2}\, \zeta^2\, (\zeta \tilde{\psi}_1'''+\tilde{\psi}_1'')={-}\frac{1}{\mathcal{G}^2}\,\zeta^2 \tilde{\varphi}_1', \end{equation}

which can be rewritten as

(5.19)\begin{equation} \left(\frac{\bar{d}_e^2}{\ell^2} +\frac{\rho_s^2}{\ell^2} \frac{\zeta^2 }{\mathcal{G}^2} \right) \chi'' -2 \,\frac{\bar{d}_e^2}{\ell^2} \frac{\chi'}{\zeta} - \left(1+ \frac{\zeta^2 }{\mathcal{G}^2}\right) \chi ={-}\frac{\zeta^2}{\mathcal{G}^2}\, \chi_{\infty}. \end{equation}

In the following, we will take different limits of this equation, depending on the tearing regime and wavelength limit considered. We emphasise that to give to the right-hand side term $\chi _\infty$ of (5.19) the geometrical (and physical) meaning we have previously discussed in terms of the boundary layer theory, the length $\ell$ must correspond to the layer width, $\delta _2$ or $\delta _1$ (depending on whether we are in the case of figure 4), to which the solution $\psi _{\text {out}}$ obtained in the ideal region is matched.

5.5 Wavelength limit and choice of integration via the equations for $\chi$ or for $\psi$ and $\varphi$

The choice of performing the integration by using (5.4)–(5.5) for $\psi _1$ and $\varphi _1$, rather than (5.19) for $\chi$, clearly depends on analytical convenience. In general, we note that this choice can be biased by the ease by which the boundary conditions imposed by the solution obtained in the outer, ideal region, are ‘transferred’ to the solution in the innermost region. The asymptotic matching generally poses some constraints on the integration constants in each subdomain, but it is easy to see that, operationally speaking, these conditions take different and quite ‘appealing’ forms in the small- and large-$\varDelta '$ limits, depending on whether one looks at the equation for $\psi _1$ or at the auxiliary equation for $\chi$, respectively. In particular, the following are observed.

1. (i) In the small- $\varDelta '$ limit, it is usually more convenient to perform the integration on the equation for $\psi _1$, for which the boundary conditions imposed from the ideal region are directly implemented by making explicit the dependence on $\varDelta '$ via (3.7) and (3.10). In this case, the constraint imposed by the solution $\psi _{\text {out}}$ in the ideal region can be directly used on the integration of the innermost equation in the form of the constant-$\psi$ condition. Indeed, if we define $c_0\equiv \lim _{x\to 0}\psi _{\text {out}}(x)$ and we heuristically assume the constant-$\psi$ condition to be valid, the latter implies that $\psi _1\simeq c_0$ both in the innermost region $|x|\leq \delta _1$ and in the possible intermediate region $\delta _1 < |x|\leq \delta _2$. Therefore, also when two boundary layers are present, the constant-$\psi$ condition can be directly implemented in the integration of the innermost equation, by thus practically allowing one to ‘bypass’ the procedure of matching between the solutions in the ideal and in the intermediate region, and then between the solutions in the intermediate and in the innermost regionFootnote ² (at least for the purpose of establishing the dispersion relation, i.e. to find the asymptotic scalings of the eigenvalue). However, as the matching occurs in the boundary layer with the ideal region, it is generally more convenient to maintain the normalisation to the intermediate layer width, also when two layers are present. This approach will be used in §§ 6.2 and 7.5.

2. (ii) In the large- $\varDelta '$ limit, it is instead usually more convenient to perform the integration on the auxiliary equation, since in this case, the constraint imposed by the outer solution $\psi _{\text {out}}$ on the solutions in the non-ideal region is simply expressed as $\chi _\infty = 0$ (cf. § 5.3.2). From an operational point of view, this realises the matching between the ideal and the intermediate region. The dispersion relation is therefore directly obtained by solving the non-ideal equation for $\chi _\infty = 0$ in case a single boundary layer is present, or by matching the innermost solution and the intermediate solution (both sought under the constraint $\chi _\infty = 0$) if two boundary layers are present. This approach will be used in §§ 6.1 and 7.4.

6.1 Solution for $\bar {d}_e^2 / \rho _s^2 \gg 1$: large-$\varDelta '$ limit

Assuming heuristically the validity of (5.7), we neglect the second term in parentheses in (5.9). The restriction to the large $\varDelta '$-limit, in which the condition in (5.17) holds, suggests that we make use of the auxiliary function $\chi (\zeta )$ defined by (5.11), which leads us to consider the appropriate limit of (5.19):

(6.1)\begin{equation} \chi''-\frac{2}{\zeta}\chi'- \frac{\ell^2}{\bar{d}_e^2}\left(1+\frac{\zeta^2}{\mathcal{G}^2}\right)\chi = 0. \end{equation}

The structure of the equation suggests that we take $\ell \equiv \bar {d}_e$ as the normalisation length to be used in the definitions of (5.8). This means that we postulate

(6.2)\begin{equation} \delta_1=\bar{d}_e, \end{equation}

which, at least for $\bar {d_e}=d_e$, is evidently consistent with the asymptotic condition $\delta _1\ll 1$. A solution to (6.1) can be sought in the form

(6.3)\begin{equation} \chi_1 = e^{-\alpha \zeta^2}. \end{equation}

Substituting it in (6.1) yields

(6.4)\begin{equation} 2\alpha = 1, \quad 4 \alpha^2 = \frac{1}{\mathcal{G}^2}. \end{equation}

Using the first part of (5.8) and combining the two conditions above, gives

(6.5)\begin{equation} \gamma = k {\bar{d}_e} J_0. \end{equation}

This result can be specialised both to the fully collisionless ($\bar {d}_e^2=d_e^2$) and to the fully collisional ($\bar {d}_e^2=S^{-1}/\gamma$) internal-kink regime. We recall in this regard that numerical analysis (Betar et al. Reference Betar, Del Sarto, Ottaviani and Ghizzo2020) confirms that the parameter space interval where both inertial and resistive effects non-negligibly combine is very narrow, so that tearing modes can be usually considered either as fully collisionless or as fully resistive.

In the collisionless limit, we thus recover the result of (Porcelli Reference Porcelli1991)

(6.6)\begin{equation} \gamma = k d_e J_0\,\quad \delta_1\sim d_e. \end{equation}

In the resistive case, (6.6) specialises to the result of (Coppi et al. Reference Coppi, Galvao, Pellat, Rosenbluth and Rutherford1976)

(6.7)\begin{equation} \gamma = k^{{2}/{3}} J_0^{{2}/{3}} S^{-{1}/{3}}, \quad \delta_1 \sim k^{-{1}/{3}} J_0^{-{1}/{3}} S^{- {1}/{3}}. \end{equation}

Substitution of the results of (6.6) and (6.7) for $x\sim \delta _1$ in (5.7) makes it equivalent to the condition $\rho _s^2/\bar {d}_e^2\ll 1$ from which we started, thus verifying a posteriori its correctness.

The scalings in (6.6) and (6.7)), in contrast to those of the large-$\varDelta '$ limit of the warm-collisionless regime, can be obtained using a heuristic-type derivation based on dimensional estimates, as we are going to prove in § 9.3.

Finally, some words about the second solution of (6.1), which we have neglected, since it can be shown that it is not of interest to the tearing instability problem. To this purpose, we can look for a solution of the form $\chi _2(\zeta ) = u(\zeta )\, \chi _1(\zeta )$. Substituting $\chi =\chi _2$ in (6.1) and using the fact that $\chi _1$ already solves it, one obtains, after a little algebra,

(6.8)\begin{equation} u''\chi_1 + 2 u'\left(\chi_1' -\frac{\chi_1}{\zeta}\right)+u\chi_1''=0. \end{equation}

Further substitution of (6.3) into the equation above yields

(6.9)\begin{equation} u''-\left(4\alpha \zeta + \frac{2 }{\zeta}\right)u'=0, \end{equation}

which, defining $w = u'$, can be brought to a directly integrable form:

(6.10)\begin{equation} \frac{{\rm d}w}{w} = 4 \alpha \zeta \,{\rm d}\zeta + \frac{2 \,{\rm d}\zeta}{\zeta} \implies w(\zeta) = w_0\, \zeta^2 e^{2\alpha\zeta^2}, \end{equation}

$w_0$ being a constant of integration. Specialising the result for $\alpha =1/2$, integrating $w$ once to find $u$ and then substituting the result in $\chi _2 = u \chi _1$, one obtains

(6.11)\begin{equation} \chi_2 = {\rm i} \sqrt{\frac{\rm \pi}{2}} {\rm erf}\left( \frac{{\rm i} \zeta}{\sqrt{2}} \right) e^{-({\zeta^2}/{2})} + \zeta, \end{equation}

which is odd and unbounded when $\zeta \to \infty$. Therefore, for our purposes, the solution $\chi _2$ can be disregarded since its spatial parity is not compatible with that of tearing modes (i.e. it is not compatible with the formation of $X$-points).

6.2 Solution for $d_e^2/ \rho _s^2 \gg 1$: small-$\varDelta '$ limit

In the small-$\varDelta '$ limit, the boundary layer analysis mirrors that of the paper by Furth et al. (Reference Furth, Killeen and Rosenbluth1963) if the resistive limit $\bar {d}_e^2\to S^{-1}/\gamma$ is considered. In general, regardless of the form of $\bar {d}_e$, we can here use the ‘constant-$\psi$’ approximation stating that $\tilde {\psi }_1(\zeta )\to c_0$ as $\zeta \to 0$ (cf. (5.15)). However, instead of using the auxiliary equation (5.19), in this case, it is more convenient to consider (5.9) in which $\tilde {\psi }_1''$ has been eliminated by using (5.10) since, after combination with the constant-$\psi$ limit, this yields an equation for the function $\tilde {\varphi }_1(\zeta )$ that can be brought to an integrable form more easily than the equation for $\chi (\zeta )$:

(6.12)\begin{equation} \tilde{\varphi}_1'' - \frac{\zeta^2}{\mathcal{G}^2} \frac{\ell^2}{\bar{d}_e^2}\tilde{\varphi}_1 = c_0\frac{\ell^2}{\bar{d}_e^2}\, \zeta. \end{equation}

Note that in this equation, the length $\ell$ is so far left unspecified. In practice, at this stage of the calculations, it can be taken $\ell =1$, which would be still consistent with the constant-$\psi$ limit with respect to which the substitution $\tilde {\psi }_1(0)\simeq c_0$ has been done. The identification of a proper normalisation length defining the boundary layer width can be done after a further change of variable that leads (6.12) to an integrable form in which the boundary layer integration can be more easily performed. Equation (6.12) for $\ell =1$ can be indeed brought to the form

(6.13)\begin{equation} \varPhi'' - z^2 \varPhi = z \end{equation}

after the change of variables and coordinates

(6.14)\begin{equation} z=\left({\mathcal{G}}{\bar{d}_e}\right)^{-{1}/{2}}\zeta,\quad {\varPhi} = \left(\frac{\bar{d}_e}{\mathcal{G}^3}\right)^{{1}/{2}} \left( \frac{\tilde{\varphi}_1}{c_0}\right). \end{equation}

The normalisation of Eq. (6.12) is detailed in Appendix D as a didactical example of how to generally proceed to find the normalisation coefficients of a differential equation. In the limit $\bar {d}_e^2\to S^{-1}/\gamma$, (6.13) is the equation of the resistive tearing mode studied by Furth et al. (Reference Furth, Killeen and Rosenbluth1963), whose solution has been provided by Basu & Coppi (Reference Basu and Coppi1976), Coppi et al. (Reference Coppi, Galvao, Pellat, Rosenbluth and Rutherford1976) and Ara et al. (Reference Ara, Basu, Coppi, Laval, Rosenbluth and Waddell1978) (cf. (III.26) and (III.27) and appendix A in the latter) in the closed integral form:

(6.15)\begin{equation} \varPhi={-} \frac{1}{2}z\int_0^1 (1-t^2)^{-{1}/{4}}e^{-({1}/{2})tz^2}\,{\rm d}t. \end{equation}

The scaling of the eigenvalue $\gamma$ and its dependence on the amplitude of the instability parameter $\varDelta '$ can be made explicit by comparing the eigenfunction $\varPhi (z)$ to ${\varphi }_1(x/\delta _2)$ in the matching condition expressed by (3.7), which also allows us to identify the layer width $\delta _2$: combining (3.7) and (3.10) (where we re-introduce the quantity $J_0$) yields

(6.16)\begin{equation} \psi_{\text{out}}(0) \varDelta' ={-}{\rm i}\frac{\gamma}{kJ_0}\int_{-\bar{X}'}^{+\bar{X}'} \frac{{\varphi}_1''(x')}{x'}\, {{\rm d} x}', \end{equation}

where $\varDelta '$ is defined by (3.8), $x'=x/\delta _1$ and $\bar {X}'=\bar {X}/\delta _1$ is indicative of some point of the matching region such that $\delta _1 \ll \bar {X}\ll 1$ in units of $L_0=a$. Connection with (6.13) follows from noticing that (6.14) implies:

(6.17)\begin{equation} \left.\frac{{\rm d}^2\tilde{\varphi}_1}{{\rm d}\zeta^2}\right|_{\ell=1} = c_0 \left.\left(\frac{1}{\bar{d}_e\mathcal{G}^3}\right)\right|_{\ell=1}^{{1}/{2}} \left( \frac{{\rm d}z}{{\rm d}\zeta} \right)^2 \frac{{\rm d}^2 \varPhi}{{\rm d}z^2} = c_0 \left.\left(\frac{\mathcal{G}}{\bar{d}_e^3}\right)\right|_{\ell=1}^{{1}/{2}} \frac{{\rm d}^2 \varPhi}{{\rm d}z^2}. \end{equation}

This authorises us to identify

(6.18)\begin{equation} x'=\frac{x}{\delta_1} \longrightarrow z=\frac{x}{(\mathcal{G}\bar{d}_e)|_{\ell=1}^{{1}/{2}}}, \quad -{\rm i}\frac{\gamma}{kJ_0}\frac{{\rm d}^2{\varphi(x')}_1}{{{\rm d} x}'^2} \longrightarrow c_0 \left.\left(\frac{\mathcal{G}}{\bar{d}_e^3}\right)\right|_{\ell=1}^{{1}/{2}} \frac{{\rm d}^2 \varPhi}{{\rm d}z^2}, \end{equation}

which now fixes the boundary layer width to be

(6.19)\begin{equation} {\delta_1}=\left.(\mathcal{G}\bar{d}_e)\right|_{\ell=1}^{{1}/{2}} =\sqrt{\frac{\gamma\bar{d}_e}{k}}, \end{equation}

whose consistency with the condition $\delta _1\ll 1$ must be verified once $\mathcal {G}|_{\ell =1}$ is explicitly obtained. Substituting then $\psi _{\text {out}}(0)\simeq \lim _{\zeta \to 0}\tilde {\psi }_1(\zeta )=c_0$ and using (6.18) makes (6.16) become

(6.20)\begin{equation} \varDelta' = \left.\left(\frac{\mathcal{G}}{\bar{d}_e^3}\right)\right|_{\ell=1}^{{1}/{2}} \int_{-\bar{Z}}^{+\bar{Z}} \frac{\varPhi''}{z} \,{\rm d}z, \end{equation}

where $\bar {Z} = \bar {X}/ (\mathcal {G}\bar {d}_e)|_{\ell =1}^{1/2}$. The integral in the previous equation is convergent, let us call its value $I\equiv \int _{-\bar {\infty }}^{+\infty } (\varPhi ''/{z})\, {\rm d}z \simeq \int _{-\bar {Z}}^{\bar {Z}} (\varPhi ''/{z}) \,{\rm d}z$. Then, using the definition in (5.8) for $\ell =1$, one obtains

(6.21)\begin{equation} \gamma = k{\bar{d}_e}^3 (\varDelta')^2 \left(\frac{J_0}{I^2}\right). \end{equation}

In the inertia-dominated regime, $\bar {d}_e=d_e$, one recovers the well-known result by Coppi (Reference Coppi1964c,Reference Coppia), later re-obtained by Porcelli (Reference Porcelli1991):

(6.22)\begin{equation} \gamma = k (\varDelta')^2 {{d}_e}^3 \left(\frac{J_0}{I^2}\right), \quad \delta_1 = \frac{\varDelta' d_e^2}{I}. \end{equation}

In taking the resistive limit, $\bar {d}_e= (S^{ -1}/\gamma )^{1/2}$, one obtains instead the classical result of Furth et al. (Reference Furth, Killeen and Rosenbluth1963), later recovered by Pegoraro & Schep (Reference Pegoraro and Schep1986):

(6.23)\begin{equation} \gamma = k^{{2}/{5}}(\varDelta')^{{4}/{5}} S^{-{3}/{5}} \left(\frac{J_0}{I^2}\right)^{{2}/{5}}, \quad \delta_1 = \frac{k^{-{2}/{5}}(\varDelta')^{{1}/{5}} S^{-{2}/{5}}}{({J_0}{I^2})^{{2}/{5}}}. \end{equation}

In both cases (6.22) and (6.23), the asymptotic condition $\delta _1\to 0$ as $d_e\to 0$ or $S^{-1}\to 0$ is verified.

7.1 Equations in the non-ideal region for $\rho _s\gtrsim d_e$

As we enter the non-ideal region from the ideal one, when $x< 1$, the first non-ideal characteristic scale that is encountered under the hypotheses previously done is $\rho _s > d_e$. This is the first scale length with respect to which we can try to identify a boundary layer: aiming to single out a subdomain in the $x\ll 1$ range where only $\rho _s$ possibly matters, we start by fixing $\ell =\rho _s$ in (5.8), so that (5.9) and (5.10) become

(7.1)\begin{equation} \tilde{\psi}_1 - \left(\frac{\bar{d}_e^2}{\rho_s^2} + \frac{\tilde{\zeta}^2}{\tilde{\mathcal{G}^2}}\right)\tilde{\psi}_1'' ={-}\frac{\tilde{\zeta}}{\tilde{\mathcal{G}^2}}\tilde{\varphi}_1,\quad \tilde{\varphi}_1'' = \tilde{\zeta} \tilde{\psi}_1'', \end{equation}

with

(7.2)\begin{equation} \tilde{\zeta}=\frac{x}{\rho_s}, \quad\tilde{\mathcal{G}}=\frac{\gamma}{k\rho_s J_0}, \quad \tilde{\varphi}_1={-}i\tilde{\mathcal{G}}{\varphi}_1(\tilde{\zeta}), \quad \tilde{\psi}_1={\psi}_1(\tilde{\zeta}). \end{equation}

The auxiliary equation (5.19) specialises to

(7.3)\begin{equation} \left(\frac{\bar{d}_e^2}{\rho_s^2} + \frac{\tilde{\zeta^2} }{\tilde{\mathcal{G}}^2} \right) \tilde{\chi}'' -2 \frac{\bar{d}_e^2}{\rho_s^2} \frac{\chi'}{\tilde{\zeta}} - \left(1+ \frac{\tilde{\zeta}^2 }{\tilde{\mathcal{G}}^2}\right) \tilde{\chi} ={-}\frac{\tilde{\zeta}^2}{\tilde{\mathcal{G}}^2}\, \tilde{\chi}_{\infty}. \end{equation}

Equation (7.3) can be shown to be equivalent to (B42) of Zocco & Schekochihin (Reference Zocco and Schekochihin2011), although the latter has been obtained from a four-field system which also accounts for the dynamics of the reduced electron guiding centre contribution and for the parallel electron temperature perturbation ($g_e$ and $\delta T_{\|,e}$, respectively, therein). These authors also solved the boundary layer problem in the coordinate space, but – also because of the different system of equations they started from – they took a different path with respect to the one we take here, in which we tackle the integration of the inner equation directly in the different wave length limits of interest. In this sense, the comparison between the boundary layer integration procedure of the aforementioned article and the one we are going to outline here is not immediate, although possible, since some different operational assumptions in the definition of the boundary layers and in the analytical technique of integration have been used in that work and in the present one. In their work, Zocco & Schekochihin (Reference Zocco and Schekochihin2011) identified two non-ideal integration regions where ion physics and electron physics were respectively dominating. In these two regions, they sought the solutions of their auxiliary equation, from which they computed the magnetic flux function $\psi$ in the electron and ion region ($A_{\|,e}$ and $A_{\|,i}$ in their notation, respectively, which, in this work, map into the solution $\psi$ that we will compute in the innermost and outermost non-ideal regions) by assuming in the large- and small-wavelength limit some heuristic orderings (a posteriori verified) in terms of $\varDelta '$ and what they named the ‘reconnection region’ ($\delta _{\text {in}}$, in their notation). As we will later discuss (§ 8), the characteristic width of their ‘reconnection region’, $\delta _{\text {in}}$, differs from that of the innermost boundary layer that we are going to identify ($\delta _1$). However, in some regimes, $\delta_{in}$ can be interpreted in the light of a further characteristic scale length, which we will introduce in § 10 (namely, what we will call $(\varDelta _{v_y}')^{-1}$). In the following, we are going to adopt an integration procedure based on the integral representation of hypergeometric functions, whereas, in their work, (Zocco & Schekochihin Reference Zocco and Schekochihin2011) used a perturbative method to obtain a closed form solution for $A_{\|,e}$. Note that, so far, the boundary layers have not yet been identified here.

7.2 Boundary layers in the non-ideal region for $\rho _s^2/d_e^2\gg 1$

We now specialise the condition $\rho _s\gtrsim d_e$ by making the assumption $\rho _s^2/\bar {d}_e^2\gg 1$, which requires a sufficiently large scale separation between $\bar {d}_e$ and $\rho _s\geq \bar {d}_e$. Note that this is henceforth the strongest assumption about the relative ordering of the non-ideal parameters at play, which is done in this analysis for tearing modes in the warm-electron regime. Nevertheless, the quantitative results obtained in this way are usually applied also for $\rho _s\gtrsim \bar {d}_e$: numerical calculations in a wide interval of the parameter space (Betar et al. Reference Betar, Del Sarto, Ottaviani and Ghizzo2020) suggest that the condition $\rho _s> \bar {d}_e/10$ can be assumed for the validity of the results we are going to obtain below. In the interval $\rho _s/10\lesssim \bar {d}_e \lesssim 10 \rho _s$, a transition from power-law scalings of the cold limit of § 6 to those of the warm limit that we are going to discuss here is observed.

7.2.1 Outermost boundary layer

Assuming $\rho _s^2/\bar {d}_e^2\gg 1$, (7.1) becomes

(7.4)\begin{equation} \tilde{\psi}_1 - \frac{\tilde{\zeta}^2}{\bar{\mathcal{G}^2}}\tilde{\psi}_1'' ={-}\frac{\tilde{\zeta}}{\bar{\mathcal{G}^2}}\tilde{\varphi}_1,\quad \tilde{\varphi}_1'' = \tilde{\zeta} \tilde{\psi}_1'' \end{equation}

and (7.3) reads

(7.5)\begin{equation} \tilde{\chi}'' - \left(1+\frac{\tilde{\mathcal{G}}^2 }{\tilde{\zeta}^2}\right) \tilde{\chi} ={-} \tilde{\chi}_\infty. \end{equation}

These equations must be assumed to be valid in an interval where spatial scales are comparable to $\rho _s$ but larger than $\bar {d}_e$, that is, $\bar {d}_e/\rho _s\lesssim \tilde {\zeta }\ll 1$, where the first inequality is meant in the sense discussed above. Therefore, restriction to this interval formally corresponds to take $\bar {d}_e = 0$ in (5.4). The integration of (7.8) and the matching of their solutions with $\psi _{\text {out}}$ of the ideal region via (3.5) and (3.6) formally identifies the outermost boundary layer of width

(7.6)\begin{equation} \delta_2 =\rho_s. \end{equation}

7.2.2 Innermost boundary layer

Of course, since the formation of the $X$-point via reconnection requires $\bar {d}_e\neq 0$, a second, innermost boundary layer can be defined in a region where both scales $\bar {d}_e$ and $\rho _s$ matter. The width $\delta _1$ of this second boundary layer does not necessarily correspond to the scale $\bar {d}_e$ (if it were so, the condition $\delta _1\ll \delta _2$ would require $\bar {d}_e\ll \rho _s$, which is an even stronger restriction than that defining the validity of (7.6)). To identify this inner boundary layer, we can look back to (7.1) and (7.3) or to ((5.9) and (5.10) and (5.19)) and notice that the two terms in parentheses, in which there is an explicit dependence on both $\rho _s$ and $\bar {d}_e$, are comparable (regardless of the definition of $\ell$ and $\mathcal {G}$) if $(\zeta \rho _s)^2/(\mathcal {G}\bar {d}_e)^2\sim 1$. In dimensional units, this means $x k\rho _sJ_0/(\gamma \bar {d}_e)\sim 1$, which suggests that we introduce the normalisation scale $\ell = {\gamma \bar {d}_e}/(k\rho _s J_0).$ Using then in (5.9) and (5.10) the change of coordinates and variables expressed by

(7.7)\begin{equation} \zeta=\frac{x }{\ell}= \frac{x k\rho_sJ_0}{\gamma \bar{d}_e}, \quad{\mathcal{G}}=\frac{\rho_s}{\bar{d}_e }, \quad {\varphi}_1={-}{\rm i}{\mathcal{G}}{\varphi}_1(\zeta), \quad {\psi}_1={\psi}_1(\zeta), \end{equation}

we write

(7.8)\begin{equation} {\psi}_1 - \frac{\bar{d}_e^2}{\ell^2}\left(1 + {\zeta^2}\right){\psi}_1'' ={-}\frac{\zeta}{{\mathcal{G}^2}}\varphi_1,\quad {\varphi}_1'' = \zeta {\psi}_1''. \end{equation}

Using the definition of ${\mathcal {G}}$ in (7.7) and the fact that $\rho _s^2/\bar {d}_e^2\gg 1$, one sees that in this interval, $\zeta {\varphi }_1 /{\mathcal {G}}^2\sim \zeta ^2 {\psi }_1 \bar {d}_e^2/(J_0^2)\rho _s^2 \ll {\psi }_1$ so that (7.8) can be approximated as

(7.9)\begin{equation} {\psi}_1 - \frac{\bar{d}_e^2}{\ell^2} \left(1 + \zeta^2\right) {\psi}_1'' =0,\quad {\varphi}_1'' = \zeta {\psi}_1'', \end{equation}

and with the same normalisation, (7.3) (or (5.19)) becomes

(7.10)\begin{equation} (1+\zeta^2)\chi''-\frac{2}{\zeta}\chi'-\frac{\ell^2}{\bar{d}_e^2}\chi={-}\zeta^2\chi_\infty. \end{equation}

Note that in the definitions above, we have used the same symbols for $\ell$, $\mathcal {G}$ and $\zeta$, previously introduced in (5.8), since the different regime of interest here should not induce any confusion with those symbols used in § 6. Also note that to avoid unnecessary burdens to notation in the following calculations, we have used for $\psi _1(\zeta )$ and $\varphi _1(\zeta )$ the same symbol that we use everywhere else for $\psi _1(x)$ and $\varphi _1(x)$. This, also, should not cause any confusion as the argument of the function will be later specified in the few cases where the two meanings could be confounded.

Later in this section, we will proceed with the integration of these solutions and with their matching. It is therefore reasonable to postulate the innermost layer width to be

(7.11)\begin{equation} \delta_1= \frac{\gamma \bar{d}_e}{k\rho_s J_0}, \end{equation}

which we will later show to be the case, indeed, although we remind that the identification $\delta _1\to \ell$ formally requires us to a posteriori verify that $\delta _1\ll \rho _s\ll 1$. From now on, we will therefore write $\delta _1$ defined by (7.11) in place of $\ell$ defined by (7.7).

7.3 Integration strategy in the warm-electron regime for $\rho _s^2/\bar {d}_e^2\gg 1$

We can now better establish the integration strategy for this problem, which has been generally outlined at the end of § 3. We have identified three overlapping spatial intervals of the integration domain in which different limits of the eigenvalue equations hold. For practical reasons, these can be identified in terms of the variable $\bar {\zeta }$ introduced in (7.1) and (7.2).

1. (i) The innermost region $x\leq \delta _1$, where $\delta _1={\gamma \bar {d}_e}/(k\rho _s J_0)$, to be yet quantified, and (7.9)- and (7.10) are valid. This integration domain applies for $\tilde {\zeta }\ll 1$, i.e. for $x\ll \rho _s$.

2. (ii) The intermediate region $\delta _1 \ll x\leq \delta _2$ where $\delta _2=\rho _s$ and (7.4) and (7.5) hold. This integration domain applies for $\tilde {\zeta }\sim 1$, i.e. for $\delta _1 \ll x \ll a$.

3. (iii) The outer region valid for $x$ such that $\delta _2\ll x$, whose governing equations are those of ideal MHD. Here, (4.2) is used. This integration domain applies for $\tilde {\zeta }\gg 1$, i.e. for $x \sim a$.

Two matching regions can be therefore identified.

1. (a) In the interval $\delta _1 < x< \delta _2$, the solution ${\chi }({\zeta })$ of (7.10) is matched with solution $\tilde {\chi }(\tilde {\zeta })$ of (7.5) according to $\lim _{{\zeta }\to \infty }{\chi }({\zeta })= \lim _{\tilde {\zeta }\to 0}\tilde {\chi }(\tilde {\zeta })$.

2. (b) In the interval $\delta _2< x< a$, the solution $\tilde {\chi }(\tilde {\zeta })$ of (7.5) (corresponding to $\tilde {\psi }_1(\tilde {\zeta })$ and $\tilde {\varphi }_1(\tilde {\zeta })$ of (7.4)) is matched with the solutions $\psi _{\text {out}}(x)$ and $\varphi _{\text {out}}(x)$ of (4.1) and (4.2) according to $\lim _{\tilde {\zeta }\to \infty }\tilde {\psi }_1(\tilde {\zeta }) = \lim _{{x}\to 0}{\psi _{\text {out}}}(x)$ and $\lim _{\tilde {\zeta }\to \infty }\tilde {\varphi }_1(\tilde {\zeta }) = \lim _{{x}\to 0}{\varphi _{\text {out}}}(x)$.

Notice that the matching conditions also generally depend on the wavelength range, that is, on the value of $\varDelta '$, which rules the solution in the ideal region. This can influence the integration strategy, as outlined in § 5.5.

Before solving (7.1) or (7.3) in both the inner and the outer regions, and then matching the corresponding solutions to obtain the general solution of the problem, it is useful to discuss (see § 7.3.1) a general method that can be used to solve equations having the form of (7.9). This method is based on the integration in the complex plane via representation by means of hypergeometric functions.

7.3.1 Solution of the equation in the innermost interval ($x\ll \delta _1$) via integral representation of hypergeometric functions

Equation (7.9) is a second-order ordinary differential equation that in the complex plane has three singular points at $\zeta = \pm i$ and $\zeta = \infty$. Therefore, employing some transformations, it can be written in the standard form of a hypergeometric equation. This suggests that we look for the solution of (7.9) and (7.10) in terms of hypergeometric functions. Inspired by the integral representation of hypergeometric functions, we can thus look for a solution of (7.9) of the form

(7.12)\begin{equation} \psi_1 = \frac{1}{2{\rm \pi} {\rm i}} \oint_C F(s) (1 + \zeta^2)^s \,{\rm d}s, \end{equation}

where $C$ is a closed contour in the complex plane which starts at $-\infty$ of the real axis, goes around the singularities of the function $F(s)$ (it will turn out that these singularities are given by an infinite number of simple poles of $F(s)$) and then returns back to $-\infty$. This representation will allow us to employ the residue theorem which states that $\psi _1$ equals the sum of residues enclosed by the contour $C$. The second derivative of $\psi _1$ becomes

(7.13)\begin{equation} \psi_1'' = \oint_C 2s (2s - 1) F(s) (1 + \zeta^2)^{s - 1} \,{\rm d}s - \oint_C 4s(s - 1) F(s) (1 + \zeta^2)^{s - 2} \,{\rm d}s. \end{equation}

Shifting the operator in the second term by taking $s' = s + 1$, (7.13) reads

(7.14)\begin{equation} \psi_1'' = \oint_C 2s (2s - 1) F(s) (1 + \zeta^2)^{s - 1} \,{\rm d}s - \oint_C 4s(s + 1) F(s + 1) (1 + \zeta^2)^{s - 1} \,{\rm d}s. \end{equation}

Substituting the previous relation in (7.9), one obtains an expression stating that the integral of $(1+\zeta ^2)^s$ times a function of $s$ equals zero. A sufficient condition for its validity is to set to zero the function of $s$, which reads

(7.15)\begin{equation} [2s (2s - 1) - (\delta_1/\bar{d}_e)^2] F(s) - 4s (s + 1) F(s + 1) = 0. \end{equation}

To get an insight on the form of $F(s)$, let us write $F(s+k)$ in terms of $F(s)$ by using (7.15). After a little algebra, we obtain

(7.16)\begin{equation} F(s + 1) = \frac{\left(s - \dfrac{\alpha}{2}\right) \left(s - \dfrac{\beta}{2}\right)}{s (s + 1)} F(s), \end{equation}

where $\alpha =\frac {1}{2} - \nu$, $\beta = \frac {1}{2} + \nu$ and $\nu = ({1}/{4} + (\delta _1/\bar {d}_e)^2)^{{1}/{2}}$. Then,

(7.17)\begin{equation} F(s + j) = \frac{\left(s - \dfrac{\alpha}{2}\right)_j \left(\dfrac{s - \beta}{2}\right)_j}{(s)_j (s + 1)_j} F(s), \end{equation}

where

(7.18)\begin{equation} (q)_j \equiv q (q + 1) \cdots (q + j - 1) =\frac{\varGamma(q + j)}{\varGamma(q)} \end{equation}

is the so-called ‘Pochhammer symbol’ defined in terms of the $\varGamma$-function (see Abramowitz & Stegun Reference Abramowitz and Stegun1964, § 6.1.22). Motivated by the form of (7.17) one can look for a function $F(s)$ expressed as a combination of Gamma functions, that is,

(7.19)\begin{equation} F(s) = \frac{\varGamma(s + a) \varGamma(s + b)}{\varGamma(s + c) \varGamma(s + d)}. \end{equation}

If, by comparison with (7.17), we choose $c = 0$ and $d = 1$, (7.15) becomes

(7.20)\begin{equation} \left(a + b + \frac{1}{2}\right) s + ab + \left(\frac{\delta_1}{2 \bar{d}_e}\right)^2 = 0. \end{equation}

The previous equation should be satisfied for any value of $s$. Therefore, one obtains

(7.21)\begin{equation} a + b + \frac{1}{2} = 0, \quad ab + \left(\frac{\delta_1}{2 \bar{d}_e}\right)^2 = 0, \end{equation}

which has the following solutions:

(7.22)\begin{equation} a ={-}\frac{\alpha}{2}, \quad b ={-}\frac{\beta}{2}, \quad \alpha = \frac{1}{2} - \nu, \quad \beta = \frac{1}{2} + \nu, \quad \nu = \left( \frac{1}{4} + \frac{\delta_1^2}{\bar{d}_e^2} \right)^{{1}/{2}}. \end{equation}

Here, $a$ and $b$ are real numbers because $(\delta _1/\bar {d}_e)^2$ is also a real number, whence it is evident that the singularities of $F(s)$ are points on the real axis. Using (7.22), $F(s)$ reads

(7.23)\begin{equation} F(s) = \frac{\varGamma\left(s - \dfrac{\alpha}{2}\right) \varGamma\left(s - \dfrac{\beta}{2}\right)}{\varGamma(s) \varGamma(s + 1).} \end{equation}

One sees that $F(s)$ has two series of simple poles at

(7.24)\begin{equation} s = \frac{\alpha}{2} - m, \quad s = \frac{\beta}{2} - m, \quad m = 0, 1, 2, 3,\ldots. \end{equation}

It is also clear that the previous poles are associated with two series of residues. The first one is obtained by taking the limit $s = {\alpha }/{2} - m$ in the evaluation of the residue of $F(s) (1 + \zeta ^2)^s$,

(7.25)\begin{align} \psi_\alpha &= \text{Res}_{\alpha} \equiv \sum_{m = 0}^{\infty} \lim_{s \to {\alpha}/{2} - m} \left[ (s + m) \varGamma\left(s - \frac{\alpha}{2}\right) \right] \frac{\varGamma(-\nu - m) (1 + \zeta^2)^{{\alpha}/{2} - m}}{\varGamma\left(\dfrac{5}{4} - \dfrac{\nu}{2} - m\right) \varGamma\left(\dfrac{1}{4} - \dfrac{\nu}{2} - m\right)} \nonumber\\ &= \sum_{m = 0}^{\infty} \frac{({-}1)^{m}}{m!} \frac{\varGamma(-\nu - m)}{\varGamma\left(\dfrac{5}{4} - \dfrac{\nu}{2} - m\right) \varGamma\left(\dfrac{1}{4} - \dfrac{\nu}{2} - m\right)} (1 + \zeta^2)^{{\alpha}/{2} - m}, \end{align}

where the limit in the previous equation is obtained using the property $\varGamma (s + 1) = s\varGamma (s)$, which leads to

(7.26)\begin{align} &\lim_{s\to -m} (s+m)\varGamma(s) = \lim_{s\to -m} (s+m)\frac{\varGamma(s + 1)}{s} \nonumber\\ &\qquad = \lim_{s\to -m} (s+m)\frac{\varGamma(s + m + 1)}{s (s + 1) \cdots. (s + m - 1)(s+m)} = \frac{({-}1)^{m}}{m!}. \end{align}

The second series of the residues can be calculated taking the limit $s = {\beta }/{2} - m$. Then,

(7.27)\begin{equation} \psi_\beta = \text{Res}_\beta = \sum_{m = 0}^{\infty} \frac{({-}1)^{m}}{m!} \frac{\varGamma(\nu - m)}{\varGamma\left(\dfrac{5}{4} + \dfrac{\nu}{2} - m\right) \varGamma\left(\dfrac{1}{4} + \dfrac{\nu}{2} - m\right)} (1 + \zeta^2)^{{\beta}/{2} - m}. \end{equation}

Therefore, using the residue theorem ($\oint _C f(z) \,{\rm d}z = 2 {\rm \pi}i \sum _{z \in C} \text {Res}[ f(z)]$), and substituting in (7.12), $\psi _1$ becomes

(7.28)\begin{equation} \psi_1 = \psi_\alpha + \psi_\beta. \end{equation}

We recall that this solution is valid as long as $x\ll \delta _2$, i.e. $x \ll \rho _s$.

The same result could be found by using the Frobenius method (see, e.g. Bender & Orszag Reference Bender and Orszag1978, § 3.3) by assuming a power series solution of the form $\psi _1 = \sum _{m = 0}^{\infty } g_m (1 + \zeta ^2)^{\lambda - m}$: after substituting $\psi _1$ in (7.9) and using $g_0\neq 0$, the indicial equation leads to $\lambda _1 = {\alpha }/{2}$ and $\lambda _2 = {\beta }/{2}$, where $\alpha$ and $\beta$ are given by (7.22). The coefficients $g_m$ for both $\lambda _{1}$ and $\lambda _{2}$ can be found by equating the coefficients of the powers of $\zeta$ to zero. In this way, one obtains again (7.25)–(7.28).

A discussion about both the convergence and the independence of solutions (7.25) and (7.27) is in Appendix E.

7.4 Solution of the auxiliary equation for $\rho _s^2/d_e^2\gg 1$: large-$\varDelta '$ limit

We now solve the auxiliary equation (7.10) in the large-$\varDelta '$ limit. Its detailed analytical investigation becomes essential in this wavelength limit since, as we will see in § 9.3, the ‘standard’ heuristic derivation here (surprisingly) fails to recover the correct scaling laws.

In this wavelength limit, the solution of the auxiliary equation in the innermost region can be matched directly to the solution obtained in the outermost non-ideal region by setting $\chi _\infty =0$, as discussed in § 5.5. From this matching, we will obtain the scaling laws for both the growth rate and the reconnecting layer width.

7.4.1 Solution in the innermost layer $x\leq \delta _1$

We can solve (7.10) by following the same line of thought used to solve (7.9), as described in § 7.3.1: we thus assume a solution of the form

(7.29)\begin{equation} \chi(\zeta) = \frac{1}{2{\rm \pi} i} \oint_C G(s) (1 + \zeta^2)^s \,{\rm d}s. \end{equation}

After calculating the first and second derivative of (7.29) with respect to $\zeta$, and substituting them into (7.10), one obtains the condition

(7.30)\begin{equation} [2s(2s - 1) - (\delta_1/\bar{d}_e)^2 ] G(s) - 4(s + 1)^2 G(s + 1) = 0. \end{equation}

By noticing the similarity of the previous relation with (7.15), we look for a $G(s)$ of the form

(7.31)\begin{equation} G(s) = \frac{\varGamma(s + a) \varGamma(s + b)}{\varGamma^2(s + 1)} = \frac{F(s)}{s}. \end{equation}

We recall that the previous equation has been written by assuming $c = d = 1$ in (7.19), and this allows us to eliminate the $(s+1)^2$ coefficient in (7.30). In (7.31), we have also used $\varGamma ^2(s+1) = s \varGamma (s)\varGamma (s+1)$. This choice leads us again to (7.21) and, therefore, to the same values for $a$ and $b$ which are given by (7.22), meaning that both (7.19) and (7.31) have the same infinite series of poles, as can be expected from the definition of $\chi$ expressed by (7.29). Therefore, the general solution of (7.10) reads

(7.32)\begin{equation} \chi_{\text{in}}(\zeta) = \chi_\alpha(\zeta) + \chi_\beta(\zeta), \end{equation}

where

(7.33)\begin{gather} \chi_\alpha(\zeta) = \sum_{m=0}^{\infty} \frac{({-}1)^m}{m!} \frac{\varGamma(-\nu - m)}{\varGamma^2\left(\dfrac{5}{4} - \dfrac{\nu}{2} - m\right)} (1 + \zeta^2)^{{\alpha}/{2} - m}, \end{gather}

(7.34)\begin{gather} \chi_\beta(\zeta) =\sum_{m = 0}^{\infty} \frac{({-}1)^{m}}{m!} \frac{\varGamma(\nu - m)}{\varGamma^2\left(\dfrac{5}{4} + \dfrac{\nu}{2} - m\right)} (1 + \zeta^2)^{{\beta}/{2} - m}. \end{gather}

We close the discussion in this section by noticing that the leading terms in (7.32) are $(1 + \zeta ^2)^{{\alpha }/{2}}$ and $(1 + \zeta ^2)^{{\beta }/{2}}$. Therefore, we can look for an approximate asymptotic solution of the form

(7.35)\begin{equation} \chi_{\text{in}}(\zeta) \approx \chi_{\text{in}}^{\text{appr}}(\zeta) = A (1 + \zeta^2)^{{\alpha}/{2}} + B (1 + \zeta^2)^{{\beta}/{2}} + C, \end{equation}

where $A$, $B$ and $C$ are constants. Their values will be estimated by matching the solutions in the different regions and by applying the constraints at $\zeta =0$, which we are going to obtain in § 7.4.2. Note that for limit $\zeta \to \infty$, (7.35) can be further approximated to

(7.36)\begin{equation} \chi_{\text{in}}(\zeta) \approx \chi_{\text{in}}^{\text{appr}}(\zeta) = A \zeta^\alpha + B \zeta^\beta + C. \end{equation}

7.4.2 Constraints at the origin

Any solution of the auxiliary equation must satisfy the constraints below that follow from the definition of $\chi (\zeta )$ in (5.11) and from (7.9) at $\zeta = 0$:

(7.37)\begin{equation} \psi_1''(0) = \frac{\delta_1^2}{\bar{d}_e^2} \psi_1(0) = \chi''(0) ={-} \frac{\delta_1^2}{\bar{d}_e^2} \chi(0), \end{equation}

(7.38)\begin{equation} \psi_1^{iv}(0) = \frac{\delta_1^2}{\bar{d}_e^2} \psi_1''(0) - 2 \frac{\delta_1^2}{\bar{d}_e^2} \psi_1(0) ={-} 2 \left(1 - \frac{1}{2}\frac{\delta_1^2}{\bar{d}_e^2}\right) \psi_1''(0) ={-}2\frac{\delta_1^2}{\bar{d}_e^2} \left(1 - \frac{1}{2}\frac{\delta_1^2}{\bar{d}_e^2} \right) \psi_1(0), \end{equation}

(7.39)\begin{equation} \chi^{iv}(0) = 3 \psi_1^{iv}(0) ={-}6\left(1 - \frac{1}{2}\frac{\delta_1^2}{\bar{d}_e^2}\right)\psi_1''(0) ={-}6\left(1 - \frac{1 }{2}\frac{\delta_1^2}{\bar{d}_e^2}\right) \chi''(0) = 6\frac{\delta_1^2}{\bar{d}_e^2} \left(1 - \frac{1}{2}\frac{\delta_1^2}{\bar{d}_e^2}\right) \chi(0). \end{equation}

These expressions will be later used in the asymptotic matching, so as to determine the free coefficients in the expression of $\chi ^{\text {aprr}}_{\text {in}}$ in (7.35).

7.4.3 Solution in the intermediate region, $\delta _1 \ll x \leq \delta _2$

In this region, we look for a solution $\tilde {\chi }(\tilde {\zeta })$ of (7.3). Being interested in the large-$\varDelta '$ limit, we can pose $\tilde {\chi }_\infty = 0$. Equation (7.3) then becomes the homogeneous equation

(7.40)\begin{equation} \tilde{\chi}'' - \left(1 + \frac{\tilde{\mathcal{G} }^2}{\tilde{\zeta}^2}\right) \tilde{\chi} = 0, \end{equation}

where ‘$'$’ refers to the derivative with respect to $\tilde {\zeta }=x/\rho _s$. This equation belongs to the family of modified Bessel equations (see Abramowitz & Stegun Reference Abramowitz and Stegun1964, § 9.6.1). In fact, it can be transformed into the standard modified Bessel equation via the change of variable $\chi \tilde {\zeta } = \tilde {\zeta }^{{1}/{2}}u(\tilde {\zeta })$ that leads to

(7.41)\begin{equation} \tilde{\zeta}^2 u'' + \tilde{\zeta} u' - \left(\tilde{\zeta}^2 + \nu^2\right) u = 0, \end{equation}

where $\nu$ is given by the last of (7.22). The only solution of the previous equation that is bounded when $\tilde {\zeta } \to \infty$ (or equivalently, when $x \to \infty$) is the modified Bessel function of the second kind, that is, $u(\tilde {\zeta }) = \mathcal {K}_\nu (\tilde {\zeta })$. Therefore, (7.40) has the solution

(7.42)\begin{equation} \tilde{\chi}_{\text{mid}}(\tilde{\zeta}) = \tilde{\zeta}^{{1}/{2}} \mathcal{K}_\nu(\tilde{\zeta}), \end{equation}

the label ‘$_{\text {mid}}$’ standing for the intermediate region domain $\delta _1 \ll x \leq \delta _2$ where it is defined. Its asymptotic behaviour as $\tilde {\zeta } \to 0$ is given by

(7.43)\begin{align} \tilde{\chi}^{\text{appr}}_{\text{mid}} & \equiv \lim_{\tilde{\zeta}\ll 1} \tilde{\chi}_{\text{out}}(\tilde{\zeta}) = \lim_{\tilde{\zeta} \ll 1} \tilde{\zeta}^{{1}/{2}} \mathcal{K}_\nu(\tilde{\zeta}) = \lim_{\tilde{\zeta} \ll 1} \tilde{\zeta}^{{1}/{2}} \left(\frac{\rm \pi}{2}\right) \frac{\mathcal{I}_{-\nu}(\tilde{\zeta}) - \mathcal{I}_{\nu}(\tilde{\zeta})}{\sin({\rm \pi} \nu)}\nonumber\\ &\sim \frac{2^\nu}{\varGamma(1 - \nu)} \tilde{\zeta}^\alpha - \frac{2^{-\nu}}{\varGamma(1 + \nu)} \tilde{\zeta}^\beta, \end{align}

where we used the relation in § 9.6.2 of Abramowitz & Stegun (Reference Abramowitz and Stegun1964). Here, $\alpha$ and $\beta$ are once more given by (7.22), and $\mathcal {I}_{\nu }(\tilde {\zeta })$ is the modified Bessel function of the first kind,

(7.44)\begin{equation} \mathcal{I}_{\nu}(\tilde{\zeta})\equiv \left(\frac{\tilde{\zeta}}{2}\right)^\nu \sum_{m = 0}^{\infty} \frac{\displaystyle{\left(\frac{\tilde{\zeta}^2}{4}\right)^m}}{m! \varGamma(\nu + m + 1)}. \end{equation}

7.4.4 Matching of solutions: scaling laws and eigenfunctions

The scaling laws for both the growth rate $\gamma$ and the width of the innermost layer, $\delta _1$, are obtained from matching the solutions in the innermost and intermediate regions of the domain. This formally closes the eigenvalue problem.

To accomplish this, we start by estimating the second and the fourth derivative of $\chi _{\text {in}}^{\text {appr}}(\zeta )$ at $\zeta = 0$. Then we will apply the constraints at the origin written in § 7.4.2 and the condition $\lim _{\zeta \to \infty }\chi (\zeta )=\lim _{\tilde {\zeta }\to 0}\tilde {\chi }(\tilde {\zeta })$, that is, $\lim _{\zeta \to \infty }\chi ^{\text {appr}}_{\text {in}}(\zeta )=\lim _{\tilde {\zeta }\to 0}\tilde {\chi }^{\text {appr}}_{\text {mid}}(\tilde {\zeta })$, so as to obtain the equations that will allow us to find $A$, $B$ and $C$, and, finally, the scaling laws.

First of all, from the matching of (7.43) with (7.36), we obtain

(7.45)\begin{equation} A = \frac{2^\nu}{\varGamma(1 - \nu)} \left(\frac{\delta_1}{\rho_s}\right)^\alpha, \quad B ={-}\frac{2^{-\nu}}{\varGamma(1 + \nu)} \left(\frac{\delta_1}{\rho_s}\right)^\beta. \end{equation}

Taking the second and the fourth derivatives of $\chi _{\text {in}}^{\text {appr}}$ (we here omit superscript ‘${\text {appr}}$’ to simplify the notation), one finds

(7.46)\begin{gather} \chi_{\text{in}}(0) = A + B + C, \end{gather}

(7.47)\begin{gather} \chi_{\text{in}}''(0) = A \alpha + B \beta, \end{gather}

(7.48)\begin{gather} \chi_{\text{in}}^{iv}(0) = 6 A \alpha \left(\frac{\alpha}{2} - 1\right) + 6 B \beta \left(\frac{\beta}{2} - 1\right). \end{gather}

Then, using the constraint on the second derivative given by (7.37) and (7.47),

(7.49)\begin{equation} C ={-} A \left(1 + {\alpha}\left(\frac{\bar{d}_e}{\delta_1}\right)^2\right) - B \left(1 + {\beta}\left(\frac{\bar{d}_e}{\delta_1}\right)^2\right). \end{equation}

The constraints on the fourth derivative in (7.39) and on (7.48) instead lead to

(7.50)\begin{equation} A \alpha \left(\alpha - \left(\frac{\delta_1}{\bar{d}_e}\right)^2\right) + B \beta \left(\beta - \left(\frac{\delta_1}{\bar{d}_e}\right)^2 \right) = 0. \end{equation}

Combining the previous equation with (7.45), one obtains

(7.51)\begin{equation} -\frac{A}{B} = \frac{\beta \left(\beta - \left(\delta_1/\bar{d}_e\right)^2\right)}{\alpha \left(\alpha - \left(\delta_1/\bar{d}_e\right)^2\right)} = 2^{2\nu} \left(\frac{\delta_1}{\rho_s}\right)^{{-}2\nu} \frac{\varGamma(1 + \nu)}{ {\varGamma(1 - \nu)}}, \end{equation}

where we used $\alpha - \beta = - 2\nu$. For small values of $(\delta _1/\bar {d}_e)$, (7.22) gives $\nu \approx 1/2$, $\alpha \approx - (\delta _1/\bar {d}_e)^2$ and $\beta \approx 1$; therefore, (7.51) becomes

(7.52)\begin{equation} \frac{1}{2}\left(\frac{\bar{d}_e}{\delta_1}\right)^4 \left(1- \left(\frac{\delta_1}{\bar{d}_e}\right)^2\right) = 2 \frac{\varGamma\left(\dfrac{3}{2}\right)}{\varGamma\left(\dfrac{1}{2}\right)} \left(\frac{\rho_s}{\delta_1}\right). \end{equation}

Using $\varGamma (3/ 2)=\varGamma (1/ 2)/ 2 = \sqrt {{\rm \pi} }/ 2$ and $(\delta _1/\bar {d}_e)^2\ll 1$ (to be verified a posteriori), the equation above implies

(7.53)\begin{equation} \delta_1 \sim \rho_s^{-{1}/{3}} \bar{d}_e^{{4}/{3}} \end{equation}

whence, using the definition in (7.11), one obtains

(7.54)\begin{equation} \gamma \sim k\rho_s^{{2}/{3}} \bar{d}_e^{{4}/{3}}. \end{equation}

Specialising (7.53) and (7.54) to the purely resistive regime, $\bar {d}_e=(S^{-1}/\gamma )^{1/2}$, one recovers the result of Pegoraro & Schep (Reference Pegoraro and Schep1986) for the growth rate (cf. (33) therein and also (20) in Pegoraro et al. Reference Pegoraro, Porcelli and Schep1989) and of Zocco & Schekochihin (Reference Zocco and Schekochihin2011) for $\delta _1$ (cf. $\delta _\eta$ in (B100), therein),

(7.55)\begin{equation} \gamma \sim k^{{6}/{7}} \rho_s^{{4}/{7}} S^{-{1}/{7}},\quad \delta_1 \sim k^{-{4}/{7}} \rho_s^{-{5/}{7}} S^{-{4}/{7}}. \end{equation}

Specialising the same equations to the purely collisionless regime, $\bar {d}_e=d_e$, one recovers instead the result of Porcelli (Reference Porcelli1991) for both for the growth rate ((8), therein) and for the innermost layer width (quantity $\sigma$, defined therein above (7); see also Bhattacharjee, Germaschewski & Ng Reference Bhattacharjee, Germaschewski and Ng2005),

(7.56)\begin{equation} \gamma \sim k \rho_s^{{2}/{3}} d_e^{{1}/{3}},\quad \delta_1 \sim \rho_s^{-{1}/{3}} d_e^{{4}/{3}}. \end{equation}

In both the purely resistive and the collisionless case, the conditions $\delta _1^2\ll \bar {d}_e^2$ and $\delta _1\ll \delta _2\sim \rho _s$ we had previously heuristically assumed, are a posteriori verified.

Finally, by looking at the definition of $\chi$, it is obvious that $\psi _1 \approx \chi _{\text {in}}^{\text {appr}}$ when $x \ll \delta _1$, since $x\psi '_1 \sim {x}\psi _1/{\delta _1} \ll 1$. Therefore, substituting the values of $\nu$, $\alpha$, $\beta$, $\mathcal {G}$ and $\delta _1$ in (7.35), one obtains for the inner layer

(7.57)$$\begin{gather} \lim_{x\ll \delta_1}\psi_{1} \simeq \frac{2^{1/2}}{\varGamma(1/2)} \left( \frac{\delta_1}{\rho_s} \right)^{- {\delta_1^2}/{\bar{d}_s^2}} \left( 1 + \frac{x^2}{\delta_1^2} \right)^{- {\delta_1^2}/{\bar{d}_s^2}} - \frac{2^{{-}1/2}}{\varGamma(3/2)} \displaystyle{\left( \frac{\delta_1}{\rho_s} \right) \left( 1 + \frac{x^2}{\delta_1^2} \right)^{{1}/{2}} } \nonumber\\ + \left( \frac{\delta_1}{\rho_s} \right) \left( 1 + \frac{\bar{d}_s^2}{\delta_1^2}\right). \end{gather}$$

Integrating (5.12) and using (7.35), one finds for $x \ll \delta _1$,

(7.58)$$\begin{gather} \lim_{x\ll \delta_1} \varphi_{1} \simeq \sum_{m = 0}^{\infty} \left[ \frac{2^{1 /2}}{\varGamma({1}/{2})} \left( \frac{\delta_1}{\rho_s} \right)^{- {\delta_1^2}/{\bar{d}_s^2} } \left( \begin{matrix} -{\delta_1^2}/{\bar{d}_s^2} \\ m \end{matrix} \right)_b - \frac{2^{{-}1 /2}}{\varGamma({3}/{2})} \left( \frac{\delta_1}{\rho_s} \right) \left( \begin{matrix} {-1}/{2} \\ m \end{matrix} \right)_b \right] \left( \frac{x}{\delta_1} \right)^{2\,m + 1} \nonumber\\ + \left( \frac{\delta_1}{\rho_s} \right) \left( 1 + \frac{\bar{d}_s^2}{\delta_1^2} \right) \left( \frac{x}{\delta_1} \right), \end{gather}$$

where

(7.59)\begin{equation} \left( \begin{matrix} r \\ m \end{matrix} \right)_b = \displaystyle{\frac{r (r - 1) \cdots (r - m + 1)}{m!}} = \displaystyle{\frac{\left( r \right)_m}{m!}} \end{equation}

is the usual binomial coefficient, generalised for real values of the argument $r$.

Notice that for ${x} \leq \delta _1$, the dominant term in (7.58) is the one with $m=0$. Therefore, to the lowest order, $\varphi$ can be approximated as $\varphi \sim x$. Figure 8 shows the spatial profile of $\psi _1$, $\chi _{\text {in}}^{\text {appr}}$ and $x\psi _1'$ in blue, black and orange colours, respectively, in a sub-interval of the integration domain. These profiles have been computed numerically by solving the eigenvalue problem for $\psi _1$ and $\varphi _1$ with the solver of Betar et al. (Reference Betar, Del Sarto, Ottaviani and Ghizzo2020) and using the definitions in (7.36). We can this way numerically verify that $\chi _{\text {in}}^{\text {aprr}} \approx \psi _1$ within the innermost layer with $|x|\ll \delta _1$, in agreement with the previous discussion.

Figure 8.

Profiles of $\psi _1$ (blue), $\chi _{\text {in}}^{\text {appr}}$ of (7.46) (black-dashed) and $x \psi _1'$ (orange) in the reconnection layer around the neutral line: these profiles are obtained by numerical integration of the complete eigenvalue problem and should be compared to the approximated analytical solutions ((7.57) and (7.58)). The width of reconnection layer, evaluated as it will be specified in § 7, is $\delta = \delta _1 \simeq 10^{-6}$. The physical parameters of the numerical integration performed with the solver of Betar et al. (Reference Betar, Del Sarto, Ottaviani and Ghizzo2020) in the large-$\varDelta '$ limit are shown in the figure.

The eigenfunctions (7.57) and (7.58) agree with the behaviour of the perturbed current density profile provided by Pegoraro & Schep (Reference Pegoraro and Schep1986) in the inner layer for the warm-resistive regime ($\bar {d}_e=(S^{-1}/\gamma )^{1/2}$), which had been already numerically verified by Betar et al. (Reference Betar, Del Sarto, Ottaviani and Ghizzo2020).