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

HAL
open science

# Microscopic scales of linear tearing modes: a tutorial on boundary layer theory for magnetic reconnection

H. Betar, D. Del Sarto, M. Ottaviani, A. Ghizzo

► To cite this version:

H. Betar, D. Del Sarto, M. Ottaviani, A. Ghizzo. Microscopic scales of linear tearing modes: a tutorial on boundary layer theory for magnetic reconnection. Journal of Plasma Physics, 2022, 88 (6), pp.925880601. 10.1017/S0022377822001088. hal-03809727v2

HAL Id: hal-03809727

https://hal.science/hal-03809727v2

Submitted on 27 Jun 2023

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# TUTORIAL

## Microscopic scales of linear tearing modes: a tutorial on boundary layer theory for magnetic reconnection

H. Betar<sup>1,2,†</sup>, D. Del Sarto<sup>2,†</sup>, M. Ottaviani<sup>3</sup> and A. Ghizzo<sup>2</sup>

(Received 31 January 2022; revised 8 October 2022; accepted 12 October 2022)

We revise in detail and in a pedagogical way the analysis of the boundary layer theory of warm tearing modes in slab, reduced magnetohydrodynamics (MHD), when magnetic reconnection is driven by electron inertia and/or resistivity, and ion-sound Larmor radius effects are included. By comparison with the numerical solution of the corresponding eigenvalue problem, we interpret these results by means of a heuristic approach, which in the warm-electron regime, we show to be in general not feasible without knowledge of the scaling of the gradient of the magnetic flux function, differently from what happens in the cold-electron regimes. We put in evidence for a non-trivial relation between the first derivative of the magnetic flux function and of the velocity parallel to the neutral line, evaluated in its proximity, by thus providing insight to the multiple boundary layer analysis that Pegoraro & Schep (Plasma Phys. Control. Fusion, vol. 28, 1986, p. 647) first showed to be required in warm-tearing regimes. In this way, we also suggest and justify a general operational definition of the reconnecting layer width and we discuss the linear appearance of microscopic scales related to the gradients of the eigenfunctions of the tearing modes.

**Key words:** fusion plasma, astrophysical plasmas, plasma instabilities

## Contents

1 Introduction 4

2 Model equations and linear problem: general features 9

2.1 Magnetic reconnection in tearing modes and formation of magnetic islands 10

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

<footnote n="1">¹Laboratoire M2P2, UMR CNRS 7340 – Université Aix-Marseille, F-13451 Marseille, France</footnote>

<footnote n="2">²Institut Jean Lamour, UMR 7198 CNRS – Université de Lorraine, F-54000 Nancy, France</footnote>

<footnote n="3">³CEA, IRFM, F-13108, Saint-Paul-lez-Durance, France</footnote>

<footnote>† Email addresses for correspondence: homam.betar@univ-amu.fr, daniele.del-sarto@univ-lorraine.fr</footnote>
8.2 Operational definition of the reconnecting layer width, $\delta$, and its scalings

8.3 Further micro-scales related to the gradients of the eigenmodes on the neutral line

8.3.1 Asymptotic scalings of $\delta_{\psi}^{(2)}$ and of $\delta_{\psi}^{(4)}$

8.3.2 Asymptotic scalings of $\delta_{\psi}^{(2)}$

8.3.3 Asymptotic scalings of $\delta_{\psi}^{(4)}$

9 Heuristic derivation of the scaling laws of tearing modes

9.1 General hypotheses in the heuristic approach to the scaling estimate

9.2 Heuristic derivation of the scaling laws in the cold-electron regime ($\rho_s^2 \ll \bar{d}_e^2$)

9.3 Heuristic derivation of the scaling laws in the warm-electron regime ($\rho_s \gtrsim \bar{d}_e$)

10 An ansatz about the ‘generalisation’ of the heuristic estimates: the role of the velocity gradient in the non-ideal region

10.1 Velocity gradient in the innermost non-ideal region

10.2 A heuristic generalisation of the definition of the scale $l_c$

10.3 Numerical evaluation of $\Delta'$ and $\Delta'_{vy}$

10.4 Role of $\Delta'_{vy}$ in heuristic-type estimates

10.5 Significance of the inverse spatial scale $\Delta'_{vy}$: coherence with boundary layer calculations and comparison with previous work

11 Summary and conclusions

A Appendix A. Derivation of the model equations for tearing modes in slab RMHD

B Appendix B. About alternative definitions of the reconnection rate

C Appendix C. Integration of the boundary layer equation via Fourier transformation

C.1 A brief historical review on the boundary layer approach to tearing-type equations in the Fourier space

C.2 Comparison of the analysis of §§ 4–6 with the Fourier approach

D Appendix D. An example of renormalisation of a differential equation

E Appendix E. Convergence and independence of the solutions of the hypergeometric form of the inner equation

E.1 Convergence of $\psi_{\alpha}$ and $\psi_{\beta}$

E.2 Linear independence of $\psi_{\alpha}$ and $\psi_{\beta}$
F Appendix F. Heuristic estimation of the scaling laws by making use of the inverse gradient scale $l_c$ and of $\Delta'_{vy}$ 94

References 98

1. Introduction

Magnetic reconnection is one of the chief processes of conversion between electromagnetic energy and particle kinetic energy in magnetised plasmas (see, e.g. Biskamp 2000; Yamada, Kulsrud & Ji 2010). At the basis of extreme energy releasing phenomena that naturally occur in space plasmas, such as solar flares, coronal mass ejections or magnetic substorms, magnetic reconnection is of fundamental importance also in magnetically confined thermonuclear fusion plasma experiments: in fusion devices like tokamaks, it can cause disruptions (see, e.g. Wesson 1990), that is, the sudden loss of the magnetic confinement. This is usually due to magnetic island formation via tearing-like instabilities, and can affect the transport of matter and energy. Among the different scenarios, which have been devised since the first formulations of the ‘magnetic reconnection’ concept in the works by Giovanelli (1946), Hoyle (1949) and Dungey (1950), tearing-type modes (Furth, Killeen & Rosenbluth 1963) are the prototypical example of spontaneous reconnecting instabilities, that is, of magnetic perturbations that grow exponentially in time by inducing magnetic reconnection. Their ‘explosive’ behaviour in time make them the chief candidates for the explanation of several fast energy releasing events, which, both in Nature and in laboratory plasmas, are attributed to magnetic reconnection processes. Their theoretical modelling is however not trivial, since the mathematics required to solve the corresponding eigenvalue problem is complicated by some technical features that make its integration typically more complex than that of most linear instabilities. Notably, due to the multi-scale nature of the linear problem in which linear differential operators intervene, a boundary layer integration procedure is required.

Since Furth’s remark (Furth 1963, 1964) of the relevance of electron inertia in allowing magnetic reconnection and since the early model by Furth (1964) and Coppi (1964c,a) for inertia-driven reconnecting instabilities, kinetic scale effects have attracted increasing attention in the attempt to model the rapid magnetic reconnection phenomena observed in low collisionality plasmas. After the seminal paper by Furth et al. (1963), in which the theory of resistive tearing modes was first formulated, several attempts have been made to include collisionless and/or kinetic physics in the linear tearing mode theory. Even if some works exist in which a full kinetic or girokinetic analytical treatment has been considered (Coppi, Laval & Pellat 1966; Laval, Pellat & Vuillemin 1966; Hazeltine & Ross 1975; Drake & Lee 1977; Cowley, Kulsrud & Hahm 1986; Daughton 1999; Daughton & Karimabadi 2005; Zocco & Schekochihin 2011; Connor, Hastie & Zocco 2012b), the complexity of the boundary layer analysis and of the heuristic assumptions on the ordering of the different microscopic scales makes the identification of the different reconnection regimes a quite difficult task, especially when two or more non-ideal magnetohydrodynamic (MHD) parameters enter in the reconnection rate. This is true even in a relatively simple fluid description extended to include non-ideal MHD effects.

For this reason, heuristic approaches have been developed to tackle the boundary layer problem with a simplified dimensional-like analysis, based on estimations about the characteristic gradients and about the balance of the terms of the linear equations. These techniques have been successfully used to recover the asymptotic scalings of the growth rate and of the reconnection layer width in the purely resistive and purely inertial regimes of reduced MHD (RMHD) instabilities. They are also often presented in textbooks on
magnetic reconnection as a ‘short-cut’ procedure to obtain these results without carrying out the full boundary layer integration of the eigenmodes. Providing an estimate of the characteristic temporal and spatial scales of spontaneous reconnecting instabilities, and comparing them with the values inferred from experimental measures of reconnection events occurring in the laboratory or in Nature, are indeed among the elements of principal interest, in this context. Heuristic methods based on dimensional analysis have also proven to work for both small and large values of the $\Delta'$ instability parameter (Ottaviani & Porcelli 1995) and for the fastest growing mode in a large aspect ratio current sheet (Bhattacharjee et al. 2009; Comisso et al. 2013; Del Sarto et al. 2016). Moreover, they are generally used to obtain insight on the physics of the problem and on the interpretation of some non-trivial results of the boundary layer analysis (see, e.g. Drake & Lee 1977; Cowley et al. 1986). However, as we are going to discuss in this work, heuristic methods ‘fail’ to obtain the scalings first computed by Pegoraro & Schep (1986) and Porcelli (1991) with a boundary layer approach, when electron temperature effects are included, unless some further careful assumptions are made, which are related to the boundary layer decomposition of the spatial domain: whether it is possible to develop heuristic arguments which allow one to obtain these results without relying on further information from boundary layer analysis (or from numerical integration of the linear problem) is, at the moment, an open question.

One of the main purposes of this work is to revise in a tutorial way the linear analysis of tearing-type modes in these regimes, by discussing their analytical solutions in the coordinate space and by showing in this way the advantages and limitations of heuristic-type derivations. To the best of our knowledge, indeed, details of the analytical theory for this double-boundary layer approach have been seldom discussed in the existing literature. When this has occurred, the analysis was based on the solution of the eigenvalue problem in the Fourier space (Pegoraro & Schep 1986), or by using a different linear model based on a reduced fluid-kinetic approach (Zocco & Schekochihin 2011) and by taking a somewhat different analytical approach based on perturbation methods to solve the boundary layer equations that here we will tackle by direct integration instead.<fnref n="1" /> It should be nevertheless recalled that other analytical works on this subject exist. In these, the boundary layer analysis of tearing modes is however complicated by further ingredients related, for example, to the geometry of the magnetic equilibrium profile (see, e.g. Militello et al. 2004; Connor et al. 2012a; Zocco, Helander & Weitzner 2020) or to the inclusion of full kinetic effects (see the references previously cited in this regard). Moreover, also in the few works that addressed these subjects with different approaches with respect to the one we develop here, most details of the complex analysis involved in the calculations were often not reported and are difficult to track.

To simplify the analysis and to show in a pedagogical way the essential features of the double boundary layer separation, we focus here on the limit of cold ions, which allows us to take a fluid closure for the latter, differently from the semi-kinetic models considered in practically all previous works that already treated this problem (see, e.g. Cowley et al. 1986; Pegoraro & Schep 1986; Pegoraro, Porcelli & Schep 1989; Porcelli 1991; Zocco & Schekochihin 2011). In particular, we rely on the equations of the two-fluid reduced-MHD model extended to include ion-sound Larmor radius effects, for which different derivations exist (Zank & Matthaeus 1992; Schep, Pegoraro & Kuvshinov 1994; Bergmans 2001; Del Sarto, Califano & Pegoraro 2006; Bian & Tsiklauri 2009), which are based on two different

<footnote n="1">¹After acceptance of this manuscript the authors have become aware of another recent work (Granier et al. 2022), in which an integration of the boundary layer equations in the coordinate space has been performed in the warm-electron, small-$\Delta'$ regime of tearing modes.</footnote>
types of ordering between fluctuations of the ion density and of the guide field magnetic component (cf. Appendix A). Several more recent and refined gyrofluid models exist, which are based on the evolution of more than two scalar fields. However, the two-field model we focus on, which has been used in several nonlinear studies of different magnetic reconnection scenarios (see Kleva, Drake & Waelbroeck 1995; Cafaro et al. 1998; Grasso et al. 1999; Bergmans & Schep 2001; Del Sarto, Califano & Pegoraro 2003; Wang et al. 2011 to cite the earliest) or even of turbulence (Milosevich, Morrison & Tassi 2018), contains the whole essential physics of the problem. More specifically, it yields the same linear system to which all other, more refined, cold-ion models converge in the isothermal electron limit.

To accomplish our pedagogical purpose, we then give a step-by-step presentation of the boundary layer integration procedure in the coordinate space, with the aim of keeping it ‘self-contained’, that is, by trying to provide all the analytical tools useful for the purpose (e.g. element of complex analysis), when they result necessary for the algebra. By then comparing the analytical results with those obtained with a numerical eigen-solver (Betar et al. 2020) based on a multi-precision toolbox (Holoborodko 2012), we analyse the spatial behaviour of the corresponding eigenfunctions and we discuss the limitations and delicate points of heuristic estimations that are sometimes used as a ‘quicker’ alternative to full boundary layer calculations. Although the main focus of this article is on the warm-electron regime, we also consider, for comparison, the cold-electron limit, in which a single boundary layer analysis suffices. In particular, we consider in a unified way both the resistive and collisionless regimes, and we separately address the cold- and warm-electron regimes. Each regime is associated with some characteristic non-ideal parameter: the Lundquist number $S^{-1}$, related to resistivity; the electron skin depth $d_e$, related to electron inertia in the collisionless limit; and the ion-sound Larmor radius $\rho_s$, related to the electron temperature.

We thus discuss the results of the boundary layer analysis, which yields the asymptotic scaling (in terms of the non-ideal parameters ruling the reconnection process) of the growth rates and of some characteristic spatial scales associated with the eigenfunctions, as well as the approximated profile of the latter in some regions of the domain of integration. We then discuss the asymptotic estimates of the first derivatives of the current and velocity field on the neutral line. We also provide (to the best of our knowledge, for the first time) a formal quantitative definition of the reconnection layer width, which we argue to be valid in all reconnection regimes and which is given in terms of numerically measurable quantities that are related to the spatial profile of the eigenfunctions.

Then, in interpreting these quantities in terms of a heuristic approach, whose limitations in the warm-electron case we point out, we provide evidence that the introduction of the scale $\rho_s$ makes a further characteristic ‘mixed’ scale appear in the collisionless large wavelength limit, which is smaller than $d_e$ when $\rho_s > d_e$, since it scales like $\sim \rho_s^{-1/3} d_e^{4/3}$. Such a characteristic length, which in previous works (Porcelli 1991) had been already noted, is identifiable as the reconnecting layer width, in agreement both with the operational definition here proposed and with previous works that already recognised it as such (Zocco & Schekochihin 2011).

