Collisionless kinetic theory of oblique tearing instabilities

Tearing instabilities are magnetic reconnection events that occur in transition layers where the frozen flux condition of ideal magnetohydrodynamics (MHD) is violated.¹ In tokamaks, they contribute to the basic magnetic field topology,² transport,^3–6 rotation,⁷ and may trigger disruptive events.⁸ In space, they are responsible for initiating large-scale magnetic reconnection in Earth's magnetosphere,⁹ as well as flares in the solar corona.¹⁰ If the plasma is sufficiently dense or cool, tearing instabilities are well described by resistive,^1,13–17 or visco-resistive¹⁸ 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.³³ Recent particle-in-cell (PIC) simulations, together with numerical solutions of the linearized Vlasov equation,³⁴ have revealed the surprising result that the spectrum of oblique modes is much more narrow than standard analytic theories predict.^24,25 Here, $θ=arctan(ky/kz)$ is the angle of obliquity and k the wavevector; see Fig. 1. The resonant surface locations $xs=λarctanh(μ)$⁠, where $μ≡ tan θBoy/B¯oz$⁠, 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,³⁹ 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)²⁴ and Galeev et al.,²⁵ 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⁵⁰ and finite-element expansions^9,51 of the eigenfunction. In this paper, the results from the finite-element code⁵¹ 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.⁵² Simulations and collisionless kinetic theory, similar to what are discussed here, were found to agree well for a force-free equilibrium.³⁹ Similarly, two-fluid theory was found to accurately capture the dispersion of oblique modes in this case.⁵³ 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.⁵⁴ 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.²⁹ 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.

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=∂/∂t+v·∇+(qs/ms)(Eo+v×Bo/c)·∇v$ 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)ŷ+Bzo(x)ẑ$ and that the lowest-order distribution functions are Maxwellians

where the flow velocity is $Us=Us,y(x)ŷ+Us,z(x)ẑ$⁠. We write $δE$ and $δ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⁵⁰ 

where

$τ=t−t′, δϕ̃=δϕ(k,x′,t)$ and $δÃ=δA(k,x′,t)$⁠. The single particle characteristics (⁠$x−x′$⁠) 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 (⁠$δρ=∑sqsδns$⁠, where $δns=∫d3v δfs$⁠) and current density (⁠$δJ=∑sqsδΓs$⁠, where $δΓs=∫d3v vδfs$⁠), which are then used in Maxwell's equations in the Lorenz gauge

As a next level of approximation, gyroradius effects are neglected (⁠$ρs→0$⁠). In this case, the single particle characteristics are $x−x′=v∥τb̂+O(ρs)$⁠, where $b̂=Bo/|Bo|$⁠, and Eq. (5) reduces to

assuming $ℑ{ω}>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 $∇2δϕ=−4πδρ$ and $∇2δA=−(4π/c)δJ$⁠.

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

in which $λDs2=∑s(4πqs2nso)/Ts$ is the square of the Debye length for species s, and $λD−2=∑sλDs−2$ gives the total Debye length

The coordinate rotation $(x,y,z)→(x,b,η)$ has been applied, where $b̂=byŷ+bzẑ$ is parallel to the equilibrium magnetic field and $η̂=x̂×b̂=−bzŷ+byẑ$⁠. The plasma dispersion function⁵⁵ 

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

where $ωps=4πqs2nso/ms$ is the plasma frequency of species s. Also, note that $δBx=i(k∥δAη−kηδA∥)$ and $δE∥=−ik∥δϕ+i(ω/c)δA∥$⁠.

Explicit evolution equations for $δAy$ and $δAz$ can be obtained by projecting Eq. (12) along $ŷ$ and $ẑ$⁠. 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̂$⁠, which is simply $k̂:b̂s×k̂=x̂$⁠. Thus, $(x̂,b̂s,k̂)$ provides a convenient orthogonal coordinate system. Applying this to Eq. (12) gives

and

in which

$Vbs=b̂s·V$ and $Vk=k̂·V$⁠. Here, $δAk=k·δA$

and

The set of three coupled Eqs. (9), (14), and (15) provides a closed description for the small gyroradius limit.

Although $δE∥$ in Eq. (14) depends on $δϕ, δÂ$ and $δAk$⁠, consideration of two-scale features of a tearing layer can be used to justify an approximate decoupling of $δÂ$ from $δϕ$ and $δAk$⁠. In the outer region, far from a resonant surface, magnetic flux is frozen into the plasma, so the parallel electric field vanishes $δE∥→0$⁠. In this ideal MHD-like outer region, the right side of Eq. (14) is negligible and $δÂ$ is the only remaining dependent variable. In the inner region near a resonant surface, $k∥≪k$⁠, 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 $ω2/(kvTs)2≪1$⁠. The right side contributes only in the inner region where $k∥$ is small. Thus, Eq. (22) is an approximate equation spanning the layer that is decoupled from $δAk$ and $δϕ$⁠. The validity of this scale separation will be evaluated in Sec. V.

We utilize the following properties of the Harris equilibrium³³ 

where an additional uniform background density (n_b) for each species is included, $n(x)=n¯(x)+nb$⁠. Noting that $Us,η=−bzUs$⁠, 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·Bo$⁠. The individual components of interest are then

