The linear dispersion relation for collisionless kinetic tearing instabilities is calculated for the Harris equilibrium. In contrast to the conventional 2D geometry, which considers only modes at the center of the current sheet, modes can span the current sheet in 3D. Modes at each resonant surface have a unique angle with respect to the guide field direction. Both kinetic simulations and numerical eigenmode solutions of the linearized Vlasov-Maxwell equations have recently revealed that standard analytic theories vastly overestimate the growth rate of oblique modes. We find that this stabilization is associated with the density-gradient-driven diamagnetic drift. The analytic theories miss this drift stabilization because the inner tearing layer broadens at oblique angles sufficiently far that the assumption of scale separation between the inner and outer regions of boundary-layer theory breaks down. The dispersion relation obtained by numerically solving a single second order differential equation is found to approximately capture the drift stabilization predicted by solutions of the full integro-differential eigenvalue problem. A simple analytic estimate for the stability criterion is provided.

## I. INTRODUCTION

Tearing instabilities are magnetic reconnection events that occur in transition layers where the frozen flux condition of ideal magnetohydrodynamics (MHD) is violated.^{1} In tokamaks, they contribute to the basic magnetic field topology,^{2} transport,^{3–6} rotation,^{7} and may trigger disruptive events.^{8} In space, they are responsible for initiating large-scale magnetic reconnection in Earth's magnetosphere,^{9} as well as flares in the solar corona.^{10} If the plasma is sufficiently dense or cool, tearing instabilities are well described by resistive,^{1,13–17} or visco-resistive^{18} MHD. At lower densities or higher temperatures, two-fluid or Hall MHD effects modify the resistive reconnection rate.^{19–23} At even lower densities or higher temperatures, resistivity plays no role and reconnection is determined entirely by collisionless two-fluid or kinetic effects.^{24–28} Space plasmas are frequently found in the latter collisionless regime, and fusion plasmas can also reach this regime.^{29–32}

In this paper, we revisit the linear kinetic dispersion relation for collisionless tearing instability of a Harris current sheet.^{33} Recent particle-in-cell (PIC) simulations, together with numerical solutions of the linearized Vlasov equation,^{34} have revealed the surprising result that the spectrum of oblique modes is much more narrow than standard analytic theories predict.^{24,25} Here, $\theta =arctan(ky/kz)$ is the angle of obliquity and **k** the wavevector; see Fig. 1. The resonant surface locations $xs=\lambda arctanh(\mu )$, where $\mu \u2261\u2009tan\u2009\theta Boy/B\xafoz$, locate the chain of magnetic islands associated with angle *θ*. This narrower spectrum has significant consequences for the subsequent nonlinear reconnection dynamics. Flux ropes on adjacent resonant surfaces can overlap, which leads to stochastic magnetic field generation and even turbulence.^{35,36} The spectrum of unstable tearing modes determines the plasma volume that is ultimately susceptible to this form of turbulence. It also has consequences for theories of particle acceleration by contracting magnetic islands, which relies on a broad spectrum of unstable modes (i.e., volume filling of magnetic islands).^{37,38} Resonant surfaces unstable to linear tearing modes also become the sites of electron layers in the nonlinear regime,^{39} which can lead to multiple electron layers in domains long enough to support multiple resonant surfaces. Understanding the spectrum of linear tearing modes is critical.

This work focuses on understanding the stabilization of oblique modes. There are many theories of linear collisionless tearing modes in the literature,^{24,25,40–49} but none of them captures this effect. We focus on comparing with Drake and Lee (DL)^{24} and Galeev *et al.*,^{25} as these are the most commonly used. In a kinetic plasma, the linear response is determined by orbit integrals, leading to an eigenvalue problem with complex integro-differential operators. All of the theories make assumptions to reduce the complexity of these equations, including various combinations of a small gyroradius expansion, the constant-*ψ* approximation, neglect of electrostatic terms, etc. Each of these theories also assumes a scale separation between an “inner” region where reconnection occurs and an “outer” (ideal) region.

To identify the approximation responsible for missing the effect, we apply a series of three successive approximations and discuss the influence of each on the growth rate. The most general description is a normal mode analysis of the Vlasov-Maxwell integro-differential (VMID) equations that includes the full orbit integrals. This full system of equations has been solved numerically for the Harris sheet geometry using both Hermite^{50} and finite-element expansions^{9,51} of the eigenfunction. In this paper, the results from the finite-element code^{51} are used to evaluate and refine simplified theories. As a first step, the full system of equations is simplified by taking the small gyroradius limit. The next step exploits a scale separation that allows one to eliminate electrostatic terms from the system of equations. We find that each set of equations captures the narrower mode spectrum, in qualitative agreement with the previous kinetic simulations and the full numerical solutions of the tearing eigenmode equations. Quantitative differences do arise at each level of approximation.

We find that the physical process responsible for the stabilization of oblique modes is the uniform diamagnetic drift in the Harris equilibrium. The drift broadens the inner layer of oblique modes. At large angles, the inner region becomes broad enough that the scale separation that boundary layer theories are predicted on breaks down. In this regime, both kinetic simulations and VMID calculations indicate that tearing modes are damped. In this work, we demonstrate that this stabilization effect can be described from the numerical solution of a single second order differential equation. A comparison with previous boundary layer theories is provided, which details asymptotic limits and where they agree and disagree with the more complete theory.^{24,25,45,49} A simple analytic formula for the cutoff angle at which stabilization occurs is provided.

It is noteworthy that previous work on other models, or for other equilibria, have not seen the same stabilization of oblique modes. For example, modes were found to span the current sheet in resistive MHD.^{52} Simulations and collisionless kinetic theory, similar to what are discussed here, were found to agree well for a force-free equilibrium.^{39} Similarly, two-fluid theory was found to accurately capture the dispersion of oblique modes in this case.^{53} These results are consistent with what is presented here because the force-free equilibrium does not have a diamagnetic drift. We also note that Cowley *et al.*^{54} have shown that temperature gradients can stabilize tearing modes, but the density gradient did not cause stabilization in that analysis and it produced only logarithmic corrections to the growth rates calculated without a density gradient.^{29} There is no temperature gradient in the Harris equilibrium.