By discussing the boundary layer results in light of a heuristic dimensional-type approach, we also show a new, non-trivial relation between the first derivative of the magnetic flux function and the first derivative of the parallel velocity, which holds close to the neutral line. This velocity gradient displays characteristic scalings that depend on the wavelength regime. Information about these scalings is a priori not evident, and the analysis we provide on this subject at the end of this work suggests indeed that such information can not be obtained via simple dimensional analysis, if not in the cold-electron
regimes. Nevertheless, introducing this scale length allows, in any reconnection regime, the writing of the scaling laws in a way that results to be perfectly symmetric between the small and large wavelength limit, provided the characteristic instability parameter of the small wavelength limit (i.e. the $\Delta'$ parameter of Furth et al. 1963) is replaced in the large wavelength limit by the gradient of the velocity (which we have here called $\Delta'_{v_y}$, by analogy). This fact, and the correspondence between some characteristic scale length associated with this velocity gradient and scale lengths that in other works based on a kinetic approach have been interpreted in terms of inherently kinetic features (Drake & Lee 1977; Cowley et al. 1986; Ottaviani & Porcelli 1993; Zocco & Schekochihin 2011), suggests that the role of this quantity in the tearing mode linear dynamics deserves future investigation.

The structure of this article is as follows.

(i) In § 2, we discuss the model equations and some general aspects of magnetic reconnection which include: a brief historical review of the notion of magnetic reconnection associated with that of magnetic topology and some general features of tearing modes and of the magnetic structures they induce (§ 2.1); the notion of reconnection rate and its relation with the growth rate of the eigenvalue problem (§ 2.2).

