2.1. Preliminary considerations

Let the unperturbed reconnecting magnetic field be directed along the $y$ axis, $\boldsymbol {B}=B_0f(x)\hat {\boldsymbol {y}}$, where $f(x)$ is a magnetic field profile that depends only on the $x$ coordinate (see figure 1a). In such a geometry, the magnetic field is represented by the vector potential $\boldsymbol {A}=(0,0,A_z(x))$. In equilibrium, there is no electric field, and electrons and positrons drift in the $z$ direction at a velocity ${\pm }U$ (‘$+$’ for positrons and ‘$-$’ for electrons, $U>0$). The pair plasma under consideration is ultrarelativistic; therefore, the approximate dispersion law for plasma particles $\mathcal {E}\approx pc$ is used.

Figure 1. Magnetic field configuration within the current sheet: (a) unperturbed; (b) with the tearing mode perturbation.

We assume that the current sheet thickness is much greater than the Larmor radius of the particles; therefore, with the exception of the neutral layer at small $x$, the magnetohydrodynamic (MHD) approximation is valid throughout the entire space, regardless of the form of the distribution function. Namely, the static equilibrium is determined by the balance of the magnetic pressure and the effective particle pressure found as the appropriate moment of the distribution function (Alfven & Falthammar Reference Alfven and Falthammar1963, pp. 217–218). One can use a prepared magnetic field profile, and it is assumed that the distribution function of the plasma is adjusted to this field. For simplicity, we consider the current sheet with magnetic field profile $f(x)=\tanh (x/L)$ (Harris Reference Harris1962). Near the neutral layer, the MHD approximation is violated, and kinetic effects are important. Thus, one can find solutions in the inner and outer regions and then match them together, which ultimately gives us an expression for the growth rate of the tearing instability.

The general picture of collisionless tearing instability has long been known (e.g. Galeev & Sudan Reference Galeev and Sudan1984, pp. 309–315; Kadomtsev Reference Kadomtsev1987). Let us consider the equilibrium of the current sheet as a whole. The plasma pressure near the neutral plane $x=0$ is balanced with the magnetic pressure at infinity, i.e. $P_\text {tot}(0)=B_0^2/8{\rm \pi}$. The plasma pressure can be written as $P_\text {tot}(0)=(2/3)n_0\langle \mathcal {E}\rangle$, where $\langle \mathcal {E}\rangle$ is the average energy of particles. We also can consider the same balance but in terms of forces, i.e. $\boldsymbol {\nabla } P_\text {tot}=(1/c)\boldsymbol {j}\times \boldsymbol {B}$. Roughly, one can estimate $\boldsymbol {\nabla } P_\text {tot}\sim P_\text {tot}(0)/L$, where $L$ is the characteristic thickness of the considered layer and the total current density can be estimated as $j_z\sim 2en_0 U$. As a result, we obtain the following two equations that describe the current sheet equilibrium:

(2.1a,b)\begin{equation} B_0^2=\frac{16{\rm \pi}}{3} n_0\langle\mathcal{E}\rangle,\quad \frac{U}{c}\sim \frac{\langle\mathcal{E}\rangle}{eB_0 L}. \end{equation}

The second equation is just a definition of the characteristic diamagnetic drift velocity of the plasma near the neutral plane (see Appendix A for additional details). It is important to note that high-energy particles play a significant role only at $\alpha <2$ since, in this case, the average energy is $\langle \mathcal {E}\rangle \sim \mathcal {E}_\text {max}$, and at $\alpha \geq 2$ we have $\langle \mathcal {E}\rangle \sim \mathcal {E}_\text {min}$ (i.e. the high-energy tail is suppressed).

The current sheet that is described above can be considered as a set of straight filaments with a current. Since the currents in filaments flow in the same direction, they are attracted to each other, so such an equilibrium system with a current sheet is unstable to pinching (Coppi et al. Reference Coppi, Laval and Pellat1966). In fact, this is possible only if the resistivity of the plasma is not zero. Even small non-zero resistivity allows the currents to move, and any small displacement of a filament gives rise to a force imbalance. As a result, the individual current filaments merge, leading to a reconnection of the magnetic field and the formation of magnetic islands (see figure 1b).

Effective resistivity arises in the collisionless plasma because the particles are not magnetised near the plane $x=0$, so they can gain energy in the electric field. The periodic electric field, $\delta \boldsymbol {E}$, is induced near the neutral plane $x=0$ due to the magnetic reconnection. Indeed, according to Lenz's law, the electric field is directed in the negative $z$ direction in the vicinity of X-points, where plasma flows towards the neutral plane, increasing the magnetic flux, and in the positive $z$ direction between the X-points, where plasma flows away from the neutral plane. In turn, the particles are not magnetised in the layer with characteristic thickness (Parker Reference Parker1957)