The remainder of this paper is organized as follows. The three sets of equations obtained from successive approximations of the Vlasov-Maxwell system are derived in Sec. II. Numerical solutions are provided in Sec. III. Section IV gives a comparison of these results with previous theories. Section V shows an analysis of how the boundary layer theories break down at oblique angles and provides a simple analytic estimate for the cutoff angle.

## II. BASIC EQUATIONS

### A. Vlasov-Maxwell

This section describes the Vlasov-Maxwell integro-differential equations in which the full orbit integrals are solved numerically. We will refer to this as the VMID model. It starts from the linearized Vlasov equation

where $d/dt=\u2202/\u2202t+v\xb7\u2207+(qs/ms)(Eo+v\xd7Bo/c)\xb7\u2207v$ is the convective derivative. We will later focus on the Harris equilibrium, but in this section, the requisite is only that the equilibrium magnetic field can be written in the form $Bo=Byo(x)y\u0302+Bzo(x)z\u0302$ and that the lowest-order distribution functions are Maxwellians

where the flow velocity is $Us=Us,y(x)y\u0302+Us,z(x)z\u0302$. We write $\delta E$ and $\delta B$ in terms of scalar and vector potentials

and apply a normal mode analysis of the form

Inserting Eq. (2) into (1) and integrating along the characteristics, which are the unperturbed particle orbits, provides^{50}

where

$\tau =t\u2212t\u2032,\u2009\delta \varphi \u0303=\delta \varphi (k,x\u2032,t)$ and $\delta A\u0303=\delta A(k,x\u2032,t)$. The single particle characteristics ($x\u2212x\u2032$) are quite complicated in the sheared electric and magnetic fields of a general neutral sheet configuration. In the VMID model, *S* was computed numerically from Eq. (5) for the exact single particle characteristics using the orbit integration technique described in Refs. 50 and 51.

Equations (4) and (5) are used to calculate the perturbed charge density ($\delta \rho =\u2211sqs\delta ns$, where $\delta ns=\u222bd3v\u2009\delta fs$) and current density ($\delta J=\u2211sqs\delta \Gamma s$, where $\delta \Gamma s=\u222bd3v\u2009v\delta fs$), which are then used in Maxwell's equations in the Lorenz gauge

^{51}The boundary conditions on $\delta \varphi $ and $\delta A$ are that they asymptote to a constant value (zero) as $x=\xb1\u221e$.

### B. Small gyroradius limit

As a next level of approximation, gyroradius effects are neglected ($\rho s\u21920$). In this case, the single particle characteristics are $x\u2212x\u2032=v\u2225\tau b\u0302+O(\rho s)$, where $b\u0302=Bo/|Bo|$, and Eq. (5) reduces to

assuming $\u2111{\omega}>0$. To the lowest-order in this gyroradius expansion, Eq. (4) can then be written as

Since the phase speed of the waves of interest are much smaller than *c*, we drop the displacement current terms in Maxwell's equations, so $\u22072\delta \varphi =\u22124\pi \delta \rho $ and $\u22072\delta A=\u2212(4\pi /c)\delta J$.

Evaluating the velocity-space integrals for the perturbed density and putting the result into Gauss's law gives

in which $\lambda Ds2=\u2211s(4\pi qs2nso)/Ts$ is the square of the Debye length for species *s*, and $\lambda D\u22122=\u2211s\lambda Ds\u22122$ gives the total Debye length

The coordinate rotation $(x,y,z)\u2192(x,b,\eta )$ has been applied, where $b\u0302=byy\u0302+bzz\u0302$ is parallel to the equilibrium magnetic field and $\eta \u0302=x\u0302\xd7b\u0302=\u2212bzy\u0302+byz\u0302$. The plasma dispersion function^{55}

has also been used where $Z(w)\u2261Z0(w)$. Evaluating the integrals for the perturbed current and putting the result into Ampère's law gives

where $\omega ps=4\pi qs2nso/ms$ is the plasma frequency of species *s*. Also, note that $\delta Bx=i(k\u2225\delta A\eta \u2212k\eta \delta A\u2225)$ and $\delta E\u2225=\u2212ik\u2225\delta \varphi +i(\omega /c)\delta A\u2225$.

Explicit evolution equations for $\delta Ay$ and $\delta Az$ can be obtained by projecting Eq. (12) along $y\u0302$ and $z\u0302$. However, more compact expressions can be obtained by projecting Eq. (12) in the direction along $Bo$ *at a resonant surface*

as well as the direction that is perpendicular to both this and $x\u0302$, which is simply $k\u0302:b\u0302s\xd7k\u0302=x\u0302$. Thus, $(x\u0302,b\u0302s,k\u0302)$ provides a convenient orthogonal coordinate system. Applying this to Eq. (12) gives

and

in which

$Vbs=b\u0302s\xb7V$ and $Vk=k\u0302\xb7V$. Here, $\delta Ak=k\xb7\delta A$

and

### C. Decoupled magnetic field approximation

Although $\delta E\u2225$ in Eq. (14) depends on $\delta \varphi ,\u2009\delta A\u0302$ and $\delta Ak$, consideration of two-scale features of a tearing layer can be used to justify an approximate decoupling of $\delta A\u0302$ from $\delta \varphi $ and $\delta Ak$. In the outer region, far from a resonant surface, magnetic flux is frozen into the plasma, so the parallel electric field vanishes $\delta E\u2225\u21920$. In this ideal MHD-like outer region, the right side of Eq. (14) is negligible and $\delta A\u0302$ is the only remaining dependent variable. In the inner region near a resonant surface, $k\u2225\u226ak$, which implies

and

so,

If we also note that the first term on the right side of Eq. (14) is small compared to the other two, we are left with

to describe the inner region.

Although Eq. (22) was derived for the inner region, it is expected to hold in the outer region as well. This is because the left side is the ideal MHD outer region equation, and in the outer region, the right side is negligibly small since it is proportional to $\omega 2/(kvTs)2\u226a1$. The right side contributes only in the inner region where $k\u2225$ is small. Thus, Eq. (22) is an approximate equation spanning the layer that is decoupled from $\delta Ak$ and $\delta \varphi $. The validity of this scale separation will be evaluated in Sec. V.