(ii) In § 3, we introduce the notion of boundary layer, we discuss its relevance to tearing mode analysis and we outline the key points of the corresponding integration strategy. Then we introduce the notion of instability parameter ($\Delta'$) and the way it intervenes in the matching of the solutions between the ideal and non-ideal regions of the domain (§ 3.1). We discuss the wavelength regimes of the eigenmode solution, in terms of the amplitude of $\Delta'$, by pointing out similarities and differences between tearing modes in slab and in cylindrical geometry (§ 3.2). In § 3.3, we briefly review previous works in which a boundary layer analysis similar to the one we develop here has been discussed.

(iii) In § 4, we address the boundary layer integration by starting from the ideal region, where the hypotheses of ideal MHD hold, and we detail the integration procedure which allows the analytical evaluation of $\Delta'(k)$ in terms of the instability wavenumber $k$ for a specific equilibrium profile, which exemplifies a larger class of equilibria (§ 4.1).

(iv) In § 5, we introduce the notions of ‘generalised’ resistivity and electron inertia, which allow a unified treatment of the resistive and collisionless case altogether, and we discuss the combination of the two non-ideal parameters. We then discuss the strategy for the identification of the integration layers in the non-ideal region via the general approximations valid in the non-ideal region (§ 5.1) and the criteria of normalisation (§ 5.2). We then introduce the auxiliary function, useful for the integration of the boundary layer problem in some regimes, by discussing its small and large wavelength limits (§ 5.3) and then its normalisation (§ 5.4). Finally, we outline the criteria that may make one prefer to perform the integration of the non-ideal equation by using the auxiliary equations rather than the equations for the two scalar fields (§ 5.5).

(v) In § 6, we solve the boundary layer problem in the cold-electron regimes. First we find the solution in the large-$\Delta'$ limit (§ 6.1) and then in the small-$\Delta'$ limit (§ 6.2).

(vi) In § 7, we address the warm-electron regimes: after discussing the general form of the equations in the non-ideal region (§ 7.1), we identify the two boundary layers of interest in this case (§ 7.2) and then we outline the integration strategy that will be pursued, by relying on the integral representation of hypergeometric functions
(§ 7.3). First we consider the solution in the large-$\Delta'$ (§ 7.4) and then in the small-$\Delta'$ limit (§ 7.5).

(vii) In § 8, we address the problem of identifying the characteristic width of the reconnecting layer by relying on further hypotheses of physical character and by starting from the characteristic scales obtained from the boundary layer integration. After reviewing different notions of the reconnecting layer that have been adopted in the literature over the years (§ 8.1), we propose an operational definition of its width, related to the distance of the local maxima (or minima) of the current density from the neutral line, which can be useful for both experimental and numerical application, and which is shown via numerical integration to provide the same scalings of the inner layer width (§ 8.2). The asymptotic scalings of further characteristic spatial scales associated with the derivatives of the magnetic field and of the velocity profile on the neutral line are then estimated and compared to the scalings of the reconnecting layer width in different regimes (§ 8.3).

(viii) In § 9, we address the problem of the heuristic derivation of the scaling laws of tearing modes via dimensional analysis. After discussing the relevance and usefulness of the approach and after having outlined its general hypotheses (§ 9.1), we apply it to the ‘textbook’ example of the cold-electron regimes, where its efficiency is well established (§ 9.3). Then we show and discuss its ‘failure’ in the warm-electron regimes, where it yields wrong estimates with respect to those obtained from the boundary layer analysis (§ 9.3).

(ix) In § 10, starting from the physical insight brought by the analysis of the logical points which may lead to the wrong dimensional estimates in the warm-electron regime, we introduce a new characteristic scale length associated with the gradient of the velocity component parallel to the neutral line and evaluated close to it: by analogy with the instability parameter, we call it $\Delta'_{v_y}$ (§ 10.1). We then postulate a further estimate for the gradient of the magnetic flux function, which, depending on the reconnection regime, can be related to $\Delta'_{v_y}$ (§ 10.2). After having shown the convergence of a numerical procedure which allows the quantification of $\Delta'_{v_y}$ (§ 10.3), we show its relevance for the heuristic estimates: they can now lead to the correct results – although in the warm-electron regime, the procedure is not ‘closed’, since it is strongly suggested that the scale $\Delta'_{v_y}$ can not be obtained by simple dimensional analysis – and they can be cast in a ‘symmetric form’, in which the scalings in the large-$\Delta'$ limits mirror those of the small-$\Delta'$ limits provided the substitution $\Delta' \leftrightarrow \Delta'_{v_y}$ (§ 10.4). The significance of the inverse scale length $\Delta'_{v_y}$ is further discussed by comparison with results obtained in previous works on boundary layer calculations, in which a characteristic length displaying similar scalings in some regimes had been already noted (§ 10.5).

(x) Conclusions follow in § 11.

(xi) Further technical details are reported in the Appendices. These include: calculations related to the derivation of the model and a discussion of the relevant hypotheses (Appendix A); a brief discussion of alternative definitions found in literature for the notion of ‘reconnection rate’, which are not directly relevant to the problem we consider here (Appendix B); a discussion on the general strategy of integration of the boundary layer problem of tearing modes in the Fourier space, and a comparison of the equations we have integrated with those appearing in previous works available in the literature (Appendix C); a didactical example of normalisation of a differential
equation (Appendix D); the proof of convergence and linear independence of the solution found in § 7.4 (Appendix E); and the detailed discussion of the logical steps, which allow one to generalise the heuristic estimate discussed in § 10.4 in terms of $\Delta'_{v_y}$ (Appendix F).

## 2. Model equations and linear problem: general features

We consider the reduced MHD limit, from now on noted as RMHD (Strauss 1976, 1977; Zank & Matthaeus 1992), and we consider the standard tearing equations in slab geometry ($\partial/\partial z = 0$) for a strongly magnetised plasma whose guide field component is along the $z$ direction. Fluid incompressibility, valid for the bulk plasma at the leading $E \times B$-drift ordering and possibly with inclusion of drift-diamagnetic corrections (see Appendix A), allows us to reduce the number of independent vector components and therefore to consider a limited set of scalar quantities.

Choosing perturbations of the form $\sim f(x)e^{iky+\gamma t}+\text{c.c.}$ we linearise the equations around equilibria with an odd $B_y^0(x)$ component, which vanishes at $x=0$, and with a null in-plane fluid velocity at equilibrium, $U^0(x, y) = 0$. Here we have already used knowledge of the fact that the eigenmodes we are considering have low frequencies, so that $\omega/\gamma \sim 0$. This can be taken here as a verifiable heuristic assumption.

It can be shown that only two scalar fields are necessary. The relevant equations can be then represented by the equation for one of the in-plane components of $B$ and by the equation for one of the in-plane components of the fluid velocity. These can be cast in the form of two coupled equations for the stream functions $\psi$ and $\varphi$. Details about all these features, and bibliographical references as well, are provided in Appendix A.

In particular, $\psi(x, y, t)$ and $b(x, y, t)$ are the normalised scalar functions defining the magnetic field components, $B = \nabla\psi \times e_z + be_z$, whereas $\varphi(x, y, t)$ is the stream function for the in-plane $E \times B$-drift velocity, $U_{\perp} = -\nabla\varphi \times e_z$, expressing the normalised gradient of the electrostatic potential in RMHD. Here, the large guide field ordering allows us to completely neglect the dynamics related to the scalar function $b$.

Labelling with indices '0' and '1' the equilibrium quantities and first-order perturbations, respectively, the linearisation is then performed around equilibria of the form $\psi(x, y, 0) = \psi_0(x)$, even with respect to $x=0$, with $\varphi(x, y, 0) = 0$. Normalised as indicated below and using, for simplicity of notation, the same symbols for the fields defined in the coordinate space and for their Fourier transform with respect to the $y$ variable, the equations we will work with, read:

$$ \gamma[\psi_1 - d_e^2(\psi_1'' - k^2\psi_1)] - ik\varphi_1(\psi_0' - d_e^2\psi_0'') = -ik\rho_s^2(\varphi_1'' - k^2\varphi_1)\psi_0' + S^{-1}(\psi_1'' - k^2\psi_1), \quad (2.1) $$

$$ \gamma[\varphi_1'' - k^2\varphi_1] = ik\psi_0'(\psi_1'' - k^2\psi_1) - ik\psi_1\psi_0'''. \quad (2.2) $$

We can also cast (2.1)–(2.2) in a matrix form (useful for the eigenvalue solution) as

$$ \gamma \begin{pmatrix} \mathcal{F} & 0 \\ 0 & \mathcal{L} \end{pmatrix} \begin{pmatrix} \psi_1 \\ \varphi_1 \end{pmatrix} = \begin{pmatrix} S^{-1}\mathcal{L} & \mathcal{B} \\ \mathcal{A} & 0 \end{pmatrix} \begin{pmatrix} \psi_1 \\ \varphi_1 \end{pmatrix}, \quad (2.3) $$

where we have introduced the differential operators

$$ \mathcal{L} \equiv \frac{\partial^2}{\partial x^2} - k^2, \quad \mathcal{F} \equiv 1 - d_e^2\mathcal{L}, \quad (2.4) $$

$$ \mathcal{A} \equiv ik(\psi_0'\mathcal{L} - \psi_0'''), \quad \mathcal{B} \equiv ik[(\psi_0'' - d_e^2\psi_0''') - \rho_s^2\psi_0'\mathcal{L}], \quad (2.5) $$

with $\mathcal{C}$ being the $\rho_s = 0$ limit of $\mathcal{B}$.
<figure><img src="image_3.png" /></figure>

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 ↔ B and C ↔ 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 ↔ D and B ↔ C.

The non-ideal parameters of the model are normalised to the reference scale length $L_0$. They correspond to the electron skin-depth, $d_e$, to the ion-sound Larmor radius, $\rho_s$, and to the Lundquist number, $S$. The reference length $L_0$ is chosen to be the equilibrium shear length $a$ of $\psi_0(x)$ at $t=0$, whereas times are normalised to the reference Alfvén time $\tau_A$ computed with respect to the in-plane equilibrium magnetic field component, evaluated sufficiently far from the neutral line.

Several derivations exist of (2.1)–(2.2). The parameter $\rho_s$, which violates the Lagrangian invariance of the parallel electron canonical momentum $F \equiv \psi - d_e^2 \nabla^2 \psi$, is related to electron parallel compressibility effects (Grasso et al. 1999). It is also considered to be a finite-Larmor-radius (FLR) effect, since in a strong guide field limit, it can be shown to be due to a component of the non-isotropic electron pressure tensor that enters in the equation of $\psi$ via the diamagnetic drift component of the fluid electron velocity (see Schep et al. (1994) and Appendix A).

## 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.

Non-ideal effects allow a local violation of the topological conservations implied by the ideal Ohm's law (Cowling 1933; Alfvén 1942; Batchelor 1950; Elsasser 1950a,b; Truesdell 1950; Lundquist 1951; Newcomb 1958), 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,

$$ \frac{\partial B}{\partial t} = \nabla \times (U \times B), \qquad (2.6) $$

is equivalent, thanks to a well-known vector identity $(\nabla \times (U \times B) = U(\nabla \cdot B) - B(\nabla \cdot U) + (B \cdot \nabla)U - (U \cdot \nabla)B)$ combined with the continuity equation (Truesdell 1950), to a vector expression that corresponds to the null Lie-derivative of the vector field $B/n$ dragged by the velocity field $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 1933; Alfvén 1942), ‘Woltjer invariants’ (Woltjer 1958) and ‘connection theorem’ (Newcomb 1958) – see (Tur & Yanovsky 1993; Kuvshinov & Schep 1997; Del Sarto et al. 2006) and references therein for a more detailed discussion; see, e.g. Dubrovin, Nivikov & Fomenko (1991, § III.23) for a definition of Lie derivative in tensor notation in the coordinate representation and Schouten (1989, § II.8) for the Lie-derivative of tensor densities of arbitrary ‘weight’ (cf. also Lovelock & Rund (1989, 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 1954). 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 1950; Elsasser 1950b; Truesdell 1950; Axford 1984; Greene 1993). 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 2021). 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 $U$ according to $\varphi \to \varphi - \rho_s^2 \nabla^2 \varphi$ (Pegoraro et al. 2004). 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. 1998; Grasso et al. 2001).

The early notion of magnetic reconnection has been formalised by Dungey (1950, 1953) 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 1946); the second one conjectured that the same may have occurred also in the terrestrial magnetotail (Hoyle 1949). Resistivity is one of the non-ideal effects capable of violating alone the Lagrangian invariance of $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 $\Phi_e$ in generalised Ohm’s law, $\epsilon$ generally indicating the infinitesimally small parameter weighting the non-ideal contribution (we use the same normalisation of (2.1)–(2.2)):

$$ E + U \times B = \Phi_e. \qquad (2.7) $$

Magnetic reconnection at least requires $\nabla \times \Phi_e \neq 0$. Although ideal Ohm’s law at the MHD scale takes the form of (2.7) with $\Phi_e = 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’: $E + \mathbf{u}^e \times \mathbf{B} = 0 \Leftrightarrow E + \mathbf{u}^i \times \mathbf{B} = d_i(\mathbf{J} \times \mathbf{B})/n$. While at MHD scales, where $U \simeq \mathbf{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 $E + \mathbf{u}^e \times \mathbf{B} \simeq 0$ (or, better, as long as $E + \mathbf{u}^i \times \mathbf{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 1957; Sweet 1958), 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 (1990), for applications to the sawtooth crash in tokamaks, and by Bulanov, Pegoraro & Sakharov (1992) 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 (1964) – see also (Vasyliunas 1975) and references therein.

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.* (1963). 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 $B_0(x_s) \cdot 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 $B_0(x_s) \cdot 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_\phi r / (B_\theta R)$, where $B_\phi$ 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 $k(r) \cdot \mathbf{B}(r) = 0$ defines a rational surface for $k = m/(2\pi r)e_\theta + n/(2\pi R)e_\phi$. 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$),
<figure><img src="image_4.png" /></figure>

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 1957; Sweet 1958) or electron inertia (Wesson 1990) 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 1966a; Biskamp & Welter 1980; Kleva et al. 1995), 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 1977; Pritchett, Lee & Drake 1979; Bhattacharjee, Brunel & Tajima 1983). 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 1975; Watson & Craig 1998; Gu et al. 2019). 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 1993).
is at the centre of the island – see figure 2(b) (see also White 1983, 1986) 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 1966a; Biskamp & Welter 1980) and to the ‘coalescence’ of a chain of magnetic islands (Finn & Kaw 1977; Longcope & Strauss 1993) 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. 2004 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. 2016, 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 2001): in this regime, the Lagrangian invariance of the magnetic field is expressed with respect to the electron flow (see Fruchtman 1991; Pegoraro et al. 2004; Del Sarto et al. 2006), 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

$$ \text{X-point:} \quad \begin{cases} \psi \text{ even in } x, & \text{even in } y, \\ \varphi \text{ odd in } x, & \text{odd in } y, \end{cases} \qquad \text{O-point:} \quad \begin{cases} \psi \text{ even in } x, & \text{even in } y, \\ \varphi \text{ odd in } x, & \text{odd in } y. \end{cases} \qquad (2.8) $$

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. 1995; Del Sarto et al. 2011).

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 1979; Syrovatskii 1981; Faganello et al. 2010), and which can also make the latter compete with the Kelvin-Helmholtz instability (Hofmann 1974; Chen & Morrison 1990; Ofman et al. 1991; Chen, Otto & Lee 1997) – cf. also (Einaudi & Rubini 1986; Wang, Lee & Wei 1988; Bettarini et al. 2006; Li & Ma 2012) for a
numerical linear study; the stabilising effect of an in-plane magnetic field component orthogonal to the neutral line (see Nishikawa 1982; Somov & Verneta 1988, 1989); 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 1965; Drake & Lee 1977; Coppi et al. 1979). 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 1984; Cowley & Hastie 1988; Fitzpatrick et al. 1993), but also influence their linear evolution and stability threshold (see, e.g. Mahajan et al. 1978, 1979; Migliuolo, Pegoraro & Porcelli 1991; Grasso, Ottaviani & Porcelli 2001; Yu 2010). 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 1975; Drake & Lee 1977; Gladd et al. 1980; Connor, Cowley & Hastie 1990), which cause the formation of microscopic magnetic islands that can affect turbulent transport (see Doerk et al. 2011) and foster stochastic fluctuations of the magnetic field, which can in turn affect tearing modes (Carreras, Rosenbluth & Hicks 1981) and the magnetic confinement (see, e.g. Firpo 2015). 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 1973; Pritchett, Lee & Drake 1980; Wang et al. 2011); 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 1988; Lütjens & Luciani 2002; Wilson 2012). 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. 2009, 2011; Agullo et al. 2017a,b; Choi 2021). These latter subjects are somewhat related (see, e.g. Brennan et al. 2002) 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 1985; Wang & Bhattacharjee 1997; Fitzpatrick 2008). 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 $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

$$R_{rc} = \frac{1}{\Phi_{\Omega}(B)} \frac{d\Phi_{\Omega}(B)}{dt}, \quad \Phi_{\Omega}(B) = \int_{\Omega} B \cdot dS, \qquad (2.9)$$

where $\Omega$ 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 $\Omega = \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 (2000, 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 1986, 2004; Porcelli, Boucher & Rosenbluth 1996; Boozer 2012; Jardin, Krebs & Ferraro 2020); a theoretical model of the fast release of energy during solar flares and coronal mass ejections (see, e.g. Shibata 1998; Cassak & Shay 2012; Aschwanden, Xu & Jing 2014; Aschwanden et al. 2016; Janvier 2017; Aschwanden et al. 2019; Aschwanden 2020); the mechanism of energy release at the basis of X-rays and γ-ray emissions (gamma-ray bursts) in pulsars, magnetars and nebulae (see, e.g. Lyutikov 2003; Tavani et al. 2011; Uzdensky 2011); 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 ~ $10^6$ over the distance of just a solar radius, supposedly via reconnection in the turbulent coronal plasma (see, e.g. Leamon et al. 2000; Matthaeus & Velli 2011); the energy injected in the magnetosphere via solar wind-magnetosphere interactions and released during geomagnetic storms (see, e.g. Lakhina & Tsurutani 2016), 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 2016; Nelson 2016; Eastwood et al. 2018; Knipp et al. 2018). 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 1957; Petschek 1964; Syrovatskii 1966a,b; Bratenahl & Yeates 1970; Baum, Bratenahl & White 1971; Parker 1973; Schnack 1978; Park, Monticello & White 1984; Chen et al. 1997; Hesse, Forbes & Birn 2005; Comisso & Bhattacharjee 2016; Grasso et al. 2020). 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

$$ \mathcal{R}_{rc} = \frac{1}{\Phi_S(\mathbf{B})} \int_{\Omega} \left( \frac{\partial \mathbf{B}}{\partial t} - \nabla \times (\mathbf{U} \times \mathbf{B}) \right) \cdot d\mathbf{S} = - \frac{c}{\Phi_S(\mathbf{B})} \int_{\Omega} (\nabla \times \mathbf{F}_{\varepsilon}) \cdot d\mathbf{S}, \quad (2.10) $$

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
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 $\mathbf{B} = \nabla \psi \times \mathbf{e}_z + b\mathbf{e}_z$ of (2.1)-(2.2), we can write $dS = dy dz \mathbf{e}_x$ and therefore

$$ \Phi_{\Omega}(\mathbf{B}) = \int_{\Omega} \mathbf{B} \cdot d\mathbf{S} = \int_{\Omega} B_{x,1} \, dy \, dz = \int_{l_z} \, dz \int_{l_y} \frac{\partial \psi_1}{\partial y} \, dy. \quad (2.11) $$

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

$$ \mathcal{R}_{rc} = \frac{\int_{l_y} dy \, \partial^2 \psi_1(0, y, t) / \partial t \, \partial y}{\int_{l_y} dy \, \partial \psi_1(0, y, t) / \partial y}. \quad (2.12) $$

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)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)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’)

$$ \mathcal{R}_{TM}(k) = \gamma(k), \quad (2.13) $$

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 |\nabla \psi_1|^2 dx dy$ 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

$$ \mathcal{R}_{rc} \approx \frac{1}{2\mathcal{E}_c} \left| \frac{\partial \mathcal{E}_c}{\partial t} \right|, \quad (2.14) $$

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 1983, 1986) 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 d_e^2 |\nabla^2 \psi_1|^2 dx dy$ associated with the equilibrium current sheet (Del Sarto, Califano & Pegoraro 2005; Del Sarto et al. 2006). 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 2017) 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 (2014), Del Sarto et al. (2018), and § IX of Betar et al. (2020) 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. (2011) for an example in the collisionless limit.

### 3. Boundary layer approach for tearing modes: large- and small-$\Delta'$ regimes

In this section, we start by describing the general strategy of integration with the boundary layer approach to solve the eigenvalue problem of (2.1)–(2.2). Then, we discuss the ordering and the instability parameter, which defines the unstable spectrum of wavenumbers, and its consequences on the classification of different reconnection regimes in both slab and toroidal geometries. Finally, we close the section by giving a brief review on previous works using the boundary layer integration in the collisionless regime.

There are no exact analytical solutions available for the general eigenvalue problem of (2.1)–(2.2). Approximated analytical solutions can be obtained by using a boundary layer approach as first shown for the purely resistive tearing by Furth et al. (1963) in the constant-$\psi$ regime and by Coppi et al. (1976), Ara et al. (1978) and Basu & Coppi (1981).
<figure><img src="image_5.png" /><figcaption>FIGURE 3. Growth rates of collisionless tearing modes numerically computed from (2.1)–(2.2) with the solver of Betar et al. (2020), 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. 2016; Del Sarto & Pegoraro 2017; Franci et al. 2022). 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.</figcaption></figure>

in the internal-kink regime. When kinetic-like effects are included, it is instead easier to solve the boundary layer problem if calculations are performed after doing a Fourier transform with respect to the variable $x$: this facilitates the analytical calculations since it lowers the order of the differential equations in the ‘inner region’, as it has been first shown by Pegoraro & Schep (1986) and Pegoraro et al. (1989).

The boundary layer approach consists in solving the linear problem in distinct overlapping regions. The number of regions depends on the non-ideal parameters that exist in the problem. Moreover, for the boundary layer approach to be applied, it is necessary that the normalised non-ideal parameters (e.g. $\epsilon = S^{-1}$ or $\epsilon = (d_e/a)^2$) be much smaller than unity, i.e. $\epsilon \ll 1$. This makes it possible to perform an asymptotic analysis by expanding quantities in powers of $\epsilon$, and also grants the scale separation between the boundary layers, since the characteristic width of the innermost layers, where ideal-MHD breaks down, is a *posteriori* found to scale with a multiplication of positive powers of the parameters ‘$\epsilon$’ at play. Quantitatively speaking, comparison between theoretical estimates and numerical integration of the boundary layer problem indicates that the condition $\epsilon \ll 1$ typically means that the microscopic scales (cf. end of § 2.2) related to the non-ideal terms capable of breaking the ideal MHD conservation of the magnetic topology be not larger than a fraction $\sim 10^{-1}$ of the equilibrium shear length $a$.

For example, for purely resistive or purely inertia-driven tearing modes in RMHD, the linear problem is solved in *two* distinct regions and the solutions are asymptotically matched in an intermediate layer, as shown in figure 4(a). In this case, the first region is the ‘outer’ region in which the plasma is assumed to be at ideal MHD force-balance equilibrium at spatial scales of the order of $L_0 = a$. The second region is the ‘inner’ one, in which a solution is sought for the differential equations while retaining the non-ideal parameters in the limit $x \ll 1$ (i.e. $x/L_0 \ll 1$ in dimensional coordinates). The characteristic size of the inner region is identified to correspond to the *reconnection layer width* (cf. § 8), whose characteristic thickness we hereafter call $\delta$, with $\delta \to 0$ as the
<figure><img src="image_6.png" /><figcaption>FIGURE 4. (a) Sketch of the boundary layer decomposition for the calculation of the tearing eigenmodes with asymptotic matching techniques, when a single boundary layer is sufficient (e.g. case of the resistive tearing first considered by Furth *et al.* 1963). The width $\delta$ of the inner boundary layer is here the natural candidate for the width of the ‘reconnecting layer’, although the latter is not formally defined, i.e. we can assume $\delta = \delta_1$ (an assumption we will support in § 8 with further arguments). (b) Sketch of the boundary layer decomposition when two boundary layers are required (e.g. case of the warm-collisionless RMHD). Note that the reconnecting layer width, which we call $\delta$, is, again, a priori not defined but, as we will see (§ 8), further arguments can be found to identify it with the innermost layer width, i.e. $\delta = \delta_1$.</figcaption></figure>

non-ideal parameters tend to zero. The inner region is therefore defined by the inequality $x/\delta \lesssim 1$.

Seen from the inner region, the convergence to some point in the outer domain is expressed with the limit $x/\delta \to \infty$ with respect to the ‘stretched variable’ $x/\delta$. In the matching layer, this limit should match the limit $x/L_0 \to 0$ (or, simply, $x \to 0$ in normalised units) for quantities evaluated in the outer region.

When more than one non-ideal parameter is involved in the problem, applying the boundary layer approach becomes more difficult. This is due to the presence of at least two boundary layers and therefore of two matching regions (cf. figure 4b), whose characteristic widths, i.e. the normalisation scales for the asymptotic matching, say $\delta_2$ and $\delta_1$, depend on the spatial scales of the problem in a non-trivial way. Moreover, the notion of ‘reconnection layer width’, which we keep on naming $\delta$, in this case is not a priori related to the width of one of the boundary layers, although, as we will see (§ 8), we can make it correspond to the extension of the innermost layer, $\delta_1$.

Operationally speaking, the integration of tearing mode equations by means of the boundary layer approach, which we will detail in the next sections, can be split in the following sequence of steps.

(i) *Solving the equations in the ideal MHD limit* provides the solution in the outer layer (the ‘outer solution’) and allows one to evaluate the instability parameter (the so-called $\Delta'$ parameter, see § 3.1), which identifies the range of wavenumbers that are tearing-unstable and allows one to distinguish among different wavelength regimes of reconnection – the constant-$\psi$ or tearing regime, the large-$\Delta'$ or internal-kink regime, and the fastest growing mode in a continuum wavelength spectrum of slab tearing modes.

(ii) *Considering the equations in the non-ideal region* by assuming that the inverse spatial gradients are small enough to count the microscopic scales associated with the non-ideal parameters. This also allows some simplifications to the equations,
related to the fact that the nth derivative with respect to x dominates over the nth power of k.

(iii) Establishing, on heuristic basis, an ordering among terms of the non-ideal equations depending on the value of the non-ideal microscopic parameters. This allows one to distinguish different boundary layers (cf. figure 4b) in which the non-ideal equations are differently approximated. This task becomes more complex as more non-ideal parameters are present, and is subordered to the consistency of the solutions a *posteriori* found after integration in each sub-region (hence, the 'heuristic' nature of the ordering of the terms in the equations).

(iv) *Integrating the non-ideal equations* in each sub-region, by taking the solution in the next-most, outer domain as a boundary condition. For example, with reference to figure 4(b), the outer solution in the ideal region provides the boundary condition in the matching layer II for the solution in the sub-domain $|x| \lesssim \delta_2/2$, which in turn provides the boundary condition in the matching layer I for the solution in the region $|x| \lesssim \delta_1/2$.

(v) Finally note that, in any case, constructing the global solution in the whole domain in an explicit form is a not evident task and often it is not possible, since the integral solutions in each non-ideal region may be not obtained in closed form. Nevertheless, it is possible to obtain a quantitative estimate of the eigenvalues of the linear problem (i.e. the growth rate of the reconnecting mode) from the conditions on the solutions in the non-ideal region. Similarly, other spatial scales of interest (which, by further arguments, can be interpreted as corresponding, e.g. to the reconnecting layer width, $\delta$ - cf. § 8) can be quantitatively evaluated.

### 3.1. Orderings and instability parameter $\Delta'$

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

$$ \gamma \sim 0, \quad \mathcal{L} \sim 1, \quad \mathcal{F} \simeq 1, \quad \mathcal{A} \sim 1, \quad \mathcal{B} = \mathcal{C} = ik\psi_0' \sim 1. \qquad (3.1) $$

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,

$$ \psi_0(x)|_{x \to 0} \simeq C_0 + C_2x^2 + C_4x^4 + \dots \qquad (3.2) $$

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

$$ \psi'_{0}|_{x \to 0} \sim \psi'''_{0}|_{x \to 0} \sim x|_{\text{inner}} \simeq \delta. \qquad (3.3) $$

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 1982; Militello et al. 2004, 2011). 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:

$$ \mathcal{L} \approx \partial_x^2, \quad \mathcal{F} = 1 - d_e^2 \mathcal{L}, \quad A = -k\delta(1 - \mathcal{L}) \sim -k\delta \mathcal{L}, \quad B = ik\delta[(1 - d_e^2) - \rho_s^2 \mathcal{L}]. \quad (3.4) $$

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 a_n^{(\text{out})}(x/\delta)^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 (1978, § 9). This translates in the condition

$$ \lim_{x/\delta \to \pm\infty} a_m^{(\text{in})} (x/\delta)^m \approx \lim_{x \to \pm\infty} a_m^{(\text{out})} x^m \quad (\text{for } m \text{ that corresponds to the leading term of the series}). \quad (3.5) $$

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

$$ \lim_{x/\delta \to \pm\infty} \psi_{\text{in}}(x/\delta) = \lim_{x \to \pm\infty} \psi_{\text{out}}(x). \quad (3.6) $$

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

$$ \int_{-\infty}^{+\infty} \psi_{\text{in}}'' d(x/\delta) = \psi_{\text{out}}(0) \Delta', \quad (3.7) $$

which introduces the instability parameter (Furth et al. 1963) related to the discontinuity which is met in the first derivative of $\psi_{\text{out}}$ as $x \to \pm 0$,

$$ \Delta'(k; \psi_0) \equiv \lim_{\epsilon \to 0} \frac{\psi'_{\text{out}}(\epsilon) - \psi'_{\text{out}}(-\epsilon)}{\psi_{\text{out}}(0)}. \quad (3.8) $$

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

$$ \int_{-\infty}^{+\infty} \psi_{\text{in}}'' d(x/\delta) = \psi_{\text{in}}'(+\infty) - \psi_{\text{in}}'(-\infty) = \psi_{\text{out}}'(+0) - \psi_{\text{out}}'(-0) = \psi_{\text{out}}(0) \Delta', \quad (3.9) $$

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_{in}$ and $\varphi_{in}$. According to the inner region ordering, (2.2) becomes $\gamma\varphi_{in}'' = ikx\psi_{in}''$. Therefore,

$$ \int_{-\infty}^{+\infty} \psi_{in}'' d(x/\delta) = -(i\gamma/k) \int_{-\infty}^{+\infty} \frac{\varphi_{in}''}{x/\delta} d(x/\delta). \quad (3.10) $$

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. 1992), 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 (1986), Porcelli (1991) and Bulanov et al. (1992). 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_{norm}$, with respect to which to consider the limits $x/l_{norm} \ll 1$ and $x/l_{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_{norm}$ for $l_{norm} = L_0 = a$. At the same time, also the relations between $\Delta'$ and the spatial gradients of $\varphi_{in}$ expressed by combinations of (3.7) and (3.10) become non-trivial, as we will discuss in § 10.

## 3.2. Small- and large-$\Delta'$ limits in slab geometry and in tokamaks

For each eigenvalue problem in slab geometry, the numerical value of $\Delta'$ 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 $(\Delta')^{-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-$\Delta'$ limit for $\Delta'\delta \ll 1$ and of a large-$\Delta'$ limit for $\Delta'\delta \gg 1$.

Once the scaling of $\delta$ on the non-ideal parameters involved is known, the large- and small-$\Delta'$ 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 (2017) for the purely resistive case).

In particular, since $\Delta'$ 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

$$ \lim_{ka \to 0} \Delta'(ka) = (ka)^{-p}, \quad p > 0, \quad (3.11) $$

where $p$ depends on the initial equilibrium profile (Del Sarto et al. 2016; Pucci et al. 2018; Betar et al. 2020). While the limit above corresponds to $\Delta' \to \infty$, the marginal stability condition $a\Delta'(ka) \to 1$ is approached as $ka \to 1$.

A fundamental difference between large- and small-$\Delta'$ 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-$\Delta'$ 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-$\Delta'$ 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).

We remark indeed that the distinction between small- and large-$\Delta'$ limits has been historically introduced to characterise, in terms of the instability parameter $\Delta'$ defined by Furth *et al.* (1963), two different types of unstable modes observed in tokamaks. In tokamak physics, for positive values of $\Delta'$, a distinction is made between the *tearing mode* or *constant-$\psi$ mode*, first studied in the cylindrical geometry approximation by Furth *et al.* (1963) (and which formally corresponds to $\Delta'\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.* (1976) (and which formally has $\Delta' = \infty$) (we recall indeed that, even if $\Delta' < 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 (1990), § 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 \text{const.}$ of the plasma inside the inner region; more specifically, for the internal kink mode, $\xi \sim \text{const.}$ in the whole subdomain $r \le r_s$, where $r_s$ is the radius of the resonant surface, whereas $\xi = 0$ for $r > r_s$ (see, e.g. Porcelli 1987; Del Sarto & Ottaviani 2017 for a more detailed discussion). It has been however shown by Ara *et al.* (1978) 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-$\Delta'$ 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-$\Delta'$ 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 $\Delta'\delta \sim 1$ (see Loureiro, Schekochihin & Cowley 2007; Bhattacharjee *et al.* 2009; Del Sarto *et al.* 2016; Betar *et al.* 2020) 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 $\Delta'\delta \ll 1$ to $\Delta'\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.* 1963, appendix D, and also Biskamp 1982). 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 $\Delta'(k)$ formula has been discussed by Ottaviani & Porcelli (1993, 1995) and for which the scaling $\gamma_M \sim S^{-1/2}$ can be deduced (dashed line in the figure).
<figure><img src="image_7.png" /></figure>

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

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

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

$$q(r) = \left| \frac{B_{\varphi} r}{B_{\theta} R} \right| = - \frac{m}{n}$$

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

$m = 1, n = 1$  
$\Delta' = \infty$  
internal kink mode

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

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

large wavelength limit
$\Delta'\delta \gg 1$

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

$m = 2, n = 1$  
tearing mode  
$\Delta'\delta < 1$

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

small wavelength limit
$\Delta'\delta \ll 1$
single mode: $k_y = \frac{4\pi}{L_y}$

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

double (symmetric) mode
$k_y = \frac{2\pi}{L_y}$

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 1976). 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. 1973), 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 \epsilon_B$, in terms of the ratio $\epsilon_B$ between the in-plane and the guide magnetic field components at equilibrium (cf. (A11) and discussion in the last paragraph of Appendix A).
<figure><img src="image_17.png" /></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. (2016) the value $p=2$, which can be deduced from the $ka \ll 1$ limit of the corresponding $\Delta'(ka)$ formula (see Ottaviani & Porcelli 1993 or (8) of Ottaviani & Porcelli 1995). An 8th-order accurate scheme has been used here for integration on a uniform grid by using the solver discussed by Betar et al. (2020).

When tearing modes are destabilised in a large aspect ratio, *periodic* current sheet, so that also the $\Delta'\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 (1989) – see also figure 1 of Betar *et al.* (2020). In this scenario, the ‘soliton-like’ behaviour of magnetic islands, noted already by Biskamp (1982) in his early simulations, has earned them the name of ‘plasmoids’ (see, e.g. Biskamp & Welter 1989; Biskamp 1996; Shibata & Tanuma 2001; Loureiro *et al.* 2007; Bhattacharjee *et al.* 2009), 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 1967). 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 (1979) and Birn & Hones (1981) or later works such as Hautz & Scholer (1987), Scholer (1987) and also Lin, Cranmer & Farrugia (2008), 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 1983; Gosling & McComas 1987). In this sense, in the recent literature, the term
‘plasmoid regime’ of reconnection (Uzdensky, Loureiro & Schekochihin 2010; Huang & Bhattacharjee 2010; Comisso et al. 2018; Zhou et al. 2021) is typically used to identify both tearing reconnection in the $\Delta'$δ ~ 1 wavelength limit (cf. the early numerical study by Biskamp (1986) 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 1989; Biskamp & Bremer 1994; Wei et al. 2000; Servidio et al. 2009, 2011; Wan et al. 2013; Franci et al. 2017 and several more recent works) or by experimental measurements (see, e.g. Nishizuka et al. 2015; Kumar et al. 2019; Yan et al. 2022). 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 1997; Loureiro et al. 2007; Bhattacharjee et al. 2009). 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 (2014).

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 1971; Van Hoven & Cross 1971) and solar physics (see, e.g. Priest 1976; Velli & Hood 1989), and has been more recently been considered for applications to secondary instabilities to primary reconnection processes (see, e.g. Tajima & Shibata 1997; Shibata & Tanuma 2001; Tanuma et al. 2001; Loureiro et al. 2005; Daughton & Scudder 2006; Drake et al. 2006; Loureiro et al. 2007; Bhattacharjee et al. 2009; Landi et al. 2015; Tenerani et al. 2015; Del Sarto et al. 2016; Del Sarto & Ottaviani 2017; Papini, Landi & Del Zanna 2019b; Singh et al. 2019) and to reconnection in turbulence (see, e.g. Pucci & Velli 2014; Loureiro & Uzdensky 2016; Tenerani et al. 2016; Comisso et al. 2018; Papini et al. 2019a; Betar et al. 2020; Kowal et al. 2020; Pucci et al. 2020; Schekochihin 2020; Tenerani & Velli 2020; Franci et al. 2022) in the context of the so-called ‘turbulent (or turbulence-mediated or turbulence-driven) reconnection’ scenario (Matthaeus & Lamkin 1986; Strauss 1988; Loureiro et al. 2009; Matthaeus & Velli 2011; Schekochihin 2020) – 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. 2020).