(2.2)\begin{equation} l\sim \sqrt{2r_{g0} L}, \end{equation}

where $r_{g0}=p_\perp c/eB_0$ is the particle Larmor radius in the asymptotic magnetic field. Inside this layer, particles perform a meandering motion along the current sheet and, therefore, they are accelerated by the induced electric field. The acceleration time is ${\sim }1/kc$, because at a larger time, when the particle moves along the $y$ axis, it ‘sees’ an oscillating field and does not gain energy on average.

Thus, one can estimate the effective conductivity within the inner region by using the Drude formula with the effective collision time $\tau _\text {coll}\sim 1/kc$ (Artsimovich & Sagdeev Reference Artsimovich and Sagdeev1979, p. 242), which gives

(2.3)\begin{equation} \sigma^\text{eff}\sim \frac{2n_0 e^2}{m}\frac{1}{k c} \left\langle\frac{mc^2}{\mathcal{E}}\right\rangle. \end{equation}

Here, angular brackets denote averaging over particle energies. The factor $\langle mc^2/\mathcal {E}\rangle$ is just the average inverse Lorentz factor of particles; it reflects the fact that conductivity in a relativistic plasma is inversely proportional to the average relativistic mass $m_{\rm eff}=\mathcal {E}/c^2$. Since the averaged energy of the particles is in the denominator, the non-magnetised particles with the lowest energy make the largest contribution to the conductivity. According to this, we can assume that $\sigma ^\text {eff}$ operates only within the ‘inner’ region of size

(2.4)\begin{equation} l_\text{min}\sim \sqrt{\frac{2L}{eB_0\langle\mathcal{E}^{{-}1}\rangle}}. \end{equation}

Let us emphasise that in our case, the flux-freezing condition is violated because there is no magnetic field near the neutral layer, and the particles there are not magnetised such that they decouple from the fluid motion. This is what determines the thickness $l$ of the inner layer where the particles can be accelerated by the electric field. Here, the effective resistivity appears simply because the particles have inertia and ‘lazily’ respond to the electric field when moving within the inner region during the ‘collisional’ time $1/kc$.

2.2. The outer region

Before estimating the growth rate of tearing instability, let us consider the outer region $|x|>l_\text {min}$, where the plasma can be considered as an ideal magnetised fluid with infinite conductivity $\sigma \rightarrow \infty$. Indeed, since the Larmor radius of the particles here is much smaller than the characteristic thickness of the current sheet $L$, one can use the ideal MHD approximation. Even though local thermodynamic equilibrium is not valid in the considered problem, and strictly speaking we cannot use MHD equations, tearing instability develops slowly, so at each moment of time, we can consider the problem as quasi-static and consider the plasma velocity as a small perturbation. The validity of this statement is also confirmed by the fact that the kinetic theory gives the same equations (see Appendix A).

Slow plasma motions create the current density perturbation, which, in turn, creates a disturbance of the magnetic field, such that the total magnetic field can be represented as $\boldsymbol {B}(x,y)=B_0f(x)\hat {\boldsymbol {y}}+\delta \boldsymbol {B}(x,y)$. The magnetic field perturbation $\delta \boldsymbol {B}$ can be expressed via the vector potential perturbation $\delta \boldsymbol {A}=\delta A_z\hat {\boldsymbol {z}}$ according to $\delta \boldsymbol {B}=\text {rot}\,\delta \boldsymbol {A}$. Considering the magnetic field profile $f(x)=\tanh (x/L)$, Ampère's equation is reduced to (e.g. Sturrock Reference Sturrock1994; Boldyrev & Loureiro Reference Boldyrev and Loureiro2018)

(2.5)\begin{equation} \frac{\text{d}^2}{\text{d}\kern0.07em x^2}\delta A_z+\left[{-}k_y^2+\frac{2}{L^2\cosh^2(x/L)}\right]\delta A_z=0. \end{equation}

Therefore, the even-parity vector potential perturbation in the outer layer is (Sturrock Reference Sturrock1994; Boldyrev & Loureiro Reference Boldyrev and Loureiro2018)

(2.6)\begin{equation} \delta A_z(x)=\delta A_z(0)\left(1+\frac{1}{kL}\tanh\frac{|x|}{L}\right)\exp({-}k|x|), \end{equation}