### D. Harris equilibrium

We utilize the following properties of the Harris equilibrium^{33}

where an additional uniform background density (*n _{b}*) for each species is included, $n(x)=n\xaf(x)+nb$. Noting that $Us,\eta =\u2212bzUs$, Eq. (16) reduces to

and Ampère's law gives

where “prime” denotes a derivative. Taking a derivative of Eq. (23b) and utilizing Eqs. (23c) and (25) give

Substituting this into Eq. (24) provides

where $F=k\xb7Bo$. The individual components of interest are then

and $Vk=tan\u2009\theta Vbs$.

## III. NUMERICAL SOLUTIONS

In this section, the dispersion relation computed from the VMID model of Sec. II A is compared with that computed from the approximate layer equation of Sec. II C [Eq. (22)]. The VMID computations were calculated using the method and the code described in Refs. 50 and 51. Solutions of Eq. (22) were obtained by numerically integrating from $x=\xb1\u221e\u2192xs$ and determining the dispersion relation from the requirement that $\delta A\u0302\u2032(xs+)=\u2212\delta A\u0302\u2032(xs\u2212)$.

### A. Growth rate profiles

The VMID solutions shown in Fig. 2(c) demonstrate that at a realistic mass ratio ($mi/me=1836$), the growth rate is a maximum at $\theta =0\xb0$, corresponding to the usual parallel mode at the center of the current sheet, and it monotonically decreases with increasing angle of obliquity. The growth rate falls off slowly with the angle until near $\theta =20\xb0$, where it rapidly decreases, and quickly reaches stability for larger angles. These general features of the growth rate profiles hold at each of the current sheath widths considered ($\lambda =\rho i,\u20092\rho i$ and $5\rho i$), but details of the slope of the fall off differ in each case. The real frequency of the instability is to a good approximation linearly proportional to *θ* over the entire range of unstable angles at each sheet width.

Figure 2(b) shows the results of analogous simulations, but at a lower mass ratio ($mi/me=100$). Similar trends in the growth rate as for the proton-electron mass ratio are observed, including the maximum growth rate at $\theta =0\xb0$ and a monotonic decrease of the growth rate. The primary differences in this comparison are that the unstable range of wavenumber extends further out on the current sheet (to $\theta \u224330\xb0$) and that the values of the growth rate in units of Ω_{ci} take larger values at a lower mass ratio. Here, $\Omega ci=eB\xafoz/(mic)$ is the ion gyrofrequency based on the asymptotic magnetic field strength $B\xafoz$. In the latter comparison, an additional consideration is the mass dependence of the growth rate units being displayed. Accounting for this, $1836/100=18.36$, demonstrates that in dimensional units ($s\u22121$), the growth rate is higher at a larger mass ratio.

Figure 2(a) shows that at the mass ratio of unity $mi/me=1$, qualitative differences in the growth rate profiles are observed. The most apparent is that the growth rate is positive and near constant over the entire range of angles. This suggests that at unity mass ratio, linear tearing modes would be excited throughout the entire current sheet; a result that has also been observed in kinetic simulations.^{56} Another apparent difference is that for larger current sheet widths ($\lambda =5\rho i$ and, to a lesser extent, $\lambda =2\rho i$), the most unstable mode is oblique.

Each of the panels in Fig. 2 also shows the predictions of the approximate layer equation from Eq. (22). A comparison with the Vlasov simulation curves shows that, despite the severity of some the approximations made in Secs. II B and II C, these solutions capture key features of the growth rate profiles and also give a fair quantitative prediction. The key agreement is that at large mass ratios ($mi/me=100$ and 1836), the range of unstable angles is well approximated by this solution. From the perspective of understanding how linear tearing modes will contribute to the overall reconnection dynamics, this is a key feature because it determines the locations in the current layer that are expected to be unstable. It is also observed that at large mass ratios, Eq. (22) accurately predicts the growth rate of the central mode at $\theta =0\xb0$ and generally predicts a flat profile. One distinction observed at $\lambda =5\rho i$ and $mi/me=1836$ is that Eq. (22) predicts that the most unstable mode is oblique (near the cutoff), whereas the Vlasov solutions predict a monotonically decreasing profile. Finally, much poorer agreement between Eq. (22) and the Vlasov solutions is observed at unity mass ratio. With the exception of the larger current sheet widths and angles near the center of the current sheet, the agreement is poor. This is expected because, for a pair plasma, the electron gyroradius becomes comparable to the gradient scale of the current sheet.

The observed stabilization at sufficiently oblique angles is associated with the diamagnetic drift $Us$. This is confirmed by the dashed lines in Fig. 2, which show solutions of Eq. (22) with $Us=0$. These show that in the absence of the diamagnetic drift, the theory predicts positive and nearly constant growth rates across the entire current sheet. Recall from Eq. (23d) that the diamagnetic drift in the Harris equilibrium is associated with the density gradient, since the temperature is uniform. Cowley, Kulsrud and Hahm previously studied the influence of temperature gradients.^{54} It is also interesting to note that for force-free equilibria, which do not have diamagnetic drifts, tearing modes have been observed to span the current sheet.^{39,53} The results of Eq. (22) in the absence of diamagnetic drift agree well with the predictions of the standard kinetic tearing mode theory of Drake and Lee.^{24} In Sec. V, it is shown that the reason for the absence of the drift stabilization of oblique modes in the theory is a broadening of the inner layer at oblique angles, which causes the boundary layer approximation that the theory utilizes to break down. Because the full numerical solution of Eq. (22) accurately captures the essential feature of the unstable region of the current layer, the remainder of this work will focus on this equation.

### B. Vector potential and phase shifts

Figure 3(a) shows profiles of the real and imaginary components of the vector potential $\delta A\u0302(x)$ across the current sheet for modes centered at resonant surfaces corresponding to $\theta =0\xb0,\u20095\xb0,\u200910\xb0$, and $15\xb0$. This solution assumed a proton-electron mass ratio $mi/me=1836$, and the other parameters chosen were the same as in Fig. 2(c). The eigenfunction corresponding to the mode at the central resonant surface ($\theta =0\xb0$) is purely real and is symmetric about *x* = 0. This is the common MHD solution of the outer region familiar from the seminal Furth, Killeen, and Rosenbluth (FKR)^{1} analysis. Oblique modes are centered at resonant surfaces off-center of the current sheet, and the most pronounced feature of the eigenfunctions is the formation of an asymmetry about the resonant surface. This asymmetry was also observed in our previous reduced MHD analysis of oblique tearing modes (see Fig. 7 of Ref. 52), and is associated with the ideal MHD outer region for off-center modes.