In the remainder of the article, we just focus on the formal solution of the eigenvalue problems in a slab current sheet in the $\Delta'$δ ≪ 1 and $\Delta'$δ ≫ 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 (1986), Pegoraro et al. (1989), Schep et al. (1994), Basu & Coppi (1981) and Porcelli (1991). These calculations provide the available analytical formulae for the growth rates in the ‘internal-kink’ and ‘constant-ψ’ 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-$\Delta'$ limits can be obtained, has been written as in Del Sarto & Ottaviani (2017, appendix C) in terms of a characteristic scale length $\delta_L$ that coincides with the relevant layer width $\delta$ in the large- and small-$\Delta'$ regimes. In §§ 5–7, we will detail the analytical calculations that allow the relevant limits of these dispersion relations to be recovered.

Finally, we draw attention to the fact that the definition (3.8) in a slab geometry changes meaning as soon as the condition $\Delta'\delta \gtrsim 1$ is achieved, since application of (3.11) for $ka \ll 1$ and $\Delta'\delta \sim 1$ formally implies a dependence of $k$ on $\delta$. In general, indeed, $(\Delta')^{-1}$ becomes a microscopic scale as soon as $\Delta'L_0 \gg 1$, and for $\Delta'\delta \gtrsim 1$, the condition $(\Delta'(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 (1986) and Porcelli (1991), 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 (2011) 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 $E \times 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 1981); 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
<table><caption>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 (1991) and Ottaviani & Porcelli (1995) for $d_e^2 \gg \rho_s^2$ and that of Porcelli (1991) and Comisso et al. (2013) for $\rho_s > d_e \neq 0$ are reported. In both cases, we have adopted the notation used in Del Sarto & Ottaviani (2017, 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. 1978) and of the warm-electron resistive regimes (Pegoraro & Schep 1986) 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.</caption><thead><tr><th rowspan="2">RMHD<br>orderings</th><th colspan="4">'General' dispersion relation</th></tr><tr><th>(for Q = γ/k<sub>e</sub>, δ<sub>L</sub> = γ/k<sub>e</sub>Q<sup>1/2</sup>d<sub>e</sub>)</th><th>δ'<sub>LD</sub></th><th>γ<sub>LD</sub></th><th>δ<sub>SD</sub></th><th>γ<sub>SD</sub></th></tr></thead><tbody><tr><td>d<sub>e</sub> &gg; ρ<sub>s</sub></td><td>Δ'<sub>δ<sub>L</sub></sub> = -π/Q[Γ[(Q-1)/4]/Γ[(Q+5)/4]</td><td>d<sub>e</sub></td><td>kd<sub>e</sub></td><td>Δ'<sub>d<sub>e</sub></sub><sup>2</sup></td><td>(Δ')<sup>2</sup>kd<sub>e</sub><sup>3</sup></td></tr><tr><td>S<sup>-1</sup> > ρ<sub>s</sub>γ<sup>1/2</sup>γ<sup>1/4</sup></td><td>Same as above for<br>d<sub>e</sub> → (Sγ)<sup>-1/2</sup></td><td>k<sup>-4/7</sup> ×<br>ρ<sub>s</sub><sup>-5/7</sup>S<sup>-4/7</sup></td><td>k<sup>6/7</sup> ×<br>ρ<sub>s</sub><sup>4/7</sup>S<sup>-1/7</sup></td><td>k<sup>-2/3</sup>Δ<sup>1/3</sup> × (ρ<sub>s</sub>S)<sup>-2/3</sup></td><td>k<sup>2/3</sup>Δ<sup>2/3</sup> × ρ<sub>s</sub><sup>2/3</sup>S<sup>-1/3</sup></td></tr><tr><td>ρ<sub>s</sub> > d<sub>e</sub> &ne; 0</td><td>Δ'<sub>δ<sub>L</sub></sub> = π/d<sub>e</sub>ρ<sub>s</sub>(1 - π/d<sub>e</sub><sup>2</sup>Q<sup>3</sup>)</td><td>ρ<sub>s</sub><sup>-1/3</sup>d<sub>e</sub><sup>4/3</sup></td><td>kd<sub>e</sub><sup>1/3</sup>ρ<sub>s</sub><sup>2/3</sup></td><td>Δ'<sub>d<sub>e</sub></sub><sup>2</sup></td><td>Δ'<sub>kd<sub>e</sub>ρ<sub>s</sub></sub></td></tr><tr><td>ρ<sub>s</sub><sup>1/2</sup>γ<sup>1/4</sup> &gg; S<sup>-1</sup></td><td>Same as above for<br>for d<sub>e</sub> → (Sγ)<sup>-1/2</sup></td><td>k<sup>-1/3</sup>S<sup>-1/3</sup></td><td>k<sup>2/3</sup>S<sup>-1/3</sup></td><td>Δ<sup>1/5</sup> ×<br>(kS<sup>-2/5</sup>)</td><td>Δ<sup>4/5</sup> ×<br>k<sup>2/5</sup>S<sup>-3/5</sup></td></tr></tbody></table>
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 1992). By matching the asymptotic solutions, they then recovered the scaling laws first evaluated in the Fourier space (Porcelli 1991).

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. (1963) in the cold-electron regime and by Pegoraro & Schep (1986) in the warm-electron regime, and the collisionless scaling laws first evaluated by Coppi (1964c,a) in the cold-electron limit and by Porcelli (1991) 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. 2012b), and we compare some key points of that procedure with the integration procedure that we follow in this work.