where $k=|k_y|$. Odd-parity solutions $\delta A_z(-x)=-\delta A_z(x)$ are responsible for kink modes and are not of interest to us. It is clearly seen from the solution (2.6) that the vector potential perturbation is continuous at $x=0$ for all $kL$, but the derivative $\delta A'_z(x)$ is discontinuous in the general case. Let us introduce the notation

(2.7)\begin{equation} \varDelta'(0)=\frac{\delta A'_z(0+)-\delta A'_z(0-)}{\delta A_z(0)}\equiv \frac{\delta B_y(0+)-\delta B_y(0-)}{\text{i}(c/\omega)\delta E_z(0)}. \end{equation}

According to (2.6), we obtain

(2.8)\begin{equation} \varDelta'(0)=\frac{2(1-k^2L^2)}{kL^2}. \end{equation}

Therefore, $\varDelta '(0)=0$ and there is no discontinuity $[\![\delta B_y]\!]=0$ if and only if $kL= 1$ (see figure 2). At $kL\neq 1$ the magnetic field perturbation jump occurs. To remove the discontinuity of the magnetic field, we have to take into account the current of non-magnetised particles in the region $|x|< l_\text {min}$.

Figure 2. The normalised vector potential perturbation and the total magnetic field for the Harris current sheet at $kL=0.4$ and $kL=1$.

2.3. Tearing instability growth rate

Let us denote the current density of particles undergoing meandering motion near the neutral plane as $\delta j_z$. We assume for simplicity that $\delta j_z$ vanishes at $|x|>l_\text {min}$. According to Ampère's law, the magnetic field jump $[\![\delta B_y]\!]=\delta B_y(0+)-\delta B_y(0-)$ can be found as

(2.9)\begin{equation} [\![\delta B_y]\!]\approx \frac{4{\rm \pi}}{c} \int_{{-}l_\text{min}}^{l_\text{min}}\delta j_z\,\text{d}\kern0.07em x. \end{equation}

Therefore, $\delta j_z$ eliminates discontinuity and provides a smooth transition of the perturbed magnetic field across the plane $x=0$.

According to Ohm's law, $\delta j_z=\sigma ^\text {eff}\delta E_z$. Therefore, the current density of non-magnetised particles is directed along the electric field $\delta E_z$. On the other hand, the sign of the magnetic field jump is determined by the outer solution, and according to (2.7), one can write $[\![\delta B_y]\!]=(c/\gamma )\varDelta '(0)\delta E_z(0)$, where $\gamma =-\text {i}\omega$ is the growth rate of the instability. At $kL>1$, the direction of the current coincides with the sign of the magnetic field jump only if $\gamma <0$ and the instability is suppressed. Therefore, the tearing mode is arising unstable only for $kL<1$.

Further, it is convenient to express the electromagnetic fields via the vector potential perturbation $\delta A_z$:

(2.10a,b)\begin{equation} \delta E_z=\frac{\text{i}\omega}{c}\delta A_z\equiv{-}\frac{\gamma}{c}\delta A_z, \quad \delta B_y={-}\frac{\text{d}}{\text{d}\kern0.07em x}\delta A_z. \end{equation}

Assuming that the electric field is approximately constant within the inner region of size ${\sim }2l_\text {min}$, we obtain from (2.9) that