A distinguishing feature of the kinetic solutions, compared to MHD solutions, is that the eigenfunctions of the oblique modes are complex. Figure 3 shows that all off-center modes ($\theta >0$) have a finite $\u2111{\delta A\u0302(x)}$, and that their profile is also asymmetric with respect to the resonant surface. The values of this component are negative (the plot shows $\u2111{\u2212\delta A\u0302}$). Figure 3(b) shows that if the complex vector potential is represented as $\delta A\u0302(x)=|\delta A\u0302(x)|\u2009exp\u2009[i\varphi (x)]$, there is a phase shift of the eigenfunction across the current sheet for oblique modes. The phase factor asymptotes to a constant value far from the center of the current sheet and vanishes at the resonant surfaces. This phase shift is observed to increase with increasing angle of obliquity. This phase factor is also a distinctly kinetic feature of the eigenfunctions.

Boundary layer theory assumes that the outer region is described by ideal MHD and does not account for a complex vector potential. Thus, one may expect the larger imaginary component of $\delta A\u0302$ at oblique angles to be associated with the observed breakdown of the analytic theories. However, we find that the dominant effect is a breakdown of the boundary layer scale separation, rather than the complex eigenfunction. This aspect will be discussed further in Sec. V.

## IV. RELATION TO PREVIOUS THEORIES

### A. Drake and Lee

The seminal kinetic tearing mode theory of Drake and Lee (DL)^{24} applies a boundary-layer analysis that considers only electrons in the inner region. A dispersion relation is obtained by matching the inner layer solution with an outer layer solution through a matching parameter called the tearing stability index

which is provided by an independent ideal MHD solution. Generalizing the inner region equation to include ions gives

Equation (30) agrees with the inner region of Eq. (22) in the limit that $ws\u2192\zeta s$, which implies $\omega /k\u2225\u226bUs,\u2225$. Note that $Z1(\zeta s)=\u2212Z\u2032(\zeta s)/2=1+\zeta sZ(\zeta s)$.

The DL dispersion relation follows from applying the constant-*ψ* approximation. In this approximation, $\delta A\u0302$ on the right side of Eq. (30) is assumed constant across the inner layer. Dividing by this and integrating across the layer give $\u222bdx\u2009\delta A\u0302\u2033/\delta A\u0302(xs)\u2243\Delta \u2032$, and

The full solution of $k\u2225$ for the Harris equilibrium described in Sec. II D is

To integrate the right side of Eq. (31), an additional assumption is applied, where the space-dependent functions $k\u2225$ and *ω _{ps}* are expanded about the resonant surface location: $k\u2225\u2243k(x\u2212xs)/ls$, and $\omega ps=\omega ps(xs)$, where

Figure 4 shows a plot of the full space-dependent solution of Eq. (32) along with the linear expansion at several resonant surface locations. This shows that the width over which the linear expansion is accurate depends on the resonant surface location.

In Ref. 24, the integral in Eq. (32) is evaluated by applying the above approximations. The result is

With this, Eq. (31) provides the dispersion relation

where $ds\u2261c/\omega ps$ is the skin depth associated with species *s*. Rearranging, the real frequency and growth rate are

Applying the expressions for a Harris sheet with a background density from Sec. II D provides $\omega =\omega r+i\gamma $, where

is the real frequency of the tearing mode, and

is the growth rate. Here, $n\u0302b\u2261nb/n\xafo$ and $\alpha \u2261Teme/(Timi)$. For the plots in Fig. 2, which show $\omega /\Omega ci$, it is convenient to note that $kyUe=\Omega cik\lambda \u2009sin\u2009\theta (Te/Ti)(\rho i/\lambda )2$ and $vTede2|xs=\Omega ci\rho i3\alpha (1+Te/Ti)/(1\u2212\mu 2+n\u0302b)$.

The tearing stability index $\Delta \u2032$ must be supplied from a model for the outer region of the boundary layer. This is typically assumed to obey ideal MHD (FKR)^{1}

Recall that $Vbs=F\u2033/F$ is provided in Eq. (28). Figure 5 shows the solution of Eq. (39) throughout the domain of possible resonant surfaces. One prominent feature is that $\delta A\u0302(x)$ is an even function of *x* only for *θ* = 0. For off center modes, the characteristic peak is much higher on the low magnetic field side of the current sheet ($x<xs$) than on the high magnetic field side ($x>xs$). A model for $\Delta \u2032$ based on Eq. (39), which considers oblique modes, was developed in Ref. 52

This was validated against numerical solutions of Eq. (39) over a broad range of conditions (see Fig. 3 of Ref. 52), and a similar validation is shown in Fig. 6(b). It was used to evaluate $\Delta \u2032$ for the DL theory curves shown in Fig. 2.

Figure 2 shows that the solution of the DL theory, obtained this way, accurately predicts the growth rate of parallel ($\theta =0\xb0$) modes; the only exception being thin current sheets at unity mass ratio. However, the figure also shows that theory does not accurately capture oblique modes, particularly the stabilization. It is interesting to notice that Eq. (22) does capture the range of unstable modes, and that it is very closely related to the fundamental equation that the DL boundary layer theory is based on. The only difference is that the argument of the *Z*-function in Eq. (22) is *w _{s}*, whereas it is

*ζ*in DL theory. We have also numerically solved a version of Eq. (22) with $Z(ws)$ replaced by $Z(\zeta s)$, and have found that the solutions are very close (they also capture the stabilization of oblique modes). Thus, the discrepancy with DL theory cannot be explained by this difference.

_{s}### B. Galeev *et al.*