## 4. Boundary layer solution in the outer ideal region

We now turn our attention to solve the linearised equations in the ideal MHD region.

In the outer region and sufficiently far from the reconnection layer, a steady force balance condition in ideal MHD can be assumed. Therefore, the terms weighted by the non-ideal parameters can be neglected. In this case, using $\psi_1 = \psi_{\text{out}}$, $\varphi_1 = \varphi_{\text{out}}$, (2.1)–(2.2) become

$$ \varphi_{\text{out}} = -\frac{i\gamma}{k\psi_0'} \psi_{\text{out}}, \qquad (4.1) $$

$$ \psi_{\text{out}}'' = \left( k^2 + \frac{\psi_0'''}{\psi_0'} \right) \psi_{\text{out}}. \qquad (4.2) $$

Equation (4.1) says that $\psi_1$ and $\varphi_1$ have opposite spatial parity and are de-phased by $\pi/2$, i.e. $\psi_1 \sim \cos(ky)$ versus $\varphi_1 \sim \sin(ky)$ if $y = 0$ is assumed at an O-point.

It is also evident that solving (4.2) gives the outer solution ($\psi_1 = \psi_{\text{out}}$) which is unequivocally determined by the equilibrium profile $\psi_0$ and by $k$. Using (4.1), one can so obtain also the eigenmodes $\varphi_{\text{out}}$. The number of inner regions depends instead on the number of non-ideal parameters involved in the problem.

A class of magnetic equilibria of particular interest is provided by those fulfilling the condition $\lim_{|x| \to \infty} (\psi_0'' / \psi_0') = \text{constant}$, for which a quite general integration procedure can be devised, as we are going to discuss in the specific example below.

### 4.1. A specific example of evaluation of $\Delta'(\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

$$ \psi_0'(x) = \frac{\tanh(x)}{\cosh^2(x)}, \qquad (4.3) $$

which is the one used by Betar et al. (2020) to obtain the scalings to which we will make comparison. This profile was first proposed by Porcelli et al. (2002). Substituting
the previous relation into (4.2), one finds

$$ \psi_{\text{out}}'' - \left( \alpha^2 - \frac{12}{\cosh^2(x)} \right) \psi_{\text{out}} = 0, \quad (4.4) $$

where $\alpha^2 = k^2 + 4$ and $\psi_0''/\psi_0' = 4 - 12/\cosh^2(x)$. Grasso first evaluated the corresponding $\Delta'$, reported by Porcelli et al. (2002), 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

$$ \psi_f'' - \alpha^2 \psi_f = 0, \quad (4.5) $$

whose solution is

$$ \psi_f(x) = A e^{-\alpha x}, \quad (4.6) $$

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

$$ \psi_{\text{out}}(x) = \psi_f(x)f(x) = A e^{-\alpha x} f(x). \quad (4.7) $$

Substituting the previous equation in (4.4), one obtains

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

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

$$ (1 - z^2)f'' - 2(z + \alpha)f' + 12f = 0, \quad (4.9) $$

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 1964, § 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

$$ f = a + bz + cz^2 + dz^3, \quad (4.10) $$

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

$$ \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\}. \quad (4.11) $$

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),

$$ \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\}, \quad (4.12) $$

we obtain, using (3.8),

$$ \Delta' = 2 \frac{-k^4 + 2k^2 + 15}{k^2 \sqrt{k^2 + 4}}. \quad (4.13) $$

It follows that $\Delta' > 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 1962) can be obtained in the same manner, as already noted by Furth (1963) (§ 4, therein). In this case, the equivalent of (4.9) is a Legendre equation of degree $p = 2$. Therefore, $f(z) = a + bz$, 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 $\Delta' = 2(1/k - k)$, where the instability condition $\Delta' > 0$ is met when $k \in [0, 1]$ (Furth et al. 1963). The point we underline here is that the solution of the outer equation and the evaluation of $\Delta'$ 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 $\Delta' > 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. 1963). In general, exponentially growing solutions of the linear problem are obtained when $\Delta' > 0$ (cf. figure 7).

## 5. Equations in the non-ideal region

After having solved the linearised equations in the ideal MHD region, we are now interested in finding their solutions in the non-ideal region. We thus begin this section by introducing a ‘generalised resistivity’, which allows us to cast the linear equations in a unique form that can be applied to both the collisionless and resistive regimes. For analytical purposes that will become evident in the following, in this section, we also normalise the linear system with respect to an arbitrary characteristic length – say $\ell$ – which will be later specified in each subdomain and in each regime. Then, we will introduce an auxiliary function $\chi$ related to both scalar fields $\psi_1$ and $\varphi_1$, which allows us to combine the inner layer equations into a single equation. It is this equation that, in some reconnection regimes, will be later approximated in each subdomain of the boundary layer approach and which will be solved analytically according to the strategy detailed in § 3.

For the analytical integration of the linear problem, it is first useful to introduce a parameter that includes both the resistive and inertia-related non-ideal terms in Ohm's law by thus highlighting their ‘almost symmetric’ contributions in (2.1). This can be equivalently done in terms of the ‘generalised resistivity’, as first suggested by Furth (1963) (§ 5 therein), which, in the normalised units we use, reads

$$ \bar{S}^{-1} \equiv S^{-1} + \gamma d_e^2, \quad (5.1) $$

or in the form of a ‘generalised electron skin depth’, as done by Porcelli (1991) (note that Porcelli (1991) uses the symbol $\Delta$ in place of the $\bar{d}_e$ we adopt here), by expressing the
Lundquist number in terms of the normalised electron-ion collision rate, $v_{ei}$, as $S^{-1} = d_e^2 v_{ei}/2$, so as to write

$$ \bar{d}_e^2 = d_e^2 + \frac{S^{-1}}{\gamma} = d_e^2 \left( 1 + \frac{1}{2} \frac{v_{ei}}{\gamma} \right). \qquad (5.2) $$

In either case, the second left-hand side term and the last right-hand side term of (2.1) can be re-absorbed into a single contribution that formally accounts for both resistive and inertial effects, the latter being therefore interpretable, in the form $\gamma d_e^2$, as a ‘collisionless resistivity’. Using for example (5.2), we rewrite (2.1) as

$$ \gamma[\psi_1 - \bar{d}_e^2(\psi_1'' - k^2\psi_1)] - ik\varphi_1(\psi_0' - d_e^2\psi_0'') = -ik\rho_s^2(\varphi_1'' - k^2\varphi_1)\psi_0'. \qquad (5.3) $$

In the following, when we will analytically integrate the linear problem in the non-ideal region, it will be more convenient to start from (2.2) and (5.3), of which we will then take the appropriate collisional or collisionless limits. Specialising the parameters $\bar{S}^{-1}$ or $\bar{d}_e$ to the purely resistive or purely collisionless limits is meaningful for linear tearing modes, because of the extremely narrow region of the parameter space in which both $S^{-1}$ and $d_e$ appreciably contribute to the asymptotic scalings (Betar et al. 2020). A transition from resistive to inertia-dominated regimes can instead be relevant to the onset of secondary collisionless modes over primary resistive tearing modes (Del Sarto & Ottaviani 2017), provided the applicability of the WKB approximation we spoke of in § 3.2 for secondary tearing instabilities developing on finite length, large-aspect ratio current sheets generated by the primary reconnection event. In this case, the possibility of having a transition of regime depends on both the rescaling of the shear length and of the magnetic field amplitude of reference (see Del Sarto et al. 2018).

A further distinction will be made between the ‘cold’ regime $\bar{d}_e^2 \gg \rho_s^2$ (discussed in § 6) and the opposite ‘warm’ regime $\rho_s^2 \gg \bar{d}_e^2$ (discussed in § 7). The reason why we are going to treat the two cases separately is due to the fact that the normalisation scales that define the extent, in an asymptotic sense, of the innermost boundary layer domain (i.e. $\delta_1$ of figure 4b) turn out to be different, depending on whether $\rho_s$ is negligible or not with respect to $d_e$ (or to $(S^{-1}/\gamma)^{1/2}$). The asymptotic scaling of the layer width, $\delta_1$, also determines, via the matching conditions, the scalings of the reconnection rate $\gamma$. It follows that the ‘cold’ collisionless solution cannot be obtained as a trivial limit $\rho_s \to 0$ of the solution obtained for the ‘warm’ case.

Operationally speaking, we will solve the inner equations using the generalised electron skin depth of (5.3), so we will first formally recover the collisionless results of Coppi (1964c), Coppi (1964a) and Porcelli (1991). Then, we will comment about the correspondence of these results with those obtained by Furth et al. (1963), Coppi et al. (1976), Ara et al. (1978), Pegoraro & Schep (1986) and Pegoraro et al. (1989) when resistivity dominates over electron inertia.

### 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 2xC_2 = x\psi_0''|_{x=0} = xJ_0$, where we have named $J_0 = \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

$$ \gamma \psi_1 - \gamma \bar{d}_e^2 \psi_1'' = ikJ_0 x \varphi_1 - ikJ_0 \rho_s^2 x \varphi_1'', \qquad (5.4) $$

$$ \gamma \varphi_1'' = ikJ_0 x \psi_1''. \qquad (5.5) $$
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

$$ \gamma \psi_1 = ikJ_0x\varphi_1 + \gamma \left( \bar{d}_e^2 + \rho_s^2 \frac{k^2 J_0^2}{\gamma^2} x^2 \right) \psi_1''. \quad (5.6) $$

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. 1986; Zocco & Schekochihin 2011) 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. (2012b) 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

$$ \frac{\rho_s^2}{\bar{d}_e^2} \ll \left( \frac{\gamma}{kxJ_0} \right)^2. \quad (5.7) $$

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 $\Delta'$-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 $\Delta'$ 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_v$ – would have appeared. If the hypothesis $\lambda_v \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 $v$ alone dominates, one where both $v$ 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

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

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

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

$$ \tilde{\varphi}_1'' = \zeta \tilde{\psi}_1'', \qquad (5.10) $$

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\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. 2020). 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 'χ' and large- and small-Δ' 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-Δ' limit, consists in casting the equations for $\tilde{\psi}_1$ and $\tilde{\phi}_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. 1976; Ara et al. 1978)

$$ \chi(\zeta) = \zeta \tilde{\psi}'_1(\zeta) - \tilde{\psi}_1(\zeta). \qquad (5.11) $$

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

$$ \tilde{\phi}'_1(\zeta) = \chi(\zeta) - \lim_{\zeta \to \infty} \chi(\zeta) = \chi(\zeta) - \chi_{\infty}. \qquad (5.12) $$

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{\phi}'_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. 1978). 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-Δ' limit

In the case of the small-Δ' 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

$$ \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 \text{ ordering}). \qquad (5.13) $$

In this case, we have used the 'constant-ψ' approximation (Furth et al. 1963), i.e. the fact that in the small-Δ' 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-ψ ordering can be assumed and heuristically verified after integration of the eigenmode problem for the single boundary layer limit (cf. Furth et al. 1963; White 1983), and in the small-Δ' regime, it can be shown to be equivalent to the condition $\delta_1 \Delta' \ll 1$. In this case, we can also estimate $\lim_{x \to 0} x \psi'_1(x) \sim \lim_{x \to 0} x \Delta' \psi_{\text{out}}(x) \sim \delta_1 \Delta' 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 $\Delta' \delta_1 \ll 1$ condition, the last passage of (5.13) is thus justified, whence one deduces

$$ \chi_{\infty} = -c_0, \quad \text{for } \delta_1 \Delta' \ll 1. \qquad (5.14) $$

#### 5.3.2. Asymptotic estimate of $\chi_\infty$ in the large-Δ' limit

In the large-Δ' limit for $x \to 0$, instead, one finds that $\chi_\infty = 0$. This has been first discussed by Coppi et al. (1976) and Ara et al. (1978), 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{\phi}_1(\zeta) = -\gamma\zeta$ for $|\zeta| \le 1$ and $\tilde{\phi}_1(\zeta) = 0$ for $|\zeta| > 1$ (cf. also Porcelli 1987; Del Sarto & Ottaviani 2017), 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):

$$ \psi_{\text{out}}(x) \sim c_0 + c_1|x| \implies \psi_{\text{out}} \sim c_0 + c_1|x| \implies \Delta' \sim \frac{2c_1}{c_0}. \qquad (5.15) $$

This solution formally holds in the outermost matching layer, where the instability parameter $\Delta'$ 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-$\Delta'$ limit, in which it can be a posteriori shown, from heuristic arguments (see Ottaviani & Porcelli (1995) 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 $\Delta'(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

$$ \lim_{\zeta \to \infty} \chi(\zeta) \sim c_1|x| - c_0 - c_1|x| = -c_0 \sim \frac{c_1}{\Delta'}. \qquad (5.16) $$

The fact that $c_0 \to 0$ as $\Delta'$ 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 $\Delta'(k) = -19.19$ (green curve), $\Delta'(k) = 48.9$ (blue curve) and $\Delta'(k) = 3.2766$ (orange curve), evaluated using (4.13) for $k = 10, 0.1$ and $1.7$, respectively. Thus, in the large-$\Delta'$ limit where $\Delta' \to \infty$, (5.16) implies that

$$ \chi_{\infty} = 0 \quad \text{for } \Delta' \to \infty. \qquad (5.17) $$

The fact that condition $\Delta' \to \infty$ can be replaced by $\Delta'$ 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 $\Delta'\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

$$ (\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}{G^2} \frac{\rho_s^2}{\ell^2} \zeta^2 (\zeta \tilde{\psi}'''_1 + \tilde{\psi}''_1) = -\frac{1}{G^2} \zeta^2 \tilde{\phi}'_1, \qquad (5.18) $$

which can be rewritten as

$$ \left( \frac{\bar{d}_e^2}{\ell^2} + \frac{\rho_s^2 \zeta^2}{G^2} \right) \chi'' - 2 \frac{\bar{d}_e^2}{\ell^2} \frac{\chi'}{\zeta} - \left( 1 + \frac{\zeta^2}{G^2} \right) \chi = - \frac{\zeta^2}{G^2} \chi_\infty. \qquad (5.19) $$

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
<figure><img src="image_18.png" /><figcaption>FIGURE 7. Profile of the eigenfunction $\psi_1$ in the purely resistive regime, computed via the numerical solver of Betar et al. (2020) 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 ($\Delta' = -19.19 < 0$, green curve); of the small-$\Delta'$ wavelength limit ($\Delta' = 3.2766$, red curve); and of the large-$\Delta'$ wavelength limit ($\Delta' = 48.9$, blue curve). In the stable case ($\Delta' < 0$), the eigenfunction does not change concavity as $x$ approaches the neutral line, whereas for the unstable tearing-type solutions obtained for $\Delta' > 0$, the value $c_0 = \lim_{x\to 0} \psi_1(x)$ decreases approaching zero as the value of $\Delta'(k)$ increases approaching $+\infty$.</figcaption></figure>

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-$\Delta'$ 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.

(i) In the small-$\Delta'$ 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 $\Delta'$ 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| \le \delta_1$ and in the possible intermediate region $\delta_1 < |x| \le \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 region<fnref n="2" /> (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.

(ii) In the *large-$\Delta'$ 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. Solutions in the inner region: cold regime

Let us first discuss the reconnection regime with both cold electrons and ions. The linear problem has been first solved in the small-$\Delta'$ wavelength limit by Furth *et al.* (1963) in the resistive regime and by Coppi (1964c,a) in the inertial regime, using boundary layer calculations in the coordinate space. A solution of the linear problem in the large-$\Delta'$ wavelength limit has been first provided in the resistive regime by Coppi *et al.* (1976) and has been extended in the collisionless regime by Porcelli (1991). In this latter work, the calculations in the Fourier space, with which Pegoraro & Schep (1986) also recovered the resistive results of Furth *et al.* (1963) and Coppi *et al.* (1976), have been extended to the purely collisionless regime. We will come back to these more recent works relying on the integration in the Fourier space in § 7, where we will discuss the ‘warm’ reconnection regime including FLR effects related to parallel electron compressibility (i.e. a kind of diamagnetic effect).

The reason for which we start from the case of cold species is that, from the point of view of the boundary layer theory, the ‘cold’ regime in which we neglect $\rho_s$-related terms in (2.1)–(2.2) is both conceptually and analytically simpler than the ‘warm’ tearing regime (cf. figure 4). Since we consider the cold limit $\bar{d}_e^2 \gg \rho_s^2$ which encompasses the case $\rho_s = 0$, to ensure convergence of the solution with the $\rho_s = 0$ limit, we assume the validity of the condition (5.7) in the whole integration domain. We can verify a *posteriori* its validity, once the eigenvalue problem is solved.

We also mention here that the analytical results we are going to obtain below have been numerically verified by Betar *et al.* (2020) to be valid for $\bar{d}_e^2 > \rho_s^2/100$, i.e. for $\bar{d}_e \gtrsim \rho_s/10$.

<footnote n="2">²In (Granier *et al.* 2022) it has been pointed out that in the warm small-$\Delta'$ regime the equation for $\phi_1$ changes of form while moving from the layer of width $\delta_2$ to the innermost layer of width $\delta_1$: while this would formally require to carry out the full boundary layer calculation for its integration, doing so only yields negligible additive corrections of order $O(\delta_1/\rho_s) \ll 1$ to the growth rate scalings.</footnote>
6.1. Solution for $\bar{d}_e^2/\rho_s^2 \gg 1$: large-$\Delta'$ limit

Assuming heuristically the validity of (5.7), we neglect the second term in parentheses in (5.9). The restriction to the large $\Delta'$-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):

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

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

$$ \delta_1 = \bar{d}_e, \qquad (6.2) $$

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

$$ \chi_1 = e^{-\alpha \zeta^2}. \qquad (6.3) $$

Substituting it in (6.1) yields

$$ 2\alpha = 1, \quad 4\alpha^2 = \frac{1}{\mathcal{G}^2}. \qquad (6.4) $$

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

$$ \gamma = k \bar{d}_e J_0. \qquad (6.5) $$

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. 2020) 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 1991)

$$ \gamma = k d_e J_0 \quad \delta_1 \sim d_e. \qquad (6.6) $$

In the resistive case, (6.6) specialises to the result of (Coppi et al. 1976)

$$ \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}. \qquad (6.7) $$

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-$\Delta'$ 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,

$$u''\chi_1 + 2u'\left(\chi_1' - \frac{\chi_1}{\zeta}\right) + u\chi_1'' = 0. \qquad (6.8)$$

Further substitution of (6.3) into the equation above yields

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

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

$$\frac{dw}{w} = 4\alpha\zeta \, d\zeta + \frac{2 \, d\zeta}{\zeta} \implies w(\zeta) = w_0 \zeta^2 e^{2\alpha\zeta^2}, \qquad (6.10)$$

$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

$$\chi_2 = i\sqrt{\frac{\pi}{2}} \operatorname{erf}\left(\frac{i\zeta}{\sqrt{2}}\right) e^{-(\zeta^2/2)} + \zeta, \qquad (6.11)$$

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-$\Delta'$ limit

In the small-$\Delta'$ limit, the boundary layer analysis mirrors that of the paper by Furth et al. (1963) 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{\phi}_1(\zeta)$ that can be brought to an integrable form more easily than the equation for $\chi(\zeta)$:

$$\tilde{\phi}_1'' - \frac{\zeta^2 \ell^2}{G^2 \bar{d}_e^2} \tilde{\phi}_1 = c_0 \frac{\ell^2}{\bar{d}_e^2} \zeta. \qquad (6.12)$$

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

$$\Phi'' - z^2 \Phi = z \qquad (6.13)$$

after the change of variables and coordinates

$$z = (G\bar{d}_e)^{-1/2} \zeta, \quad \Phi = \left(\frac{\bar{d}_e}{G^3}\right)^{1/2} \left(\frac{\tilde{\phi}_1}{c_0}\right). \qquad (6.14)$$

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. (1963), whose solution has been provided by Basu & Coppi (1976), Coppi et al. (1976) and Ara et al. (1978) (cf. (III.26) and (III.27) and appendix A in the latter) in the closed integral form:

$$ \Phi = -\frac{1}{2}z \int_{0}^{1} (1-t^2)^{-1/4} e^{-(1/2)t z^2} dt. \qquad (6.15) $$

The scaling of the eigenvalue $\gamma$ and its dependence on the amplitude of the instability parameter $\Delta'$ can be made explicit by comparing the eigenfunction $\Phi(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

$$ \psi_{\text{out}}(0) \Delta' = -i \frac{\gamma}{kJ_0} \int_{-\bar{X}'}^{+\bar{X}'} \frac{\varphi''(x')}{x'} dx', \qquad (6.16) $$

where $\Delta'$ 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:

$$ \left. \frac{d^2 \tilde{\varphi}_1}{d\zeta^2} \right|_{\ell=1} = c_0 \left( \frac{1}{\bar{d}_e \mathcal{G}^3} \right) \left. \left( \frac{dz}{d\zeta} \right)^2 \right|_{\ell=1} \frac{d^2 \Phi}{dz^2} = c_0 \left( \frac{\mathcal{G}}{\bar{d}_e^3} \right) \left. \left( \frac{d^2 \Phi}{dz^2} \right) \right|_{\ell=1} \frac{d^2 \Phi}{dz^2}. \qquad (6.17) $$

This authorises us to identify

$$ x' = \frac{x}{\delta_1} \to z = \frac{x}{(\mathcal{G}\bar{d}_e)|_{\ell=1}^{1/2}}, \quad -i \frac{\gamma}{kJ_0} \frac{d^2\varphi(x')_1}{dx'^2} \to c_0 \left(\frac{\mathcal{G}}{\bar{d}_e^3}\right) \left.\left(\frac{d^2\Phi}{dz^2}\right)\right|_{\ell=1}^{1/2} \frac{d^2\Phi}{dz^2}, \qquad (6.18) $$

which now fixes the boundary layer width to be

$$ \delta_1 = (\mathcal{G}\bar{d}_e)|_{\ell=1}^{1/2} = \sqrt{\frac{\gamma \bar{d}_e}{k}}, \qquad (6.19) $$

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

$$ \Delta' = \left( \frac{\mathcal{G}}{\bar{d}_e^3} \right) \left. \int_{-\bar{Z}}^{+\bar{Z}} \frac{\Phi''}{z} dz, \qquad (6.20) $$

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_{-\infty}^{+\infty} (\Phi''/z) dz \simeq \int_{-\bar{Z}}^{\bar{Z}} (\Phi''/z) dz$. Then, using the definition in (5.8) for $\ell = 1$,
one obtains

$$ \gamma = k \bar{d}_e^3 (\Delta')^2 \left( \frac{J_0}{I^2} \right). \qquad (6.21) $$

In the inertia-dominated regime, $\bar{d}_e = d_e$, one recovers the well-known result by Coppi (1964c,a), later re-obtained by Porcelli (1991):

$$ \gamma = k (\Delta')^2 d_e^3 \left( \frac{J_0}{I^2} \right), \quad \delta_1 = \frac{\Delta' d_e^2}{I}. \qquad (6.22) $$

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

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

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. Solutions in the inner region: warm-electron regime

Let us now consider the boundary layer problem in the ‘warm’ reconnecting regime by limiting our attention to the case of cold ions, consistent with (2.1) and (2.2).

Extending the analysis to include ion temperature effects related to FLR corrections in the fluid description of these low frequency modes is possible, but generally it requires relying on reduced models (i.e. ‘gyrofluid’ models), in which the ion response related to the separation between the motion of particles and of gyrocentres is obtained by starting from the gyrokinetic Vlasov equation. Fluid models of reconnection, in which this has been done, typically differ because of the way the gyrokinetic operator is approximated and because of the number of scalar fields that are retained in the description (see, e.g. Pegoraro & Schep 1981, 1986; Pegoraro et al. 1989; Porcelli 1991; Schep et al. 1994; Loureiro & Hammett 2008 for essentially two-field models, and, e.g. Aydemir (1992), Zakharov & Rogers (1992), Dorland & Hammett (1993), Smolyakov, Pogutse & Hirose (1995), Pogutse, Smolyakov & Hirose (1998), Waelbroeck, Hazeltine & Morrison (2009), Zocco & Schekochihin (2011), Waelbroeck & Tassi (2012), Connor et al. (2012b) and Tassi, Sulem & Passot (2016) for other gyrofluid models retaining more than two scalar fields). A common feature of linear tearing analysis build on these gyrofluid models, as well as on other descriptions retaining kinetic features at a more fundamental level (see Drake & Lee 1977; Cowley et al. 1986; Rogers et al. 2007), is the ‘symmetric’ way by which ion and electron temperature effects enter in the dispersion relation in regimes where an estimate of the growth rate of tearing-type modes is analytically available. This symmetry is expressible via the formal substitution $\rho_s^2 \to \rho_s^2 + \rho_i^2$ in the dispersion relations which we are going to discuss in the remainder of this work. In regimes of the parameter space that go beyond the applicability of the analytical predictions by Pegoraro & Schep (1986), Pegoraro et al. (1989) and Porcelli (1991), a departure from the symmetry between $\rho_s$ and ion-Larmor-radius effects associated with the parameter $\rho_i = \rho_s \sqrt{T_i/T_e}$ has been numerically observed in the growth rate of the internal-kink mode (Del Sarto et al. 2011). Dealing with these subjects would include some complicacy about the build of the linear model, which are unnecessary with respect to the purpose of the present work. This is why we neglect them here.

We therefore start from (5.4) and (5.5). Since the problem has now two characteristic spatial scales, namely $\bar{d}_e$ and $\rho_s$, the identification of the two boundary layers will require
us to perform two subsequent normalisations and approximations of the terms of the aforementioned equations. Since we are now restricting to the parameter range where $\rho_s^2$ is not negligible with respect to $\bar{d}_e^2$, we can assume $\rho_s^2 \gtrsim \bar{d}_e^2$. A natural choice for choosing a first normalisation scale defining the ‘largest’ non-ideal region of interest is therefore to pose $\ell = \rho_s$. We note that identifying this as a boundary layer width or not is a further logical step, whose appropriateness depends on the kind of approximations of the equations and integration strategies that can be possibly performed with respect to the stretched variable $x/\rho_s$. This is what we are going to discuss below.

Also note that for economy of symbols – let us say – in this section, we will use some symbols used in the changes of variables of § 6 with a different definition that will be given in each specific case (this will concern, in particular, the definition of $\zeta$, $\mathcal{G}$ and the use of the ‘$\cdot\cdot\cdot$’ to label some normalised quantities).

### 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

$$ \tilde{\psi}_1 - \left( \frac{\bar{d}_e^2}{\rho_s^2} + \frac{\tilde{\zeta}^2}{\tilde{G}^2} \right) \tilde{\psi}''_1 = - \frac{\tilde{\zeta}}{\tilde{G}^2} \tilde{\phi}_1, \quad \tilde{\phi}''_1 = \tilde{\zeta} \tilde{\psi}''_1, \qquad (7.1) $$

with

$$ \tilde{\zeta} = \frac{x}{\rho_s}, \quad \tilde{G} = \frac{\gamma}{k\rho_s J_0}, \quad \tilde{\phi}_1 = -i\tilde{G}\phi_1(\tilde{\zeta}), \quad \tilde{\psi}_1 = \psi_1(\tilde{\zeta}). \qquad (7.2) $$

The auxiliary equation (5.19) specialises to

$$ \left( \frac{\bar{d}_e^2}{\rho_s^2} + \frac{\tilde{\zeta}^2}{\tilde{G}^2} \right) \tilde{\chi}'' - 2 \frac{\bar{d}_e^2 \chi'}{\rho_s^2 \tilde{\zeta}} - \left( 1 + \frac{\tilde{\zeta}^2}{\tilde{G}^2} \right) \tilde{\chi} = - \frac{\tilde{\zeta}^2}{\tilde{G}^2} \tilde{\chi}_{\infty}. \qquad (7.3) $$

Equation (7.3) can be shown to be equivalent to (B42) of Zocco & Schekochihin (2011), 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_{\parallel,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 (2011) 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_{\parallel,e}$ and $A_{\parallel,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 $\Delta'$ and what they named the 'reconnection region' ($\delta_{in}$, in their notation). As we will later discuss (§ 8), the characteristic width of their 'reconnection region', $\delta_{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 $(\Delta'_v)^{-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 2011) used a perturbative method to obtain a closed form solution for $A_{\parallel,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 \ge \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. 2020) 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

$$ \tilde{\psi}_1 - \frac{\tilde{\zeta}^2}{\tilde{G}^2} \tilde{\psi}_1'' = - \frac{\tilde{\zeta}}{\tilde{G}^2} \tilde{\phi}_1, \quad \tilde{\phi}_1'' = \tilde{\zeta} \tilde{\psi}_1'' \quad (7.4) $$

and (7.3) reads

$$ \tilde{\chi}'' - \left( 1 + \frac{\tilde{G}^2}{\tilde{\zeta}^2} \right) \tilde{\chi} = - \tilde{\chi}_{\infty}. \quad (7.5) $$

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_{out}$ of the ideal region via (3.5) and (3.6) formally identifies the outermost boundary layer of width

$$ \delta_2 = \rho_s. \quad (7.6) $$

### 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 $G$) if $(\zeta\rho_s)^2/(G\bar{d}_e)^2 \sim 1$. In dimensional
units, this means $xk\rho_sJ_0/(\gamma\bar{d}_e) \sim 1$, which suggests that we introduce the normalisation scale $\ell = \gamma\bar{d}_e/(k\rho_sJ_0)$. Using then in (5.9) and (5.10) the change of coordinates and variables expressed by

$$ \zeta = \frac{x}{\ell} = \frac{xk\rho_sJ_0}{\gamma\bar{d}_e}, \quad \mathcal{G} = \frac{\rho_s}{\bar{d}_e}, \quad \varphi_1 = -i\mathcal{G}\varphi_1(\zeta), \quad \psi_1 = \psi_1(\zeta), \qquad (7.7) $$

we write

$$ \psi_1 - \frac{\bar{d}_e^2}{\ell^2} (1 + \zeta^2) \psi_1'' = -\frac{\zeta}{\mathcal{G}^2} \varphi_1, \quad \varphi_1'' = \zeta \psi_1'', \qquad (7.8) $$

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

$$ \psi_1 - \frac{\bar{d}_e^2}{\ell^2} (1 + \zeta^2) \psi_1'' = 0, \quad \varphi_1'' = \zeta \psi_1'', \qquad (7.9) $$

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

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

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

$$ \delta_1 = \frac{\gamma \bar{d}_e}{k \rho_s J_0}, \qquad (7.11) $$

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 $\tilde{\zeta}$ introduced in (7.1) and (7.2).

(i) The innermost region $x \le \delta_1$, where $\delta_1 = \gamma\bar{d}_e/(k\rho_sJ_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$.

(ii) The intermediate region $\delta_1 \ll x \le \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$.
(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{\xi} \gg 1$, i.e. for $x \sim a$.

Two matching regions can be therefore identified.

(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_{\tilde{\zeta} \to \infty} \chi(\tilde{\zeta}) = \lim_{\tilde{\zeta} \to 0} \tilde{\chi}(\tilde{\zeta})$.

(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{\phi}_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{\phi}_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 $\Delta'$, 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

$$ \psi_1 = \frac{1}{2\pi i} \oint_C F(s)(1 + \zeta^2)^s \, ds, \qquad (7.12) $$

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

$$ \psi_1'' = \oint_C 2s(2s - 1)F(s)(1 + \zeta^2)^{s-1} \, ds - \oint_C 4s(s - 1)F(s)(1 + \zeta^2)^{s-2} \, ds. \qquad (7.13) $$

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

$$ \psi_1'' = \oint_C 2s(2s - 1)F(s)(1 + \zeta^2)^{s-1} \, ds - \oint_C 4s(s + 1)F(s + 1)(1 + \zeta^2)^{s-1} \, ds. \qquad (7.14) $$

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

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

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

$$ F(s + 1) = \frac{\left(s - \frac{\alpha}{2}\right) \left(s - \frac{\beta}{2}\right)}{s(s + 1)} F(s), \quad (7.16) $$

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

$$ F(s + j) = \frac{\left(s - \frac{\alpha}{2}\right)_j \left(\frac{s - \beta}{2}\right)_j}{(s)_j (s + 1)_j} F(s), \quad (7.17) $$

where

$$ (q)_j = q(q+1)\cdots(q+j-1) = \frac{\Gamma(q+j)}{\Gamma(q)}. \quad (7.18) $$

is the so-called 'Pochhammer symbol' defined in terms of the $\Gamma$-function (see Abramowitz & Stegun 1964, § 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,

$$ F(s) = \frac{\Gamma(s+a)\Gamma(s+b)}{\Gamma(s+c)\Gamma(s+d)}. \quad (7.19) $$

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

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

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

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

which has the following solutions:

$$ 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}. \quad (7.22) $$

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

$$ F(s) = \frac{\Gamma\left(s - \frac{\alpha}{2}\right) \Gamma\left(s - \frac{\beta}{2}\right)}{\Gamma(s) \Gamma(s + 1)}. \quad (7.23) $$

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

$$ s = \frac{\alpha}{2} - m, \quad s = \frac{\beta}{2} - m, \quad m = 0, 1, 2, 3, \dots \quad (7.24) $$

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,$$

$$ \begin{aligned} \psi_{\alpha} &= \mathrm{Res}_{\alpha} \equiv \sum_{m=0}^{\infty} \lim_{s \to \alpha/2-m} \left[ (s+m) \Gamma\left(s - \frac{\alpha}{2}\right) \right] \frac{\Gamma(-\nu - m)(1 + \zeta^2)^{\alpha/2-m}}{\Gamma\left(\frac{5}{4} - \frac{\nu}{2} - m\right) \Gamma\left(\frac{1}{4} - \frac{\nu}{2} - m\right)} \\ &= \sum_{m=0}^{\infty} \frac{(-1)^m}{m!} \frac{\Gamma(-\nu - m)}{\Gamma\left(\frac{5}{4} - \frac{\nu}{2} - m\right) \Gamma\left(\frac{1}{4} - \frac{\nu}{2} - m\right)} (1 + \zeta^2)^{\alpha/2-m}, \end{aligned} \quad (7.25) $$

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

$$ \begin{aligned} \lim_{s \to -m} (s+m) \Gamma(s) &= \lim_{s \to -m} (s+m) \frac{\Gamma(s+1)}{s} \\ &= \lim_{s \to -m} (s+m) \frac{\Gamma(s+m+1)}{s(s+1) \cdots (s+m-1)(s+m)} = \frac{(-1)^m}{m!}. \end{aligned} \quad (7.26) $$

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

$$ \psi_{\beta} = \mathrm{Res}_{\beta} = \sum_{m=0}^{\infty} \frac{(-1)^m}{m!} \frac{\Gamma(\nu - m)}{\Gamma\left(\frac{5}{4} + \frac{\nu}{2} - m\right) \Gamma\left(\frac{1}{4} + \frac{\nu}{2} - m\right)} (1 + \zeta^2)^{\beta/2-m}. \quad (7.27) $$

Therefore, using the residue theorem $(\oint_C f(z) dz = 2\pi i \sum_{z \in C} \mathrm{Res}[f(z)])$, and substituting in (7.12), $\psi_1$ becomes

$$ \psi_1 = \psi_{\alpha} + \psi_{\beta}. \quad (7.28) $$

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 1978, § 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-$\Delta'$ limit

We now solve the auxiliary equation (7.10) in the large-$\Delta'$ 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 \le \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

$$ \chi(\zeta) = \frac{1}{2\pi i} \oint_C G(s)(1 + \zeta^2)^s \, ds. \qquad (7.29) $$

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

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

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

$$ G(s) = \frac{\Gamma(s+a)\Gamma(s+b)}{\Gamma^2(s+1)} = \frac{F(s)}{s}. \qquad (7.31) $$

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 $\Gamma^2(s+1) = s\Gamma(s)\Gamma(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

$$ \chi_{\text{in}}(\zeta) = \chi_{\alpha}(\zeta) + \chi_{\beta}(\zeta), \qquad (7.32) $$

where

$$ \chi_{\alpha}(\zeta) = \sum_{m=0}^{\infty} \frac{(-1)^m}{m!} \frac{\Gamma(-\nu - m)}{\Gamma^2\left(\frac{5}{4} - \frac{\nu}{2} - m\right)} (1 + \zeta^2)^{\alpha/2-m}, \qquad (7.33) $$

$$ \chi_{\beta}(\zeta) = \sum_{m=0}^{\infty} \frac{(-1)^m}{m!} \frac{\Gamma(\nu - m)}{\Gamma^2\left(\frac{5}{4} + \frac{\nu}{2} - m\right)} (1 + \zeta^2)^{\beta/2-m}. \qquad (7.34) $$

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

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

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

$$ \chi_{\text{in}}(\zeta) \approx \chi_{\text{in}}^{\text{appr}}(\zeta) = A\zeta^{\alpha} + B\zeta^{\beta} + C. \qquad (7.36) $$
### 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$:

$$ \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), \quad (7.37) $$

$$ \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), \quad (7.38) $$

$$ \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). \quad (7.39) $$

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

### 7.4.3. Solution in the intermediate region, $\delta_1 \ll x \le \delta_2$

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

$$ \tilde{\chi}'' - \left(1 + \frac{\tilde{G}^2}{\tilde{\zeta}^2}\right) \tilde{\chi} = 0, \quad (7.40) $$

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 1964, § 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

$$ \tilde{\zeta}^2 u'' + \tilde{\zeta} u' - (\tilde{\zeta}^2 + \nu^2) u = 0, \quad (7.41) $$

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}) = K_\nu(\tilde{\zeta})$. Therefore, (7.40) has the solution

$$ \tilde{\chi}_{\text{mid}}(\tilde{\zeta}) = \tilde{\zeta}^{1/2} K_{\nu}(\tilde{\zeta}), \quad (7.42) $$

the label ‘mid’ standing for the intermediate region domain $\delta_1 \ll x \le \delta_2$ where it is defined. Its asymptotic behaviour as $\tilde{\zeta} \to 0$ is given by

$$ \begin{aligned} \tilde{\chi}_{\text{mid}}^{\text{aprr}} &\equiv \lim_{\tilde{\zeta} \ll 1} \tilde{\chi}_{\text{out}}(\tilde{\zeta}) = \lim_{\tilde{\zeta} \ll 1} \tilde{\zeta}^{1/2} K_{\nu}(\tilde{\zeta}) = \lim_{\tilde{\zeta} \ll 1} \tilde{\zeta}^{1/2} \left(\frac{\pi}{2}\right) \frac{\mathcal{I}_{-\nu}(\tilde{\zeta}) - \mathcal{I}_{\nu}(\tilde{\zeta})}{\sin(\pi\nu)} \\ &\sim \frac{2^{\nu}}{\Gamma(1-\nu)} \tilde{\zeta}^{\alpha} - \frac{2^{-\nu}}{\Gamma(1+\nu)} \tilde{\zeta}^{\beta}, \end{aligned} \quad (7.43) $$

where we used the relation in § 9.6.2 of Abramowitz & Stegun (1964). 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,

$$ \mathcal{I}_\nu(\tilde{\zeta}) \equiv \left(\frac{\tilde{\zeta}}{2}\right)^\nu \sum_{m=0}^\infty \frac{\left(\frac{\tilde{\zeta}^2}{4}\right)^m}{m!\Gamma(\nu + m + 1)}. \quad (7.44) $$
### 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 0} \chi(\zeta) = \lim_{\zeta \to 0} \tilde{\chi}(\tilde{\zeta})$, that is, $\lim_{\zeta \to \infty} \chi_{\text{in}}^{\text{appr}}(\zeta) = \lim_{\zeta \to 0} \tilde{\chi}_{\text{mid}}^{\text{appr}}(\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

$$A = \frac{2^{\nu}}{\Gamma(1-\nu)} \left(\frac{\delta_1}{\rho_s}\right)^{\alpha}, \quad B = -\frac{2^{-\nu}}{\Gamma(1+\nu)} \left(\frac{\delta_1}{\rho_s}\right)^{\beta}. \qquad (7.45)$$

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

$$\chi_{\text{in}}(0) = A + B + C, \qquad (7.46)$$

$$\chi''_{\text{in}}(0) = A\alpha + B\beta, \qquad (7.47)$$

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

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

$$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). \qquad (7.49)$$

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

$$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. \qquad (7.50)$$

Combining the previous equation with (7.45), one obtains

$$-\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{\Gamma(1+\nu)}{\Gamma(1-\nu)}, \qquad (7.51)$$

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

$$\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{\Gamma\left(\frac{3}{2}\right)}{\Gamma\left(\frac{1}{2}\right)} \left(\frac{\rho_s}{\delta_1}\right). \qquad (7.52)$$

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

$$\delta_1 \sim \rho_s^{-1/3} \bar{d}_e^{4/3} \qquad (7.53)$$
whence, using the definition in (7.11), one obtains

$$ \gamma \sim k \rho_s^{2/3} \bar{d}_e^{4/3}. \qquad (7.54) $$

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 (1986) for the growth rate (cf. (33) therein and also (20) in Pegoraro et al. 1989) and of Zocco & Schekochihin (2011) for $\delta_1$ (cf. $\delta_\eta$ in (B100), therein),

$$ \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}. \qquad (7.55) $$

Specialising the same equations to the purely collisionless regime, $\bar{d}_e = d_e$, one recovers instead the result of Porcelli (1991) 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 2005),

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

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, G$ and $\delta_1$ in (7.35), one obtains for the inner layer

$$ \lim_{x \ll \delta_1} \psi_1 \approx \frac{2^{1/2}}{\Gamma(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}}{\Gamma(3/2)} \left(\frac{\delta_1}{\rho_s}\right) \left(1 + \frac{x^2}{\delta_1^2}\right)^{1/2} \\ + \left(\frac{\delta_1}{\rho_s}\right) \left(1 + \frac{\bar{d}_s^2}{\delta_1^2}\right). \qquad (7.57) $$

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

$$ \lim_{x \ll \delta_1} \varphi_1 \approx \sum_{m=0}^{\infty} \left[ \frac{2^{1/2}}{\Gamma(1/2)} \left(\frac{\delta_1}{\rho_s}\right)^{-\delta_1^2/\bar{d}_s^2} \left(-\frac{\delta_1^2/\bar{d}_s^2}{m}\right)_b - \frac{2^{-1/2}}{\Gamma(3/2)} \left(\frac{\delta_1}{\rho_s}\right) \left(-\frac{1}{m}\right)_b \right] \left(\frac{x}{\delta_1}\right)^{2m+1} \\ + \left(\frac{\delta_1}{\rho_s}\right) \left(1 + \frac{\bar{d}_s^2}{\delta_1^2}\right) \left(\frac{x}{\delta_1}\right), \qquad (7.58) $$

where

$$ \binom{r}{m}_b = \frac{r(r-1)\cdots(r-m+1)}{m!} = \frac{(r)_m}{m!} \qquad (7.59) $$

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

Notice that for $x \le \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. (2020) and using the definitions in (7.36). We can this way numerically verify that $\chi_{\text{in}}^{\text{appr}} \approx \psi_1$ within the innermost layer with $|x| \ll \delta_1$, in agreement with the previous discussion.

The eigenfunctions (7.57) and (7.58) agree with the behaviour of the perturbed current density profile provided by Pegoraro & Schep (1986) 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. (2020).
<figure><img src="image_19.png" /><figcaption>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. (2020) in the large-$\Delta'$ limit are shown in the figure.</figcaption></figure>

### 7.5. Solution of the linear equations for $\rho_s^2/d_e^2 \gg 1$: small-$\Delta'$ limit

In the small-$\Delta'$ limit, we use again the constant-$\psi$ approximation, as we already did in the integration of the cold-electron limit in § 6.2. Figure 9 shows the appropriateness of this assumption also in the warm-electron regime: it displays the profile of $\psi_1$ numerically computed with the eigenvalue solver of Betar et al. (2020) in the collisionless regime for $d_e = 2 \times 10^{-2}$, $\rho_s = 5 \times 10^{-3}$ and for a value of $k = 1.7$, which corresponds to $\Delta'(k) = 3.3$ and $\delta_1 \simeq 2.4 \times 10^{-6}$ (evaluated according to the definition in (8.3) – see later).

Therefore, we write $\psi_1 = c_0$ in the whole non-ideal region and, as discussed in § 5.5, this allows us to more conveniently perform the integration by looking at the equations for $\psi_1$ and $\varphi_1$ rather than the auxiliary equation. In this way, (7.1) can be combined in the form

$$ \left( 1 + \frac{\rho_s^2 \tilde{\zeta}^2}{\tilde{d}_e^2 \tilde{G}^2} \right) \tilde{\varphi}_1'' = \frac{\rho_s^2 \tilde{\zeta}^2}{\tilde{d}_e^2 \tilde{G}^2} \tilde{\varphi}_1 + c_0 \frac{\rho_s^2}{\tilde{d}_e^2} \tilde{\zeta}. \qquad (7.60) $$

The results of § 8.3 in the small-$\Delta'$ cold-electron limit are recovered by neglecting the $\rho_s^2 \tilde{\zeta}^2 \tilde{\varphi}_1''/(\tilde{d}_e^2 \tilde{G}^2)$ contribution in the left-hand side term of (7.60). Its presence here makes an intermediate matching region also appear in this wavelength limit.

This can be seen in figure 10, where the results of a numerical integration performed in the small-$\Delta'$ limit of the warm collisionless regime ($\tilde{d}_e = d_e$) with the solver of Betar et al. (2020) are shown. The blue curve corresponds to the $\rho_s^2 \tilde{\zeta}^2 \tilde{\varphi}_1''/(\tilde{d}_e^2 \tilde{G}^2)$ term in (7.60). Its strong variation close to the neutral line is concentrated in a narrow region much
<figure><img src="image_20.png" /><figcaption>FIGURE 9. Spatial profile of $\psi_1$, numerically computed with the eigenvalue solver of Betar et al. (2020) in the collisionless regime for a choice of the parameters $d_e = 2 \times 10^{-2}$, $\rho_s = 5 \times 10^{-3}$, and for a value of $k = 1.7$. These values correspond to $\Delta'(k) = 3.3$ and to $\delta = \delta_1 \simeq 2.4 \times 10^{-6}$, evaluated according to the definition in (8.3) introduced next. The smaller box inserted at the top right-hand corner of the figure shows a zoom of the solution inside the innermost layer of width $\sim \delta = \delta_1$, in which the appropriateness of the constant-$\psi_1$ condition is evidenced.</figcaption></figure>

smaller than $\delta_2 = \rho_s$, which is of the order of $\delta_1$ (evaluated from the numerically computed eigenfunction, as it will be discussed in § 8). The green curve corresponds instead to the $\tilde{\phi}_1''$ term in (7.60), i.e. to the first left-hand side contribution between parentheses, whereas the orange curve corresponds to the first term on the right-hand side of (7.60), $\rho_s^2 \tilde{\zeta}^2 \tilde{\phi}_1 / (\bar{d}_e^2 \tilde{G}^2)$. This term can be neglected with respect to the one corresponding to the blue curve. For completeness, the contribution of the last right-hand side term of (7.60), $c_0 \rho_s^2 \tilde{\zeta}^2 / \bar{d}_e^2$, is shown as the black, dashed line. From the comparison of these contributions, one verifies the appropriateness of the hypotheses justifying (7.4).

Thus, an intermediate region, corresponding to the interval $\delta_1 < |x| \lesssim \delta_2$ of § 7.2.1, can be recognised: here, (7.4) hold, which we can combine in the form

$$ \tilde{\zeta} \tilde{\phi}_1'' = \tilde{\zeta} \tilde{\phi}_1 + c_0 \tilde{G}^2. \qquad (7.61) $$

However, as discussed in § 5.5, it is not necessary here to first perform the integration of (7.61), since the constant-$\psi$ condition holds in this case also in the innermost region.

### 7.5.1. Matching and scalings in the non-ideal region $x \ll 1$

Combining (7.8) is in practice equivalent – except for the chosen normalisation – to neglecting the $\sim \tilde{\zeta}^2 \tilde{\phi}_1$ term in (7.60). Using the latter, where we recall lengths are normalised to $\rho_s$, we write

$$ \tilde{\phi}_1'' = c_0 \frac{\rho_s^2 \tilde{\zeta}}{\bar{d}_e^2 (1 + \rho_s^2 \tilde{\zeta}^2 / (\bar{d}_e^2 \tilde{G}^2))}. \qquad (7.62) $$
<figure><img src="image_21.png" /><figcaption>FIGURE 10. Comparison of the spatial profiles of the terms $\rho_s^2/(\bar{d}_e^2\tilde{G}^2) \times \tilde{\xi}^2\tilde{\varphi}_1''$ (blue curve), of $\rho_s^2/(\bar{d}_e^2\tilde{G}^2) \times \tilde{\xi}^2\tilde{\varphi}_1$ (orange curve) and of $\tilde{\varphi}_1$ (green curve) of (7.60), in a region of the domain with width $|x| < \delta_2 = \rho_s$. The parameters of the numerical integration performed with the eigen-solver of Betar et al. (2020) in the small-$\Delta'$ limit are shown at the top of the figure. The value of $\delta = \delta_1$ in the bottom right corner has been evaluated according to the definition in (8.3).</figcaption></figure>

The equivalent of (3.10), obtained after integration of the second of (7.1) (cf. it also with (6.16)) reads

$$ \int_{-\bar{Z}}^{+\bar{Z}} \tilde{\psi}_1'' d\tilde{\xi} = \int_{-\bar{Z}}^{+\bar{Z}} \frac{\tilde{\varphi}_1''}{\tilde{\xi}} d\tilde{\xi}, \qquad (7.63) $$

where $\tilde{Z} = X/\rho_s \gg 1$ is a matching point in the interval $\delta_1 \ll X \ll \delta_2$ in units of $L_0 = a$. Use of (3.9) allows one to establish the matching conditions with the outer solution by making explicit the dependence on $\Delta'$: writing therefore $\int_{-\bar{Z}}^{+\bar{Z}} \tilde{\psi}_1'' = c_0\Delta'$ and substituting (7.62) in (7.63), one obtains

$$ c_0\Delta' = \int_{-\bar{Z}}^{+\bar{Z}} \frac{\tilde{\varphi}_1''}{\tilde{\xi}} d\tilde{\xi} = c_0 \frac{\rho_s^2}{\bar{d}_e^2} \int_{-\bar{Z}}^{+\bar{Z}} \frac{d\tilde{\xi}}{1 + \rho^2 \tilde{\xi}^2} = c_0 \frac{\rho_s^2}{\varrho \bar{d}_e^2} \arctan(\varrho \tilde{\xi}) \Big|_{-\bar{Z}}^{+\bar{Z}}, \qquad (7.64) $$

where $\varrho \equiv \rho_s/(\bar{d}_e\tilde{G})$. Taking the asymptotic limit $\tilde{Z} \to \infty$ and using $\arctan(\tilde{\rho}\tilde{\xi})|_{-\infty}^{+\infty} = \pi$, (6.16) reads

$$ c_0\Delta' = c_0\pi \frac{\rho_s^2 \tilde{G}}{\bar{d}_e} \implies \gamma = \sqrt{\pi J_0 k \rho_s \bar{d}_e \Delta'} \implies \delta_1 \sim \Delta' \bar{d}_e^2, \qquad (7.65) $$