(2.11)\begin{equation} \gamma\approx \frac{c^2}{8{\rm \pi}}\frac{\varDelta'(0)}{\sigma^\text{eff} l_\text{min}}. \end{equation}

It should be noted that this result can be used for different magnetic field profiles. One has to calculate $\varDelta '(0)$ for the given profile and substitute it into the above equation.

Substituting the effective conductivity (2.3) and expression for $\varDelta '(0)$ into (2.11) and using relations for the current sheet equilibrium (2.1a,b) we obtain

(2.12)\begin{equation} \gamma(k)\sim\frac{c}{2L}(1-k^2L^2)\left(\frac{U}{c}\right)^{3/2} [\langle \mathcal{E}\rangle \langle \mathcal{E}^{{-}1}\rangle]^{{-}1/2}. \end{equation}

In the general case, characteristic energies $\langle \mathcal {E}\rangle$ and $1/\langle \mathcal {E}^{-1}\rangle$ may not coincide. At large $\alpha \geq 2$, we expect $\langle \mathcal {E}\rangle \sim \langle \mathcal {E}^{-1}\rangle ^{-1}\sim \mathcal {E}_\text {min}$ and the growth rates are determined by the same formula as for the ultrarelativistic Maxwellian plasma (Zelenyi & Krasnoselskikh Reference Zelenyi and Krasnoselskikh1979). On the other hand, in the case of the power-law distribution with $\alpha <2$, one can expect that $\langle \mathcal {E}\rangle \sim \mathcal {E}_\text {max}$ and $\langle \mathcal {E}^{-1}\rangle ^{-1}\sim \mathcal {E}_\text {min}$, i.e. there are two characteristic energy scales. Then, the growth rate contains a small factor $\epsilon ^{1/2}$, where

(2.13)\begin{equation} \epsilon=\frac{\mathcal{E}_\text{min}}{\mathcal{E}_\text{max}}. \end{equation}

According to (2.12), the growth rate reaches a non-zero constant value when $k$ goes to zero. However, since the conductivity $\sigma ^\text {eff}\sim 1/k$ decreases with $k$, the effective resistance disappears, and there should be $\gamma \rightarrow 0$ in the limit $k\rightarrow 0$. Therefore, (2.12) is invalid at small $k$. Indeed, in this case, we cannot assume that the electric field is constant within the region of size $\sim l_\text {min}$, such that $\delta E_z(l_\text {min})\approx \delta E_z(0)$.

To find the correct asymptotic, one has to solve the Maxwell equations within the inner region. By using zero divergence of the magnetic field and Faraday's law,

(2.14a,b)\begin{equation} \frac{\partial\delta B_x}{\partial x}+\text{i}k \delta B_y=0,\quad \delta B_x=\frac{ck}{\omega}\delta E_z, \end{equation}

one can rewrite Ampère's law in terms of the electric field perturbation only:

(2.15)\begin{equation} \frac{\text{d}^2}{\text{d}\kern0.07em x^2}\delta E_z-\frac{4{\rm \pi}\gamma}{c^2}\sigma^\text{eff}\delta E_z=0, \end{equation}

where the limit $k\rightarrow 0$ and the boundary condition $\text {d}(\delta E_z)/\text {d}x=0$ at $x=0$ are assumed. Thus, the electric field perturbation profile within the inner region can be written as

(2.16)\begin{equation} \delta E^{(\text{i})}_z(x)\approx \delta E_z^{(\text{e})}(l)\frac{\cosh \eta x}{\cosh \eta l}, \end{equation}

where $\eta =\sqrt {4{\rm \pi} \gamma \sigma ^\text {eff}/c^2}$ and $l=l_\text {min}$. Performing integration in (2.9), we obtain the dispersion equation at $k\rightarrow 0$:

(2.17)\begin{equation} \varDelta'(l)=2\eta\tanh\eta l, \end{equation}

where now the tearing parameter is

(2.18)\begin{equation} \varDelta'(l)=\frac{\delta B_y(l)-\delta B_y({-}l)}{\text{i}(c/\omega)\delta E_z(l)}. \end{equation}

Taking into account the outer solution (2.6), in the limit $k\rightarrow 0$ we obtain

(2.19a,b)\begin{equation} \delta A_z(x)\approx \frac{\delta A_z(0)}{kL}\tanh\frac{|x|}{L},\quad \varDelta'(l)\approx \frac{2}{L}\frac{\text{sech}^2(l/L)}{\tanh(l/L)}\approx\frac{2}{l}. \end{equation}

This leads to the following dispersion equation:

(2.20)\begin{equation} \eta l\tanh \eta l\approx 1. \end{equation}

In order of magnitude, this equality is satisfied when the tangent argument is of the order of unity; therefore

(2.21)\begin{equation} \gamma(k)\sim\frac{c^2}{4{\rm \pi}}\frac{1}{\sigma^\text{eff}(l_\text{min})^2}\sim \frac{1}{3}ck \left(\frac{U}{c}\right). \end{equation}

This expression gives the correct asymptotic behaviour for the growth rate of instability at $k\rightarrow 0$.

Obviously, different magnetic profiles lead to different scalings of the corresponding tearing instability growth rates. In the context of MHD turbulence, a periodic magnetic field $B(x)\sim B_0\sin (x/L)$ is more appropriate (Boldyrev & Loureiro Reference Boldyrev and Loureiro2018). However, dependence on $\alpha$ and $\mathcal {E}_\text {min}/\mathcal {E}_\text {max}$ is unchanged for any magnetic field profile with the same asymptotic behaviour at $x\rightarrow 0$. Therefore, in this work, we focus only on the usual Harris layer.

3.1. Numerical procedure

Let us expand the potential $\delta A_{1z}$ in the system of basis functions

(3.7)\begin{equation} \delta A_{1z}(x)=\sum_n \alpha_n\phi_n(x). \end{equation}

The basis functions can be chosen as pyramid functions (Burkhart & Chen Reference Burkhart and Chen1989; Daughton Reference Daughton2003):

(3.8)\begin{equation} \phi_n(x)= \begin{cases} (x-x_{n-1})/(x_n-x_{n-1}), & x_{n-1}\leq x\leq x_n, \\ (x_{n+1}-x)/(x_{n+1}-x_{n}), & x_{n}\leq x\leq x_{n+1}, \\ 0, & \text{otherwise}. \end{cases} \end{equation}

The basis function at the neutral plane is

(3.9)\begin{equation} \phi_1(x)= \begin{cases} (x-x_1)/x_1, & 0\leq x\leq x_1, \\ 0, & \text{otherwise}, \end{cases} \end{equation}

and the last basis function at large $x_N\gg L$ is

(3.10)\begin{equation} \phi_N(x)= \begin{cases} (x-x_{N-1})/(x_N-x_{N-1}), & x_{N-1}\leq x\leq x_N, \\ 0, & \text{otherwise}. \end{cases} \end{equation}

Another popular choice for basis functions is the Hermite polynomials, $H_n(\xi x)$, with the exponential weight function and $\xi =1/L$ or $\xi =k$ (Daughton Reference Daughton1999; Pétri & Kirk Reference Pétri and Kirk2007). However, these functions are less appropriate for small $k$. In this case, the vector potential is practically constant throughout the entire space and varies sharply in a small region near the neutral plane. Such behaviour could hardly be fitted with the Hermite polynomials.

One can substitute the expansion of the vector potential perturbation in (3.4) and multiply on $\phi _m(x)$. Integrating over $x$, we obtain the matrix equation

(3.11)\begin{equation} \sum_n \alpha_n G_{nm}=0, \end{equation}

where the matrix elements are

(3.12)\begin{align} G_{nm}(\omega,k_y) & = k_y^2\int\phi_{n}(x)\phi_m(x)\,\text{d}\kern0.07em x-\int \frac{\text{d}^2\phi_n}{\text{d}\kern0.07em x^2}\phi_m\,\text{d}\kern0.07em x \nonumber\\ & \quad +\int V_\text{ad}(x)\phi_n(x)\phi_m(x)\,\text{d}\kern0.07em x+ \frac{\text{i}\omega}{|k_y|c}\sum_s K_{nm} \end{align}

and

(3.13)\begin{equation} K_{nm}=\frac{8{\rm \pi}^2 q_s^2}{c^3}\int\frac{\text{d}\mathcal{E}\text{d}P_z}{T(\mathcal{E},P_z)} \mathcal{E}\frac{\partial f_{0s}}{\partial \mathcal{E}}\oint\text{d}\kern0.07em x \frac{p_z(x)}{|p_x(x)|} \phi_m(x)\oint\text{d}\kern0.07em x'\frac{p_z(x')}{|p_x(x')|}\phi_n(x'). \end{equation}

We changed the order of integration in $K_{mn}$; therefore, integrating $p_z(x)/|p_x(x)|$ at fixed $\mathcal {E}$ and $P_z$ should be performed only along a particle trajectory, where $|p_x|\geq 0$.

Non-trivial solutions, $\alpha _n$, exist if and only if $\det \boldsymbol {G}(\omega,k_y)=0$. Following Burkhart & Chen (Reference Burkhart and Chen1989), we consider only symmetric solutions, i.e. we consider only symmetric trajectories and perform integration over $x$ and $x'$ in the first quadrant. This means that if we consider some trajectory at $x>0$, there is the mirrored one, which gives the same contribution to $K_{nm}$. However, one should be careful with ‘crossing’ orbits: in such an approach, we can miss ‘half’ the trajectory (which belongs to $x<0$). Therefore, when integrating ‘crossing’ orbits in the domain $x>0$ only, it is necessary to multiply the result by 2.

One can classify all orbits as ‘crossing’ and ‘non-crossing’ by using integrals of motion only. In our problem, we have three integrals of motion: $y$ component of the mechanical momentum $p_y$ (which is approximately zero for resonance particles), $z$ component of the canonical momentum $P_z$ and the energy $\mathcal {E}$. If $P_z<0$ and $\mathcal {E}^2< c^2P_z^2$, a particle has a ‘non-crossing’ orbit; at other values of $P_z$ and $\mathcal {E}$ a particle has ‘crossing’ orbits. Obviously, the contribution from electrons and positrons is the same; therefore, one can calculate the kernel only for positrons and multiply the result by 2.