Previous researchers^{25,57,58} have studied oblique collisionless tearing modes in the context of Earth's magnetopause and have obtained a different result than Drake and Lee. These papers keep a first-order finite gyroradius correction, but the small gyroradius limit of their equations can be derived from Eq. (8) by applying the following assumptions: (1) $k\u2225\u226ak$ and (2) $Us/vTs\u226a1$. These imply

so, Eq. (9) reduces to

For Ampère's equation, it is convenient to first rearrange Eq. (12) to the form

The $k\u2225\u226ak$ assumption implies $\delta A\u2243\delta A\u0302b\u0302s$ and $b\u0302s\u2243b\u0302\u2243z\u0302$. Applying these along with assumption (2), the $b\u0302s$ projection of Eq. (43) reduces to

in which

is what Galeev calls the adiabatic interaction term.

Equations (42) and (44) are the two basic equations that this line of previous collisionless oblique tearing mode research^{25,57,58} worked from.^{59} This set of two coupled equations is substantially different from any of the approximate equations derived in Sec. II. The most significant difference is in the description of the outer region. For an ideal MHD outer region, the $Vad$ is replaced by $Vbs$, which reduces to Eq. (28) for Harris equilibrium. However, the adiabatic response used in Eq. (44) is

for Harris equilibrium. This agrees with the ideal MHD response, Eq. (28), only in the limit of a standard symmetric layer $(xs\u21920)$, i.e., the parallel mode. The cause of this disagreement is that the $k\u2225\u226ak$ assumption holds only near a resonant surface. This leads to substantial errors in the outer region, which also affects the dispersion relation.

To illustrate that the primary difference between Galeev *et al.* and our approach is the outer region solution, we first apply the $i\delta E\u2225\u2243\u2212\omega \delta A\u0302/c$ approximation from Eq. (19), which was validated in Sec. III, to Eq. (44). This gives

which is the same as the $Us,\u2225\u226a\omega /k\u2225vTs$ limit of Eq. (22) except in the outer region (left side) where $Vad$ replaces *V _{bs}*.

Figure 6(a) shows solutions of the outer region of the Galeev equation

in comparison to ideal MHD from Eq. (39). For the parallel mode, the solutions are identical. However, the solutions are dramatically different for finite *θ*, showing the inconsistency of the Galeev outer region equation with ideal MHD. To further emphasize this point, Fig. 6(b) shows a comparison between $\Delta \u2032$ obtained from the numerical solutions of Eqs. (48) and (39), and the model solutions from Eq. (40) and Galeev's approximate solution^{60}

where

### C. Electrostatic effects

Previous work has considered the influence of electrostatic effects, showing that the ion acoustic wave can couple with the tearing mode, causing stabilization.^{45–48} Electrostatic effects were eliminated from our analysis in Sec. II C based on an argument of scale separation between the inner tearing layer and the outer ideal region. These theories consider cases where the inner layer can be sufficiently broad that this scale separation is not satisfied. Indeed, we will see in Sec. V that the inner layer does broaden substantially for oblique modes, calling into question the neglect of electrostatic effects. It is therefore prudent to compare with these theories.

where $k\u2225\u2032=k/ls,\u2009\rho s=cs/\Omega ci$, and $\omega o=\u211c{\omega DL}$. Here, $\omega DL$ is the DL dispersion relation from Eq. (35). The second term on the right side of Eq. (51) is the correction due to electrostatic effects. The figure shows that this theory does predict a smaller growth rate for oblique modes than DL. However, the predicted values do not agree well with the VMID solutions, or capture the trend of a sharp cutoff in the growth rate.

Electrostatic coupling does not appear to be the cause of stabilization of the oblique modes since the approximate theory neglecting electrostatic terms, Eq. (22), captures the effect, while boundary layer theories including electrostatics do not. Although the scale separation used to justify the neglect of electrostatic effects is similar to that used to justify the boundary layer theory, it is not the same. Namely, in Eq. (19), the approximation made was $k\u2225\delta \varphi \u226a\omega ck\eta k\delta A\u0302$. Although wide tearing layers at oblique angles access a range of values far from the current sheet, so $k\u2225/k$ may not be small, the real frequency of oblique modes also increases approximately as $kyUe$. Thus, these effects compete in such a way that the electrostatic terms may remain small in comparison to the vector potential terms. Indeed, the electrostatic component of the eigenfunctions from the VMID solutions have been shown to be small in comparison to the vector potential contributions even at oblique angles (see Fig. 8 of Ref. 9). This reference also showed that the ratio of the electrostatic to electromagnetic contribution to the eigenfunction scales linearly with $\omega pe/\omega ce$. Thus, one can alter the magnitude of the electrostatic terms with this parameter. Nevertheless, the computed growth rate remained independent of this parameter.^{9}

### D. Nonconstant-*ψ* approximation

The DL theory applies the constant-*ψ* approximation to obtain an analytic dispersion relation. This works best if the inner layer is very thin because $\delta A\u0302$ does not vary significantly over sufficiently short distances. For thicker inner layers, the approximation begins to break down. Methods have been developed to extend the theory, typically by keeping a linear term in the expansion of $\delta A\u0302$ near the resonant surface. These theories are called large $\Delta \u2032$, or nonconstant-*ψ* approximations. It was observed in the MHD theory that a nonconstant-*ψ* approximation was required to accurately model oblique modes.^{52} Section V will show that the inner layer of the collisionless problem becomes thicker at oblique angles. It is thus natural to investigate if nonconstant-*ψ* effects are responsible for the observed stabilization.

A large $\Delta \u2032$ theory of kinetic collisionless tearing modes has been computed by Mahajan *et al.*,^{49} which provides the dispersion relation

Here, $\gamma DL=\u2111{\omega DL},\u2009\omega o=\u211c{\omega DL}$, and $vA=B\xafoz/4\pi n\xafomi$ is the Alfvén speed. Figure 7 shows a comparison of the results of Eq. (52) with the VMID solutions. The nonconstant-*ψ* modification does reduce the growth rate substantially at oblique angles. However, it still predicts a much broader spectrum of unstable modes than either VMID or Eq. (22). Thus, the large $\Delta \u2032$ extension of DL theory does not appear to be able to explain the observed stabilization.

## V. ANALYSIS OF BOUNDARY LAYER THEORY

### A. Gyroradius effects