and $Vk=tan θVbs$⁠.

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=±∞→xs$ and determining the dispersion relation from the requirement that $δÂ′(xs+)=−δÂ′(xs−)$⁠.

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 $θ=0°$⁠, 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 $θ=20°$⁠, 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 (⁠$λ=ρi, 2ρi$ and $5ρ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 $θ=0°$ 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 $θ≃30°$⁠) and that the values of the growth rate in units of Ω_ci take larger values at a lower mass ratio. Here, $Ωci=eB¯oz/(mic)$ is the ion gyrofrequency based on the asymptotic magnetic field strength $B¯oz$⁠. 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−1$⁠), 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.⁵⁶ Another apparent difference is that for larger current sheet widths (⁠$λ=5ρi$ and, to a lesser extent, $λ=2ρ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 $θ=0°$ and generally predicts a flat profile. One distinction observed at $λ=5ρ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.⁵⁴ 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.²⁴ 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.

Figure 3(a) shows profiles of the real and imaginary components of the vector potential $δÂ(x)$ across the current sheet for modes centered at resonant surfaces corresponding to $θ=0°, 5°, 10°$⁠, and $15°$⁠. 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 (⁠$θ=0°$⁠) 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)¹ 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 (⁠$θ>0$⁠) have a finite $ℑ{δÂ(x)}$⁠, and that their profile is also asymmetric with respect to the resonant surface. The values of this component are negative (the plot shows $ℑ{−δÂ}$⁠). Figure 3(b) shows that if the complex vector potential is represented as $δÂ(x)=|δÂ(x)| exp [iϕ(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 $δÂ$ 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.

The seminal kinetic tearing mode theory of Drake and Lee (DL)²⁴ 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→ζs$⁠, which implies $ω/k∥≫Us,∥$⁠. Note that $Z1(ζs)=−Z′(ζs)/2=1+ζsZ(ζs)$⁠.

The DL dispersion relation follows from applying the constant-ψ approximation. In this approximation, $δÂ$ on the right side of Eq. (30) is assumed constant across the inner layer. Dividing by this and integrating across the layer give $∫dx δÂ″/δÂ(xs)≃Δ′$⁠, and

The full solution of $k∥$ 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∥$ and ω_ps are expanded about the resonant surface location: $k∥≃k(x−xs)/ls$⁠, and $ωps=ω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≡c/ω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 $ω=ωr+iγ$⁠, where

is the real frequency of the tearing mode, and

is the growth rate. Here, $n̂b≡nb/n¯o$ and $α≡Teme/(Timi)$⁠. For the plots in Fig. 2, which show $ω/Ωci$⁠, it is convenient to note that $kyUe=Ωcikλ sin θ(Te/Ti)(ρi/λ)2$ and $vTede2|xs=Ωciρi3α(1+Te/Ti)/(1−μ2+n̂b)$⁠.

The tearing stability index $Δ′$ must be supplied from a model for the outer region of the boundary layer. This is typically assumed to obey ideal MHD (FKR)¹ 

Recall that $Vbs=F″/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 $δÂ(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 $Δ′$ 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 $Δ′$ 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 (⁠$θ=0°$⁠) 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 ζ_s in DL theory. We have also numerically solved a version of Eq. (22) with $Z(ws)$ replaced by $Z(ζ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.

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∥≪k$ and (2) $Us/vTs≪1$⁠. 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∥≪k$ assumption implies $δA≃δÂb̂s$ and $b̂s≃b̂≃ẑ$⁠. Applying these along with assumption (2), the $b̂s$ 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.⁵⁹ 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→0)$⁠, i.e., the parallel mode. The cause of this disagreement is that the $k∥≪k$ 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δE∥≃−ωδÂ/c$ approximation from Eq. (19), which was validated in Sec. III, to Eq. (44). This gives

which is the same as the $Us,∥≪ω/k∥vTs$ 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 $Δ′$ obtained from the numerical solutions of Eqs. (48) and (39), and the model solutions from Eq. (40) and Galeev's approximate solution⁶⁰ 

where

This figure emphasizes that $Δ′$ obtained from Eq. (48) differs from the ideal MHD solution. It also shows that the approximate solution from Eq. (49) requires that $kλ$ be sufficiently large to accurately model the numerical solution.

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.

Figure 7 shows a comparison of the dispersion relation from Lee, Mahajan and Hazeline⁴⁵ 

where $k∥′=k/ls, ρs=cs/Ωci$⁠, and $ωo=ℜ{ωDL}$⁠. Here, $ω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∥δϕ≪ωckηkδÂ$⁠. Although wide tearing layers at oblique angles access a range of values far from the current sheet, so $k∥/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 $ωpe/ωce$⁠. Thus, one can alter the magnitude of the electrostatic terms with this parameter. Nevertheless, the computed growth rate remained independent of this parameter.⁹ 

The DL theory applies the constant-ψ approximation to obtain an analytic dispersion relation. This works best if the inner layer is very thin because $δÂ$ 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 $δÂ$ near the resonant surface. These theories are called large $Δ′$⁠, or nonconstant-ψ approximations. It was observed in the MHD theory that a nonconstant-ψ approximation was required to accurately model oblique modes.⁵² 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 $Δ′$ theory of kinetic collisionless tearing modes has been computed by Mahajan et al.,⁴⁹ which provides the dispersion relation

Here, $γDL=ℑ{ωDL}, ωo=ℜ{ωDL}$⁠, and $vA=B¯oz/4πn¯omi$ 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 $Δ′$ extension of DL theory does not appear to be able to explain the observed stabilization.

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 $λ=ρi$ to $5ρ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 $ρe/λ∼me/mi≪1$⁠. 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 $δi/ρ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 $δs≃|ω|ls/(kvTs)$⁠, where ω was computed from Eqs. (37) and (38).²⁴ 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 (⁠$θ=0°$⁠) 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

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 $δi/λ≪1$ 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 ζ_s in the 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∥≃k(x−xs)/ls$ and $ωps≃ω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 $λ=2ρ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.