where in the last passage, we have used the definition in (7.11) of $\delta_1 = \gamma \bar{d}_e/(k\rho_s J_0)$. We can thus explicate the scalings above in the purely collisionless and purely resistive limits.
In the former, taking $\bar{d}_e = d_e$, one obtains

$$ \gamma \sim J_0 k \Delta' \rho_s d_e, \quad \delta_1 \sim \Delta' d_e^2, \qquad (7.66) $$

that is, the result of Porcelli (1991) (cf. the discussion between (8) and (9), therein). Taking instead $\bar{d}_e = S^{-1/2}\gamma^{-1/2}$, one obtains

$$ \gamma \sim J_0^{2/3} (k\Delta')^{2/3} \rho_s^{2/3} S^{-1/3}, \quad \delta_1 \sim J_0^{-1/3} k^{-2/3} (\Delta')^{1/3} \rho_s^{-2/3} S^{-2/3}, \qquad (7.67) $$

that is, the result first obtained by Pegoraro & Schep (1986) for the growth rate (cf. (39), therein) and by Zocco & Schekochihin (2011) for $\delta_1$ (which maps into $\delta_\eta$ of (B.96), therein).

### 7.5.2. Approximated eigenmodes in the non-ideal region $x \ll 1$

Approximated expressions for the eigenmodes can be obtained by integrating (7.62) twice, and by calculating the second of these integrals by parts: one finds the following approximate formula for the eigenfunction $\tilde{\varphi}_1(x/\delta_2)$ in a spatial interval $x \lesssim \delta_2 \simeq \rho_s$ that covers the whole non-ideal region:

$$ \varphi_1(x) \approx \frac{c_0 \delta_1^3}{2 \bar{d}_e^2} \left\{ \left( \frac{x}{\delta_1} \right) \ln \left( 1 + \left( \frac{x}{\delta_1} \right)^2 \right) + 2 \arctan \left( \frac{x}{\delta_1} \right) - 2 \left( \frac{x}{\delta_1} \right) \right\}. \qquad (7.68) $$

Substituting (7.68) into (7.62), using the second part of (7.8) and integrating twice by parts, gives an approximated expression for $\psi_1$:

$$ \psi_1(x) \approx c_0 + \frac{c_0 \delta_1^2}{\bar{d}_e^2} \left\{ \left( \frac{x}{\delta_1} \right) \arctan \left( \frac{x}{\delta_1} \right) - \frac{1}{2} \ln \left( 1 + \left( \frac{x}{\delta_1} \right)^2 \right) \right\}. \qquad (7.69) $$

From the expressions above, it is straightforward to find closed forms for both the magnetic field components and the current density perturbations, using their definitions $B_x(x) \equiv k\psi_1$, $B_y(x) \equiv -\partial\psi_1/\partial x$, and $J_z(x) \equiv -\delta^2\psi_1/\partial x^2$ and (7.69). These lead to

$$ B_x(x) = -k \sin(ky) \psi_1(x), \quad B_y(x) = -\frac{k \delta_1 c_0}{\bar{d}_e^2} \arctan \left( \frac{x}{\delta_1} \right), \qquad (7.70) $$

$$ J_z(x) = -\frac{c_0}{\bar{d}_e^2} \frac{1}{1 + \left(\frac{x}{\delta_1}\right)^2}. \qquad (7.71) $$

In the appropriate limits, this equation corresponds to (74) of Pegoraro & Schep (1986) and to (31) of Betar et al. (2020).

## 8. Operational definition of the reconnecting layer width, $\delta$, and other microscopic scales related to the eigenmodes

We have seen (§§ 6 and 7) that the boundary layer calculations allow one to determine the asymptotic scaling of the width of the layers, where approximate solutions of the eigenfunctions can be analytically evaluated. These characteristic widths appear as normalisation scales $\delta_1$ and $\delta_2$, which, although sometimes can be seen as ‘natural’, being suggested by the comparison of some terms in the equations (as in the cases we previously considered), are a priori arbitrary and are essentially determined by the algebra.