A small gyroradius expansion was used in Sec. II B to reduce the full particle orbits from the Vlasov description in Eqs. (5)–(7). This is strictly valid only in the limit that both the electron and ion gyroradii are smaller than other length scales of relevance to the tearing mode. Figures 8(d)–8(f) show that this approximation is clearly violated under the conditions explored in this paper. The largest gradient scale is the current sheet width, which for the parameters chosen in Fig. 2, ranges from $\lambda =\rho i$ to $5\rho i$. Thus, ions are not strongly magnetized even on the scale of the shear of the current sheet. Electrons, however, are strongly magnetized at this scale for the large mass ratio cases, since $\rho e/\lambda \u223cme/mi\u226a1$. However, use of the approximation at the scale of the inner reconnection layer faces a more stringent requirement. Figures 8(d)–8(f) plot an estimate of the ratio of the ion inner tearing layer width to the ion gyroradius, showing that in all cases considered, the inner tearing layer width is smaller than the gyroradius. Since $\delta i/\rho i$ is independent of mass, the same result applies for electrons. Here, the inner tearing layer width was estimated as in Drake and Lee, which provides $\delta s\u2243|\omega |ls/(kvTs)$, where *ω* was computed from Eqs. (37) and (38).^{24} Although Fig. 2 has shown that this does not accurately model the growth rate of oblique modes, it does accurately model the real mode frequency. Since the inner layer width estimate is dominated by the larger real component of the frequency, we expect that it remains a reasonable estimate at oblique angles.

Although the small gyroradius expansion is not expected to be valid under these conditions, a comparison with the full VMID solutions in Fig. 2 shows that the predictions resulting from it, namely Eq. (22), capture the essential features of the growth rate profiles. The comparison reveals quantitative differences in the predicted growth rate of oblique modes, but it is quantitatively accurate for the parallel ($\theta =0\xb0$) mode (except at unity mass ratio), as well as the angle of obliquity at which stabilization occurs. The full orbit integrals are complicated, and any reduced analytic model must treat the orbits in an approximate fashion. The comparison in Fig. 2 provides confidence that the small gyroradius expansion captures the essential features of interest from linear theory.

It may seem surprising that the growth rate derived from the small gyroradius expansion captures the essential features of the full VMID solution under conditions where the expansion is not expected to be valid. In this regard, one should notice that the linear growth rate depends on the total perturbed current within the resonance layer. Previous work has shown that the electron finite Larmor radius (FLR) broadens this layer, but does not significantly influence the total perturbed current.^{11,12} Specifically, Fig. 1 of Ref. 11 shows that the perturbed current from the linear VMID calculation is smeared over a scale of a few *ρ _{e}*; it is not possible to form a narrower current channel. Although FLR effects do not seem to significantly influence the linear growth rate, they have bigger implications for the saturation. Exponential growth, as described by linear theory, is expected to saturate when the island width becomes comparable to the resonant layer thickness. Since FLR effects broaden the resonant layer, in comparison to the classical theoretical predictions, this allows modes to grow linearly to much larger amplitudes than would otherwise be expected. This issue is particularly relevant to magnetospheric layers, and has been discussed by Quest and Coroniti.

^{30,31}

### B. Inner layer width

Boundary layer theory is predicted on the assumption of an asymptotic scale separation between an “outer region” described by ideal MHD and an “inner region” where two-fluid and kinetic physics leads to violations of the frozen flux condition. Such a two-scale approximation was used not only in the theories of Sec. IV, but also in Sec. II B, which led to the decoupling of the electric and magnetic fields allowing one to describe the tearing mode from the single differential equation, Eq. (22). The possible absence of this scale separation is what led Coppi to suggest that a third electrostatic layer must be accounted for in some circumstances, as discussed in Sec. IV C.

Figures 8(a)–8(c) show the ratio of the estimated inner tearing layer width (*δ _{i}*) to the width of the current sheet (

*λ*), which characterizes the scale of the ideal MHD outer region. At unity mass ratio, the scales are not well separated. We expect that this is the cause of the poor agreement between the VMID solutions and the solutions of Eq. (22) under these conditions. But, for large mass ratio cases $\delta i/\lambda \u226a1$ near the center of the current sheet (small

*θ*), the scale separation is violated at highly oblique angles. This calls into question the validity of this assumed scale separation at oblique angles. Despite the scale separation becoming questionable at oblique angles, the fact that Eq. (22) captures the stabilization suggests that the influence of electrostatic terms is not the cause of this effect; see Sec. IV C. Regarding the subsequent application of the boundary layer approximation in developing an analytic theory, specifically the DL theory from Sec. IV A, the approximation can be tested in more detail.

To assess the boundary layer theories, we compare the left side of Eq. (22), which has a characteristic scale of *λ*, with the inner region which is described by the right side

(Here, the version of Eq. (22) with *w _{s}* replaced by

*ζ*in the

_{s}*Z*-function is used to make connection with the DL theory.) In boundary layer theory, the spatially dependent terms are expanded locally about the resonant surface: $k\u2225\u2243k(x\u2212xs)/ls$ and $\omega ps\u2243\omega ps(xs)$, leading to an approximate expression describing the inner region

Figures 9(a) and 9(b) show a comparison of the solutions of Eqs. (53) and (54) for the same parameters as the $\lambda =2\rho i$ case of Fig. 2(c). Here, the mode frequency is computed using the DL solutions from Eqs. (37) and (38). This reiterates the breakdown of a scale separation at oblique angles, which was estimated in Fig. 8. As *θ* approaches several degrees near the value where the growth rate rapidly diminishes, the scale associated with the inner region $[T(x)]$ is observed to grow to a comparable value with the scale associated with the outer region (*λ*). This explains why boundary layer theory fails to accurately model the numerical solutions of Eq. (22) at sufficiently oblique angles. The numerical solutions reveal that the growth rate rapidly diminishes under the same conditions; see Fig. 2.

### C. Phase and the outer region

Figure 9 answers another inconsistency between the boundary layer and numerical solutions: Why there is a prominent complex (phase) component of $\delta A\u0302$ for oblique modes ($\theta \u22600$) in the numerical solutions, while this does not arise in the boundary layer theory.

First, we observe that *T*(*x*) is purely real for the parallel mode ($\theta =0\xb0$) (as a result of *ω* being purely imaginary). Thus, no phase component is expected in this case, nor is one observed in the numerical solutions; see Fig. 3. The figure also shows that $T(x=\xb1\u221e)\u21920$ and that *T*(*x*) is localized if *θ* = 0. This case is amenable to a boundary layer analysis because the inner region decays to zero before reaching the scale of the outer region, and the inner and outer regions are both purely real and can be matched asymptotically. The boundary layer theory agrees with the numerical solutions in this case.

Oblique modes have fundamentally different properties. For any finite angle, *T*(*x*) becomes complex (has both real and imaginary components) implying that $\delta A\u0302$ will be complex (i.e., have a finite phase). This is consistent with what is observed in the numerical solutions at finite angles; see Fig. 3. In addition, Fig. 9 shows that at finite angles $T(x=\u2212\u221e)\u22600$, if $n\u0302b\u22600$. Since part of *T*(*x*) extends far from the resonant surfaces, this implies that at finite angles, the kinetic theory does not asymptote to the ideal MHD outer region. Instead, one should move the components of *T*(*x*) extending to asymptotically far distances $T(x\xb1\u221e)$ to the left side of Eq. (22) (the outer region). In fact, if one does not do this, the finite value of $T(x=\u2212\u221e)$ implies that the integral $\u222bdxT(x)$, which is used in the constant-*ψ* approximation, actually diverges. This divergence happens to not arise in DL theory when integrating $\u222bdxTDL(x)$ as a fortuitous consequence of the local expansion of *ω _{ps}* and $k\u2225$ (which is not valid in this case). The local expansion forces $TDL(x=\xb1\u221e)\u21920$ even at finite $n\u0302b$; see Figs. 9(a) and 9(b).

### D. Role of background density

Since *T*(*x*) extends to asymptotically far distances only when $nb\u22600$, it is interesting to compare with the situation of no background density *n _{b}* = 0. This allows one to isolate the role of the imaginary component of $\delta A\u0302$ to determine if it influences the growth rate, and in particular, if it is associated with the damping observed at oblique angles. If

*T*(

*x*) is sufficiently local near the resonant surface, one expects that $\u2111{\delta A\u0302}$ will be small compared to the real component. Figures 9(c) and 9(d) show that

*T*(

*x*) and $TDL(x)$ agree quite well in this case. Figures 10(c) and 10(d) show solutions of Eq. (22) with $n\u0302b=0$, confirming that $\u2111{\delta A\u0302}$ is small compared to the real component when

*θ*is small (see also Fig. 3), but that it grows comparable to the real component for $x>xs$ at oblique angles. Figures 10(a) and 10(b) show solutions of the dispersion relation computed from Eq. (22) with $n\u0302b=0$. These are found to damp at a similar angle as the case with $n\u0302b=0.3$ from Fig. 2(c). This provides further evidence that the damping is caused by the broadening of the inner region, not simply the finite phase component.

### E. Role of guide field

Figure 11 shows that the narrower spectrum of unstable modes predicted by Eq. (22), in comparison to DL, persists over a wide range of guide field strength. For this range of parameters, the cutoff angle is also independent of $\lambda /\rho i$. Each of these observations is consistent with the notion that stabilization is associated with the width of the inner layer exceeding the shear scale *λ*. This is further quantified in Sec. V G. As the ratio of the guide field to the asymptotic field strength, $Boy/B\xafoz$, exceeds unity, the range of unstable angles narrows substantially simply because the field is predominantly in the guide field direction. This behavior is also shown in Fig. 12, which displays the cutoff angle (where *γ* = 0) as a function of $Boy/B\xafoz$ from each theory.

### F. Role of temperature ratio

Figure 13 shows that the cutoff angle computed from Eq. (22) depends significantly on the temperature ratio $Ti/Te$. In contrast, the cutoff angle is independent of the temperature ratio in the DL model. At $Ti/Te\u223c1$, the spectrum of unstable modes from Eq. (22) is significantly more narrow than from the DL. However, as $Ti/Te$ increases, the spectrum predicted from Eq. (22) broadens substantially. At a large temperature ratio ($Ti/Te=100$), the two results agree very well; see Figs. 13(c) and 13(d). The same trend is captured in Fig. 12(b), which shows the cutoff angle as a function of temperature ratio for three values of the ratio $\lambda /\rho i$. The Drake and Lee limit is obtained at large temperature ratios. Since the scale separation between the inner region and the *λ* scale is broader at large $Ti/Te$, this also provides further evidence associating the cutoff in the growth rate with the broadening of the inner region. This will be further quantified in Sec. V G. This temperature ratio dependence is directly related to the drift stabilization. For a Harris sheet, the temperature ratio sets the ratio of flow velocities $Ui/Ue=\u2212Ti/Te$, while the ion velocity is fixed by the current sheet thickness $Ui/vTi=\rho i/\lambda $. Thus, changing the temperature ratio effectively changes the electron drift speed $Ue=\u2212UiTe/Ti$. For $Ti/Te\u226b1$, *U _{e}* is very small, and one returns the results of the DL theory valid in the limit $Ue\u21920$. Finally, Figs. 12(b) and 13 show that the cutoff angle is essentially independent of $\lambda /\rho i$, when

*T*>

_{i}*T*, but that a significant dependence on this ratio arises when

_{e}*T*<

_{i}*T*.

_{e}### G. Estimate of instability threshold

An important contribution that linear tearing theory can provide to the overall understanding of reconnection dynamics is the volume of space susceptible to tearing (or plasmoid formation), as this has significant consequences for subsequent heating processes as well as the volume susceptible to the formation of turbulence. This volume is determined by the cutoff angle (*θ _{c}*) since the resonant surface furthest from the center of the current sheet that will be unstable is given by $xc=\lambda arctanh(\mu c)$, where $\mu c=tan\u2009\theta cBoy/B\xafoz$. We have shown that the differential equation Eq. (22) provides a good description of the cutoff in comparison to the complete linear VMID solution. However, this requires solving the differential equation numerically. Here, we show that a simple analytic formula for

*θ*can be obtained based on the results described previously in this section, which have established that damping occurs when the inner tearing layer becomes as wide as the current sheet (

_{c}*λ*).

The width of the inner layer ($\lambda \u2032$) can be estimated using Eq. (54). We approximate the stability condition as $|T(x=\lambda \u2032)|=c$, where *c* is a constant that will be determined from a fit to the numerical solutions. Since $T\u221d\lambda /|x\u2212xs|$, the numerical value of the cutoff width can be absorbed into the fit parameter. The outer layer of this region is determined by the condition that the argument of *Z*_{1} becomes small, in which case $Z1\u21921$; see Fig. 9. Applying the approximation that $\omega \u2243\omega R\u2243kyUe$ at oblique angles and that $|x\u2212xs|/\lambda \u22431\u2212arctanh(\mu )\u22431\u2212\mu $, the cutoff condition is then

Finally, a simplification can be made based on the previous observation that the background density does not significantly influence the cutoff angle (see Fig. 10), so the term in parentheses can be approximated as 1. Applying this, Eq. (55) provides the following estimate for the cutoff angle

such that instability is expected when $\theta \u2272\theta c=arctan(\mu cB\xafoz/Boy)$. Alternatively, the critical temperature ratio at a given angle and layer thickness ($\lambda /\rho i$) can be related to a critical electron drift speed since $Ue=\u2212TeTiUi=\u2212TeTi\rho i\lambda $.

Figure 12 shows a comparison between the cutoff angle predicted by Eq. (56) and the numerical results of Eq. (22). The fit parameter used was *c* = 0.4, which is a reasonable number based on the fact that *T*(*x*) = *c* defines where the effective width of *T* is measured (see Fig. 9). Also shown in the figures is the cutoff from the DL dispersion relation, which is simply *μ* = 1, leading to

Figure 12(a) shows that Eq. (56) captures both the large $Boy/B\xafoz$ scaling and the value at which the deviation from this scaling occurs (for $Boy/B\xafoz<1$). Importantly, Fig. 12(b) shows that Eq. (56) accurately captures the temperature ratio dependence of the cutoff angle (when $Ti/Te$ is sufficiently large). Equation (55) shows that the thickness of the inner region scales as $1/(Ti/Te)2$ when $Ti/Te\u226b1$. That is, the inner region becomes very thin at a large temperature ratio, and one should expect in this limit that the boundary layer scale separation that the analytic dispersion relations (such as DL) are based on will be valid. Indeed, this is what is observed, as the two predictions merge at a high temperature ratio; see Fig. 12(b). For $Ti/Te\u226a1$, a dependence on the layer thickness is observed that is not captured by Eq. (56). This figure provides further evidence that broadening of the tearing layer, associated with the diamagnetic drift, is responsible for the stabilization at oblique angles. It also shows that Eq. (56) provides a reliable estimate for the cutoff angle of oblique collisionless tearing modes.

## VI. CONCLUSIONS

It was shown that oblique collisionless tearing instabilities of a Harris current sheet are stabilized due to the density-gradient-driven diamagnetic drift. The complicated Vlasov-Maxwell system of equations was simplified to a single second order differential equation for the vector potential perturbation, which captures the essential features of the dispersion relation obtained from the full numerical solution of the linear Vlasov-Maxwell equations. Theories from previous literature are not able to capture this stabilization because they are based on a boundary layer scale separation between an inner tearing layer and an outer ideal MHD region. The inner layer was observed to broaden at oblique angles, invalidating this approximation.

These results are applicable to studies of magnetic reconnection in space and fusion plasmas, where tearing and plasmoid instabilities are known to initiate and accelerate the rate of magnetic reconnection. They provide a method to compute the spectrum of oblique modes, which contributes to theories of turbulence generation and particle acceleration. They are particularly relevant to understanding the magnetopause, where reconnection occurs within thin asymmetric layers in which diamagnetic drifts are strong. In this context, Galeev coined the term “percolation” to denote the regime in which many overlapping islands formed by tearing modes give rise to turbulent particle transport.^{25} Our results show that diamagnetic drifts broaden the inner tearing layer, so that in current sheets with strong diamagnetic drifts, only resonant surfaces with small $k\xb7U$ are expected to be unstable. This leads to a narrower spectrum of instabilities, i.e., a smaller volume of space susceptible to tearing, than was previously expected. A simple analytic estimate for the stabilization threshold was provided in Eq. (56), which can be used to predict the volume of space susceptible to linear tearing instabilities. The results also highlight the importance of the equilibrium current sheet properties in understanding tearing instabilities, and by extension, the overall reconnection process. Specifically, drift-stabilization is not present in other common equilibria used to model such processes, such as the force-free current sheet.^{39}

Interesting questions remain to be answered, such as to what extent fluid theory can describe collisionless tearing and to what extent purely kinetic effects such as Landau damping play a role. For instance, Landau damping is a non-negligible, and an even prominent, contributor to the dispersion relation [arising in the plasma dispersion function of Eq. (22)] in the kinetic theory described in this work. Yet, two fluid theories have found success in modeling collisionless tearing modes in some geometries.^{53} In one case, the predicted growth rates differ only by a factor of $\pi $.^{61} It would be useful to understand if two-fluid theory captures the stabilization of oblique modes observed in this work. If so, can it provide a simplified dispersion relation? It would also be important to revisit the linear tearing stability theory in the presence of temperature gradients. Although the Harris sheet does not have them, both the magnetopause and fusion machines have density and temperature gradients.

## ACKNOWLEDGMENTS

The authors thank Dr. V. Roytershteyn for helpful discussions on this work. This material is based upon work supported by the U.S. Department of Energy, Office of Science, Office of Fusion Energy Sciences under a Postdoctoral Research Program administered by the Oak Ridge Institute for Science and Education, DOE Award No. DE-SC0016159, and by a NSF Grant No. AGS-0962698. Contributions from W.D. were supported by the Basic Plasma Science Program from the DOE Office of Fusion Energy Sciences.

## References

Reference 25 takes $b\u0302s$ to be the $\u2225$ direction in defining $\delta A\u0302\u2225$. However, the $k\u2225$ that is cited must be parallel to $Bo$ since $k\xb7b\u0302s=0$.