Collisionless Kinetic Theory of Oblique Tearing Instabilities
Scott D. Baalrud Amitava Bhattacharjee William Daughton
physics.plasm-ph



























Department of Physics and Astronomy, University of Iowa, Iowa City, Iowa 52242, USA






Department of Astrophysical Sciences, Princeton University, Princeton, New Jersey 08544, USA



Los Alamos National Laboratory, Los Alamos, New Mexico 87545, USA
















The linear dispersion relation for collisionless kinetic tearing instabilities is calculated for a 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 (θ≠ 0). 
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. 







52.35.Vd,52.55.Tn,94.30.cp










Collisionless Kinetic Theory of Oblique Tearing Instabilities
    W. Daughton
    December 30, 2023
=============================================================









§ INTRODUCTION




Tearing instabilities are magnetic reconnection events that occur in transition layers where the frozen flux condition of ideal magnetohydrodynamics (MHD) is violated. <cit.> 
In tokamaks, they contribute to the basic magnetic field topology, <cit.> transport,<cit.> rotation, <cit.> and may trigger disruptive events. <cit.> 
In space, they are responsible for initiating large-scale magnetic reconnection in Earth's magnetosphere, <cit.> as well as flares in the solar corona. <cit.> 
If the plasma is sufficiently dense or cool, tearing instabilities are well described by resistive, <cit.> or visco-resistive <cit.> MHD.
At lower densities or higher temperatures, two-fluid or Hall MHD effects modify the resistive reconnection rate. <cit.> 
At even lower densities or higher temperatures, resistivity plays no role and reconnection is determined entirely by collisionless two-fluid or kinetic effects. <cit.>  
Space plasmas are frequently in this latter collisionless regime, and fusion plasmas can also reach this regime. <cit.> 

In this paper, we revisit the linear kinetic dispersion relation for collisionless tearing instability of a Harris current sheet. <cit.> 
Recent particle-in-cell (PIC) simulations, together with numerical solutions of the linearized Vlasov equation, <cit.> have revealed the surprising result that the spectrum of oblique modes is much more narrow than standard analytic theories predict.<cit.> 
Here, θ = arctan (k_y/k_z) is the angle of obliquity and k the wavevector; see Fig. <ref>. 
The resonant surface locations x_s = λarctanh(μ), where μ≡tanθ B_oy/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. <cit.> 
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). <cit.> 
Resonant surfaces unstable to linear tearing modes also become the sites of electron layers in the nonlinear regime,<cit.> 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,<cit.> but none of them captures this effect. We focus on comparing with Drake and Lee,<cit.> and Galeev et al, <cit.> 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 <cit.> and finite-element expansions <cit.> of the eigenfunction. In this paper, the results from the finite-element code <cit.> 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 predicated 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. <cit.> 
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. <cit.> Simulations and collisionless kinetic theory, similar to what are discussed here, were found to agree well for a force-free equilibrium. <cit.> 
Similarly, two-fluid theory was found to accurately capture the dispersion of oblique modes in this case. <cit.> 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<cit.> 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 <cit.>. 
There is no temperature gradient in the Harris equilibrium.

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




§ BASIC EQUATIONS 





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

    d δ f_s/dt = - q_s/m_s( δE + v×δB/c) ·∂ f_so/∂v

where d/dt = ∂/∂ t + v·∇ + (q_s/m_s)(E_o + v×B_o/c) ·∇_v is the convective derivative. We will later focus on the Harris equilibrium, but in this section require only that the equilibrium magnetic field can be written in the form B_o = B_yo (x) ŷ + B_zo (x) ẑ and that the lowest-order distribution functions are Maxwellians 

    f_so = n_s (x)/π^3/2 v_Ts^3exp[- (v - U_s)^2/v_Ts^2]

where the flow velocity is U_s = U_s,y (x) ŷ + U_s,z (x) ẑ. We write δE and δB in terms of scalar and vector potentials 

    δE = - ∇δϕ - 1/c∂δA/∂ t   and   δB = ∇×δA

and apply a normal mode analysis of the form

    δϕ = δϕ̅ (x) exp (- i ω t + i k_y y + i k_z z) , 
    δA = δA̅(x) exp (- i ω t + i k_y y + i k_z z)  .


Inserting Eq. (<ref>) into (<ref>) and integrating along the characteristics, which are the unperturbed particle orbits, provides<cit.>

    δ f_s = - q_s f_so/T_s[ δϕ - U_s ·δA/c + i (ω - k·U_s) S ]

where

    S = ∫_0^∞ dτ ( δϕ̅ - v^'·δA̅/c ) e^i ωτ - i k· (x - x^') ,

τ = 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 is computed numerically from Eq. (<ref>) for the exact single particle characteristics using the orbit integration technique described in Refs. daug:99 and daug:03. 
Equations (<ref>) and (<ref>) are used to calculate the perturbed charge density (δρ = ∑_s q_s δ n_s where δ n_s =  ∫ d^3v δ f_s) and current density (δJ = ∑_s q_s δΓ_s where δΓ_s = ∫ d^3v vδ f_s), which are then used in Maxwell's equations in the Lorenz gauge


    ∇^2 δϕ - 1/c^2∂^2 δϕ/∂ t^2   =    - 4πδρ , 
    ∇^2 δA - 1/c^2∂^2 δA/∂ t^2   =    - 4π/cδJ .


The resulting set of three coupled integro-differential equations for δϕ and δA are then solved numerically as an eigenvalue problem for the eigenvectors and dispersion relation using a finite element approach.<cit.> The boundary conditions on δϕ and δA are that they asymptote to a constant value (zero) as x = ±∞. 




 §.§ Small gyroradius limit 


As a next level of approximation, gyroradius effects are neglected (ρ_s → 0). In this case, the single particle characteristics are x - x^' = v_∥τb̂ + 𝒪(ρ_s), where b̂ = B_o/|B_o|, and Eq. (<ref>) reduces to 

    S = i (δϕ - v_∥δ A_∥/c)/ω - k_∥ v_∥ + 𝒪(ρ_s)

assuming {ω} >0. To lowest-order in this gyroradius expansion, Eq. (<ref>) can then be written

    δ f_s = - q_s f_so/T_s[ δϕ - U_s ·δA/c - ( δϕ - v_∥δ A_∥/c) ω - k·U_s/ω - k_∥ v_∥]  .

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 

    [∂_x^2    - (k^2 + λ_D^-2)] δϕ = ∑_s λ_Ds^-2{ - U_s ·δA/c
       + ω - k·U_s/k_∥ v_Ts[ δϕ Z(w_s) - v_Ts/cδ A_∥ (1 + ζ_s Z(w_s)) ] }  ,

in which λ_Ds^2 = ∑_s (4π q_s^2 n_so)/T_s is the square of the Debye length for species s, λ_D^-2 = ∑_s λ_Ds^-2 gives the total Debye length,

    w_s≡ω - k_∥ U_s,∥/k_∥ v_Ts   and   ζ_s ≡ω/k_∥ v_Ts  .

The coordinate rotation (x,y,z) → (x,b,η) has been applied where b̂ = b_y ŷ + b_z ẑ is parallel to the equilibrium magnetic field and η̂ = x̂×b̂ = -b_z ŷ + b_y ẑ. The plasma dispersion function<cit.> 

    Z_n (w) = 1/√(π)∫_-∞^∞ dt t^n e^-t^2/t - w

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

    ∇^2    δA = ∑_s 2 ω_ps^2/c k_∥ v_Ts^2{ i δ B_x U_s,η/cU_s 
       + i δ E_∥[ U_s + (ω - k·U_s)/k_∥( U_s/v_Ts Z(w_s) + Z_1 (w_s) b̂) ] }
 
where ω_ps = √(4π q_s^2 n_so/m_s) is the plasma frequency of species s. Note also that δ B_x = i (k_∥δ A_η - k_ηδ A_∥) and δ E_∥ = -i k_∥δϕ + i (ω/c) δ A_∥. 

Explicit evolution equations for δ A_y and δ A_z can be obtained by projecting Eq. (<ref>) along ŷ and ẑ. However, more compact expressions can be obtained by projecting Eq. (<ref>) in the direction along B_o at a resonant surface 

    b̂^s = B_o (x_s) / B_o (x_s) = cosθŷ - sinθẑ

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. (<ref>) gives 

    [∂_x^2     - (k^2 + V_b^s)] δÂ = i δ E_∥∑_s 2 ω_ps^2/c k_∥ v_Ts[ b̂^s ·U_s/v_Ts
       + (ω - k·U_s)/k_∥ v_Ts( b̂^s ·U_s/v_Ts Z (w_s) + k_η/k Z_1 (w_s) ) ]

and 

    (∂_x^2     - k^2 ) δ A_k = δÂ V_k + i δ E_∥∑_s 2 ω_ps^2/c k_∥ v_Ts[ k̂·U_s/v_Ts
       + (ω - k·U_s)/k_∥ v_Ts( k̂·U_s/v_Ts Z (w_s) + k_∥/k Z_1 (w_s) ) ]

in which 

    V = ∑_s 2 ω_ps^2/c^2k/k_∥U_s,ηU_s/v_Ts^2  ,

V_b^s = b̂^s ·V and V_k = k̂·V.
Here, δ A_k = k·δA,

    δÂ = b^s ·δA =  (k_ηδ A_∥ - k_∥δ A_η)/k = i δ B_x/k ,

and

    i δ E_∥ =  k_∥δϕ - ω/c( k_η/kδÂ + k_∥/kδ A_k )   .


The set of three coupled equations (<ref>), (<ref>) and (<ref>) provide a closed description for the small gyroradius limit. 




 §.§ Decoupled magnetic field approximation


Although δ E_∥ in Eq. (<ref>) depends on δϕ, δÂ and δ A_k, consideration of two-scale features of a tearing layer can be used to justify an approximate decoupling of δÂ from δϕ and δ A_k. 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. (<ref>) is negligible and δÂ is the only remaining dependent variable. In the inner region near a resonant surface, k_∥≪ k, which implies

    i δ E_∥ = k_∥δϕ - ω/ck_η/kδÂ - ω/ck_∥/kδ A_k ≃ - ω/cδÂ ,
 
and

    b̂^s ·U_s = b_∥^s U_s,∥ + b_η^s U_s,η = k_η/k U_s,∥ - k_∥/k U_s,η≃ U_s,∥

so 

    k/k_ηb̂^s ·U_s/v_Ts Z(w_s) + Z_1(w_s) ≃ 1 + ζ_s Z(w_s) .

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

    [∂_x^2    - (k^2 + V_b^s)] δÂ =   
       - δÂ∑_s 2 ω_ps^2/c^2(ω - k·U_s)/k_∥ v_Tsζ_s [1 + ζ_s Z(w_s)]  .

to describe the inner region. 

Although Eq. (<ref>) 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/(kv_Ts)^2 ≪ 1. The right side contributes only in the inner region where k_∥ is small. Thus, Eq. (<ref>) is an approximate equation spanning the layer that is decoupled from δ A_k and δϕ. 
The validity of this scale separation will be evaluated in Sec. <ref>. 




 §.§ Harris Equilibrium 


We utilize the following properties of the Harris equilibrium<cit.>


    B_o    =    B_oyŷ + B̅_oztanh (x/λ) ẑ , 
    n̅   =   n̅_o ^2 (x/λ) , 
    B̅_oz^2    =    8πn̅_o (T_e + T_i) , 
    U_s    =    - 2cT_s/(q_s B̅_ozλ) ŷ


where an additional uniform background density (n_b) for each species is included, n(x)=n̅(x) + n_b. Noting that U_s,η = -b_z U_s, Eq. (<ref>) reduces to

    V = - 4πn̅ e^2/c^2(U_i/T_i)^2 (T_e + T_i) k B_z/k_∥ B_oŷ ,

and Ampére's law gives

    B_z^'' = - 4π e/cU_i/T_i (T_e + T_i) n̅^' ,

where “prime” denotes a derivative. 
Taking a derivative of Eq. (<ref>) and utilizing Eqs. (<ref>) and (<ref>) shows

    B_z^'' = - 4πn̅ e^2/c^2( U_i/T_i)^2 (T_e + T_i) B_z  .

Putting this into Eq. (<ref>) provides

    V = F^''/Fk/k_zŷ

where F = k·B_o. The individual components of interest are then

    V_b^s = F^''/F = - 2/λ^2^2 (x/λ ) tanh (x/λ )/tanh (x/λ) - tanh (x_s/λ) ,

and V_k = tanθ V_b^s.





§ NUMERICAL SOLUTIONS 
 

In this section, the dispersion relation computed from the VMID model of Sec. <ref> is compared with that computed from the approximate layer equation of Sec. <ref> [Eq. (<ref>)]. 
The VMID computations were calculated using the method and code described in Refs. daug:99 and daug:03. 
Solutions of Eq. (<ref>) were obtained by numerically integrating from x=±∞→ x_s, and determining the dispersion relation from the requirement that δÂ^'(x_s^+) = -δÂ^'(x_s^-).




 §.§ Growth rate profiles
 



The VMID solutions shown in Fig. <ref>c demonstrate that at realistic mass ratio (m_i/m_e = 1836), the growth rate is a maximum at θ = 0^∘, corresponding to the usual parallel mode at the center of the current sheet, and monotonically decreases with increasing angle of obliquity. The growth rate falls off slowly with 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 linear over the entire range of unstable angles at each sheet width. 

Figure <ref>b shows results of analogous simulations, but at lower mass ratio (m_i/m_e = 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 lower mass ratio. Here, Ω_ci = eB̅_oz/(m_ic) 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 larger mass ratio. 

Figure <ref>a shows that at mass ratio of unity m_i/m_e=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. <cit.> Another apparent difference is that for the 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. <ref> also show the predictions of the approximate layer equation from Eq. (<ref>). Comparison with the Vlasov simulation curves shows that, despite the severity of some the approximations made in Secs. <ref> and <ref>, 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 (m_i/m_e =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. (<ref>) accurately predicts the growth rate of the central mode at θ = 0^∘, and generally predicts a flat profile. One distinction observed at λ=5ρ_i and m_i/m_e = 1836 is that Eq. (<ref>) predicts the most unstable mode is oblique (near the cutoff), whereas the Vlasov solutions predict a monotonically decreasing profile. Finally, much poorer agreement between Eq. (<ref>) 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 U_s. This is confirmed by the dashed lines in Fig. <ref>, which show solutions of Eq. (<ref>) with U_s=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. (<ref>) 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.<cit.>
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 <cit.>
The results of Eq. (<ref>) in the absence of diamagnetic drift agree well with the predictions of the standard kinetic tearing mode theory of Drake and Lee <cit.>. 
In Sec. <ref> 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. (<ref>) accurately captures the essential feature of the unstable region of the current layer, the remainder of this work will focus on this equation.





 §.§ Vector potential and phase shifts




Figure <ref>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 m_i/m_e=1836 and the other parameters chosen were the same as in Fig. <ref>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)<cit.> 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 figure 7 of Ref. baal:12), and is associated with 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 <ref> 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 <ref>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. <ref>. 





§ RELATION TO PREVIOUS THEORIES
 




 §.§ Drake and Lee 


The seminal kinetic tearing mode theory of Drake and Lee (DL)<cit.> 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 

    Δ^'≡ [δÂ^' (x_s^+) - δÂ^' (x_s^-)]/δÂ (x_s)

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

    δÂ^'' = - δÂ∑_s 2ω_ps^2/c^2(ω - k·U_s)/k_∥ v_Tsζ_s Z_1(ζ_s) .

Equation (<ref>) agrees with the inner region of our decoupled layer equation (<ref>) in the limit that w_s →ζ_s, which implies ω/k_∥≫ U_s,∥. Note, Z_1(ζ_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. (<ref>) is assumed constant across the inner layer. Dividing by this and integrating across the layer gives ∫ dx δÂ^''/δÂ (x_s) ≃Δ^', and 

    Δ^'≃∑_s (ω - k·U_s)/ω∫_-∞^∞ dx ω_ps^2/c^2ζ_s^2 Z^' (ζ_s) .

The full solution of k_∥ for the Harris equilibrium described in Sec. <ref> is 

    k_∥/k = cosθμ + tanh (x/λ)/√( B_oy^2/B̅_oz^2 + tanh^2 (x/λ)) .

To integrate the right side of Eq. (<ref>) an additional assumption is applied, where the space-dependent functions k_∥ and ω_ps are expanded about the resonant surface location: k_∥≃ k (x-x_s) /l_s, and ω_ps = ω_ps(x_s), where

    l_s =  k/k_∥^' (x_s) = λ B_oy/ B̅_oz/cos^2θ (1 - μ^2)  .

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

In Ref. drak:77, the integral in Eq. (<ref>) is evaluated by applying the above approximations.
The result is 

    ∫_-∞^∞ dx ζ_s^2 Z^' (ζ_s) = - 2i √(π)ω l_s/kv_Ts .

With this, Eq. (<ref>) provides the dispersion relation

    Δ^' = -2i√(π)∑_s l_s (ω -k·U_s)/k v_Ts d_s^2

where d_s≡ c/ω_ps is the skin depth associated with species s. 
Rearranging, the real frequency and growth rate are 

    ω_DL = ∑_s k·U_s/(v_Ts d_s^2)/∑_s 1/ (v_Ts d_s^2) + i k Δ^'/2 √(π) l_ s ∑_s 1/(v_Ts d_s^2)  .




Applying the expressions for a Harris sheet with a background density from Sec. <ref> provides ω = ω_r + i γ, where

    ω_r = k_y U_e (1 - α T_i/T_e)/(1 + α) [1 + n̂_b/(1-μ^2)]

is the real frequency of the tearing mode, and 

    γ = k Δ^' v_Te d_e^2|_x_s/2 √(π) l_s (1+α)

is the growth rate. Here, n̂_b ≡ n_b/n̅_o and α≡√(T_em_e/(T_im_i)). For the plots in Fig. <ref>, which show ω/Ω_ci, it is convenient to note that k_y U_e = Ω_ci kλsinθ (T_e/T_i) (ρ_i/λ)^2 and v_Te d_e^2|_x_s = Ω_ciρ_i^3 α (1+T_e/T_i)/(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) <cit.>

    δÂ^'' - (k^2 + V_b^s) δÂ = 0 .

Recall that V_b^s = F^''/F is provided in Eq. (<ref>). 
Figure <ref> shows the solution of Eq. (<ref>) 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<x_s) than on the high magnetic field side (x>x_s). A model for Δ^' based on Eq. (<ref>), which considers oblique modes, was developed in Ref. baal:12 

    Δ^'_H ≃2/λ( 1 + μ^2/kλ - k λ) .

This was validated against numerical solutions of Eq. (<ref>) over a broad range of conditions (see Fig. 3 of Ref. baal:12), and a similar validation is shown in Fig. <ref>b. It was used to evaluate Δ^' for the Drake and Lee theory curves shown in Fig. <ref>. 

Figure <ref> 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. (<ref>) 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. (<ref>) is w_s, whereas it is ζ_s in DL theory. We have also solved numerically a version of Eq. (<ref>) with Z(w_s) 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.  




 §.§ Galeev et al




Previous researchers<cit.> 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. (<ref>) by applying the following assumptions: (1) k_∥≪ k and (2) U_s/v_Ts≪1. These imply 

    U_s ·δA/c≪(ω - k·U_s)/k_∥ cδ A_∥   and   w_s ≃ζ_s ,

so Eq. (<ref>) reduces to

    (∂_x^2    - k^2 - λ_D^-2) δϕ = 
       ∑_s λ_Ds^-2(ω - k·U_s)/k_∥ v_Ts[ δϕ Z(ζ_s) - v_Ts/c Z_1 (ζ_s) δ A_∥] .

For Ampére's equation, it is convenient to first rearrange Eq. (<ref>) to the form

    ∇^2 δA   = ∑_s 2 ω_ps^2/c v_Ts^2{( δϕ - U_s ·δA/c) U_s  
       + (ω - k·U_s)/k_∥[ δϕ( U_s/v_Ts Z(w_s) + Z_1 (w_s) b̂) 
       - δ A_∥v_Ts/c( U_s/v_Ts (1 + ζ_s Z(w_s)) + ζ_s Z_1 (w_s) b̂) ] } .

The k_∥≪ k assumption implies δA≃δÂb̂^s and b̂^s ≃b̂≃ẑ. Applying these along with assumption (2), the b̂^s projection of Eq. (<ref>) reduces to 

    (∂_x^2 - k^2 - V_ad ) δÂ = i δ E_∥∑_s λ_Ds^-2(ω - k·U_s)/k_∥^2 c Z_1 (ζ_s)

in which 

    V_ad = - ∑_s 2 ω_ps^2/c^2(b̂^s ·U_s)^2/v_Ts^2

is what Galeev calls the adiabatic interaction term. 

Equations (<ref>) and (<ref>) are the two basic equations that this line of previous collisionless oblique tearing mode research<cit.> worked from.<cit.> This set of two coupled equations is substantially different than any of the approximate equations derived in Sec. <ref>. The most significant difference is in the description of the outer region. For an ideal MHD outer region, the V_ad is replaced by V_b^s, which reduces to Eq. (<ref>) for a Harris equilibrium. However, the adiabatic response used in Eq. (<ref>) is 

    V_ad = - 2/λ^2cos^2 θ ^2 (x/λ)

for a Harris equilibrium. 
This agrees with the ideal MHD response, Eq. (<ref>), only in the limit of a standard symmetric layer (x_s → 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. (<ref>), which was validated in the last section, to Eq. (<ref>). This gives

    (∂_x^2 - k^2 - V_ad) δÂ = - δÂ∑_s 2 ω_ps^2/c^2(ω - k·U_s)/k_∥ v_Tsζ_s Z_1 (ζ_s)

which is the same as the U_s,∥≪ω/k_∥ v_Ts limit of Eq. (<ref>) except in the outer region (left side) where V_ad replaces V_bs. 


Figure <ref>a show solutions of the outer region of the Galeev equation 

    δÂ^''  - (k^2 + V_ad) δÂ = 0
 
in comparison to ideal MHD from Eq. (<ref>). 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. <ref>b shows a comparison between Δ^' obtained from the numerical solutions of Eq. (<ref>), (<ref>), and the model solutions from Eq. (<ref>) and Galeev's approximate solution <cit.>

    Δ^'_G ≃ - 2 (kλ + ν)/λΓ [(k λ + ν)/2] Γ [(1 + kλ - ν)/2]/Γ [ (kλ - ν)/2] Γ [(1 + kλ + ν)/2]

where 

    ν≡1/2( √(1 + 8 cos^2 θ) - 1 )   .

This figure emphasizes that the Δ^' obtained from Eq. (<ref>) differs from the ideal MHD solution. It also shows that the approximate solution from Eq. (<ref>) requires that kλ be sufficiently large to accurately model the numerical solution. 




 §.§ 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.<cit.> 
Electrostatic effects were eliminated from our analysis in Sec. <ref> based on an argument of scale separation between the inner tearing layer and 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 the next section 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 <ref> shows a comparison of the dispersion relation from Lee, Mahajan and Hazeline <cit.>

    ω = ω_DL - √(i) 0.85 √(|k_∥^'| c_s ρ_s)√(m_e/m_i)ω_o + k·U_i/√(ω_o)

where k_∥^' = k/l_s, ρ_s = c_s/Ω_ci and ω_o = {ω_DL}. Here, ω_DL is the DL dispersion relation from Eq. (<ref>). The second term on the right side of Eq. (<ref>) 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 effects do not appear to be the cause of stabilization of the oblique modes since the approximate theory neglecting electrostatic effects, Eq. (<ref>), 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 boundary layer theory, it is not the same. 
Namely, in Eq. (<ref>), 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 k_y U_e. 
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. daug:05). 
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. <cit.> 



 



 §.§ 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 δÂ does not vary significantly over a sufficiently short distance. 
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. <cit.>
The next section 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 in Mahajan et al<cit.>, which provides the dispersion relation 

    (ω - ω_o) (1 - i Δ^'ω/π k_∥^' v_A) = i γ_DL .

Here, γ_DL = {ω_DL}, and v_A = B̅_oz/√(4πn̅_o m_i) is the Alfvén speed. 
Figure <ref> shows a comparison of the results of Eq. (<ref>) 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. (<ref>). 
Thus, the large Δ^' extension of DL theory does not appear able to explain the observed stabilization. 




§ ANALYSIS OF BOUNDARY LAYER THEORY





 §.§ Gyroradius effects
 



A small gyroradius expansion was used in Sec. <ref> to reduce the full particle orbits from the Vlasov description in Eq. (<ref>) to Eq. (<ref>). 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 <ref>d-f show that this approximation is clearly violated at the conditions being explored in this paper. 
The largest gradient scale is the current sheet width, which for the parameters chosen in Fig. <ref> ranges from λ = ρ_i to 5ρ_i. Ions, then, 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/λ∼√(m_e/m_i)≪ 1. However, use of the approximation at the scale of the inner reconnection layer faces a more stringent requirement. Figures <ref>d-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 ≃ |ω| l_s/(kv_Ts) where ω was computed from Eqs. (<ref>) and (<ref>).<cit.> 
Although Fig. <ref> 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 at these conditions, comparison with the full VMID solutions in Fig. <ref> shows that the predictions resulting from it, namely Eq. (<ref>), 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 good 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. <ref> 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 at 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 Larmour radius (FLR) broadens this layer, but does not significantly influence the total perturbed current. <cit.>
Specifically, Fig. 1 of Ref. kari:05 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. <cit.>




 §.§ Inner layer width


Boundary layer theory is predicated 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. <ref>, but also in Sec. <ref> which led to the decoupling of the electric and magnetic fields allowing one to describe the tearing mode from the single differential equation (<ref>). 
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. <ref>. 


Figures <ref>a-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. (<ref>) at these conditions. 
For the large mass ratio cases δ_i/λ≪ 1 near the center of the current sheet (small θ), but 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. (<ref>) captures the stabilization suggests that the influence of electrostatic terms is not the cause of this effect; see Sec. <ref>. 
Regarding the subsequent application of the boundary layer approximation in developing an analytic theory, specifically the DL theory from Sec. <ref>, the approximation can be tested in more detail. 

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

    T(x) = - λ^2 ∑_s 2ω_ps^2/c^2(ω - k·U_s)/k_∥ v_Tsζ_s Z_1(ζ_s).

(Here, the version of Eq. (<ref>) 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-x_s)/l_s and ω_ps≃ω_ps(x_s), leading to an approximate expression describing the inner region 

    T_DL (x) = -λ^2 ∑_s 2 ω_ps^2 |_x_sω (ω-k·U_s)/c^2 k^2 v_Ts^2 (x-x_s)/l_s Z_1 [ ω l_s/kv_Ts(x-x_s)]




Figure <ref>a-b shows a comparison of the solutions of Eqs. (<ref>) and (<ref>) for the same parameters as the λ = 2ρ_i case of Fig. <ref>c. Here, the mode frequency is computed using the DL solutions from Eqs. (<ref>) and (<ref>). 
This reiterates the breakdown of a scale separation at oblique angles, which was estimated in Fig. <ref>. 

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. (<ref>) at sufficiently oblique angles. 
The numerical solutions reveal that the growth rate rapidly diminishes at the same conditions; see Fig. <ref>. 




 §.§ Phase and the outer region 
 

Figure <ref> answers another inconsistency between the boundary layer and numerical solutions: Why is there a prominent complex (phase) component of δÂ for oblique modes (θ≠ 0) in the numerical solutions, but this does not arise in the boundary layer theory? 

First, we observe that T(x) is purely real for the parallel mode (θ = 0^∘) (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. <ref>. The figure also shows that T(x=±∞) → 0 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 δÂ will be complex (i.e., have a finite phase). This is consistent with what is observed in the numerical solutions at finite angle; see Fig. <ref>. 



In addition, Fig. <ref> shows that at finite angles T(x=- ∞) ≠ 0 if n̂_b ≠ 0. 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±∞) to the left side of Eq. (<ref>) (the outer region). 
In fact, if one does not do this, the finite value of T(x=- ∞) implies that the integral ∫ dx T(x), which is used in the constant-ψ approximation, actually diverges. 
This divergence happens to not arise in Drake and Lee's theory when integrating ∫ dx T_DL(x) as a fortuitous consequence of the local expansion of ω_ps and k_∥ (which is not valid in this case). The local expansion forces T_DL (x =±∞) → 0 even at finite n̂_b; see Fig. <ref>a-b. 




 §.§ Role of background density


Since T(x) extends to asymptotically far distances only when n_b ≠ 0, 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 δÂ 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  {δÂ} will be small compared to the real component. 
Figure <ref>c-d shows that T(x) and T_DL(x) agree quite well in this case. 
Figure <ref>c-d shows solutions of Eq (<ref>) with n̂_b=0, confirming that {δÂ} is small compared to the real component when θ is small (see also Fig. <ref>), but that it grows comparable to the real component for x>x_s at oblique angles. 

Figure <ref>a-b shows solutions of the dispersion relation computed from Eq (<ref>) with n̂_b=0. 
These are found to damp at a similar angle as the case with n̂_b=0.3 from Fig. <ref>c. 
This provides further evidence that the damping is caused by the broadening of the inner region, not simply the finite phase component. 






 §.§ Role of guide field
  



Figure <ref> shows that the narrower spectrum of unstable modes predicted by Eq. (<ref>), 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 λ/ρ_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. <ref>. 
As the ratio of guide field to asymptotic field strength, B_oy/B̅_oz, 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. <ref>, which displays the cutoff angle (where γ = 0) as a function of B_oy/B̅_oz from each theory.  








 §.§ Role of temperature ratio
 

Figure <ref> shows that the cutoff angle computed from Eq. (<ref>) depends significantly on the temperature ratio T_i/T_e. In contrast, the cutoff angle is independent of the temperature ratio in the DL model. At T_i/T_e ∼ 1, the spectrum of unstable modes from Eq. (<ref>) is significantly more narrow than from the DL. However, as T_i/T_e increases the spectrum predicted from Eq. (<ref>) broadens substantially. At a large temperature ratio (T_i/T_e = 100) the two results agree very well; see Fig. <ref>c and d. The same trend is captured in Fig. <ref>b, which shows the cutoff angle as a function of temperature ratio for three values of the ratio λ/ρ_i. The Drake and Lee limit is obtained at large temperature ratio. 
Since the scale separation between the inner region and the λ scale is broader at large T_i/T_e, 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. <ref>. 
This temperature ratio dependence is directly related to the drift stabilization. 
For a Harris sheet, the temperature ratio sets the ratio of flow velocities U_i/U_e = - T_i/T_e, while the ion velocity is fixed by the current sheet thickness U_i/v_Ti = ρ_i/λ. 
Thus, changing the temperature ratio effectively changes the electron drift speed U_e = -U_i T_e/T_i. 
For T_i/T_e ≫ 1, U_e is very small, and one returns the results of the DL theory valid in the limit U_e → 0. 
Finally, Figs. <ref>b and <ref> show that the cutoff angle is essentially independent of λ/ρ_i, when T_i>T_e, but that a significant dependence on this ratio arises when T_i<T_e. 







 §.§ 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 x_c = λarctanh (μ_c) where μ_c = tanθ_c B_oy/B̅_oz.  
We have shown that the differential equation (<ref>) 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 θ_c 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 (λ). 

The width of the inner layer (λ^') can be estimated using Eq. (<ref>). 
We approximate the stability condition as |T(x=λ^')| = c, where c is a constant that will be determined from a fit to the numerical solutions. 
Since T∝λ/|x-x_s|, 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 Z_1 → 1; see Fig. <ref>. 
Applying the approximation that ω≃ω_R ≃ k_yU_e at oblique angles and that |x-x_s|/λ≃ 1-arctanh(μ) ≃ 1 - μ, the cutoff condition is then 

    T ≃4μ^2/1-μ( 1-μ^2 + n̂_b/1-μ^2) T_e/T_i/1+T_i/T_e = c .

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

    μ_c = c/8T_i/T_e(1 + T_i/T_e) [ √(1 + 16/cT_e/T_i/(1+T_i/T_e)) -1 ]

such that instability is expected when θ≲θ_c = arctan (μ_c B̅_oz/B_oy). 
Alternatively, the critical temperature ratio at a given angle and layer thickness (λ/ρ_i) can be related to a critical electron drift speed since U_e = - T_e/T_i U_i = - T_e/T_iρ_i/λ. 

Figure <ref> shows a comparison between the cutoff angle predicted by Eq. (<ref>) and the numerical results of Eq. (<ref>). 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 with of T is measured (see Fig. <ref>). Also shown in the figures is the cutoff from the DL dispersion relation, which is simply μ = 1, leading to 

    θ_c,DL = arctan(B̅_oz/B_oy)  .
 
Figure <ref>a shows that Eq. (<ref>) captures both the large B_oy/B̅_oz scaling and the value at which the deviation from this scaling occurs (for B_oy/B̅_oz < 1). Importantly, Fig. <ref>b shows that Eq. (<ref>) accurately captures the temperature ratio dependence of the cutoff angle (when T_i/T_e is sufficiently large).
Equation (<ref>) shows that the thickness of the inner region scales as 1/(T_i/T_e)^2 when T_i/T_e ≫ 1. 
That is, the inner region becomes very thin at 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 high temperature ratio; see Fig. <ref>b. 
For T_i/T_e ≪ 1, a dependence on the layer thickness is observed that is not captured by Eq. (<ref>).
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. (<ref>) provides a reliable estimate for the cutoff angle of oblique collisionless tearing modes. 






§ 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. <cit.> 
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·U 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. (<ref>), 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. <cit.>

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 even prominent, contributor to the the dispersion relation [arising in the plasma dispersion function of Eq. (<ref>)] in the kinetic theory described in this work. 
Yet, two fluid theories have found success in modeling collisionless tearing modes in some geometries. <cit.>
In one case the predicted growth rates differ only by a factor of √(π). <cit.> 
It would be useful to understand if two-fluid theory captures the stabilization of oblique modes observed in this work, and if that would lead to 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. 



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 number DE-SC0016159, and by 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.






99


furt:63 H. P. Furth, J. Killeen, and M. N. Rosenbluth, Phys. Fluids 6, 459 (1963).


hegn:93 C. C. Hegna, J. D. Callen, T. A. Gianakon, W. X. Qu, A. I. Smolyakov, and J. P. Wang, Plasma Phys. Control. Fusion 35, 987 (1993). 


doer:11 H. Doerk, F. Jenko, M. J. Pueschel, and D. R. Hatch, Phys. Rev. Lett. 106, 155003 (2011).


gutt:11 W. Guttenfelder, J. Candy, S. M. Kaye, W. M. Nevins, E. Wang, R. E. Bell, G. W. Hammett, B. P. LeBlanc, D. R. Mikkelsen, and H. Yuh, Phys. Rev. Lett. 106, 155004 (2011). 


call:77 J. D. Callen, Phys. Rev. Lett. 39, 1540 (1977).


biew:03 T. M. Biewer, C. B. Forest, J. K. Anderson, G. Fiksel, B. Hudson, S. C. Prager, J. S. Sarff, J. C. Wright, D. L. Bower, W. X. Ding, and S. D. Terry, Phys. Rev. Lett. 91, 045004 (2003). 


cole:15 A. J. Cole, J. M. Finn, C. C. Hegna, and P. W. Terry, Phys. Plasmas 22, 102514 (2015). 


butt:00 R. J. Buttery, G. Günter, G. Giruzzi, T. C. Hender, D. Howell, G. Huysmans, R. J. La Haye, M. Maraschek, H. Reimerdes, O. Sauter, C. D. Warrick, H. R. Wilson, H. Zohm, Plasma Phys. Control. Fusion 42, B61 (2000). 


daug:05 W. Daughton, and H. Karimabadi, J. Geophys. Res. 110, A03217 (2005).


inne:15 D. E. Innes, L.-J. Guo, Y.-M. Huang, and A. Bhattacharjee, Astrophysical Journal 813, 86 (2015). 


kari:05 H. Karimabadi, W. Daughton, and K. B. Quest, J. Geophys. Res. 110, A03214 (2005). 


kari:05b H. Karimabadi, W. Daughton, and K. B. Quest, J. Geophys. Res. 110, A03213 (2005). 


lour:07 N. F. Loureiro, A. A. Schekochihin, and S. C. Cowley, Phys. Plasmas 14, 100703 (2007).


bhat:09 A. Bhattacharjee, Y.-M. Huang, H. Yang and B. Rogers, Phys. Plasmas 16, 112102 (2009).


cass:09 P. A. Cassak, M. A. Shay and J. F. Drake, Phys. Plasmas 16, 120702 (2009). 


lour:12 N. F. Loureiro, R. Samtaney, A. A. Schikochihin and D. A. Uzdensky, Phys. Plasmas 19, 042303 (2012). 


huan:12 Y.-M. Huang and A. Bhattacharjee, Phys. Rev. Lett. 109, 265002 (2012). 


comi:16 L. Comisso and D. Grasso, Phys. Plasmas 23, 032111 (2016).


fitz:04 R. Fitzpatrick and F. Porcelli, Phys. Plasmas 11, 4713 (2004).


mirn:04 V. V. Mirnov, C. C. Hegna and S. C. Prager, Phys. Plasmas 11, 4468 (2004). 


baal:11 S. D. Baalrud, A. Bhattacharjee, Y.-M. Huang and K. Germaschewski, Phys. Plasmas 18, 092108 (2011). 


huan:11 Y.-M. Huang, A. Bhattacharjee and B. P. Sullivan, Phys. Plasmas 18, 072109 (2011). 


jara:16 J. Jara-Almonte, H. Ji, M. Yamada, J. Yoo, and W. Fox, Phys. Rev. Lett. 117, 095001 (2016).


drak:77 J. F. Drake and Y. C. Lee, Phys. Fluids 20, 1341 (1977).


gale:86 A. A. Galeev, M. M. Kuznetsova and L. M. Zeleny, Space Sci. Rev. 44, 1 (1986).


daug:09 W. Daughton, V. Roytershteyn, B. J. Albright, H. Karimabadi, L. Yin and K. J. Bowers, Phys. Rev. Lett. 103, 065004 (2009).


roge:07 B. N. Rogers, S. Kobayashi, P. Ricci, W. Dorland, J. Drake and T. Tatsuno, Phys. Plasmas 14, 092110 (2007).


numa:11 R. Numata, W. Dorland, G. G. Howes, N. F. Loureiro, B. N. Rogers and T. Tatsuno, Phys. Plasmas 18, 112106 (2011).


wang:90 X. Wang, A. Bhattacharjee and A. T. Y. Lui, J. Geophys. Res. 95, 15047 (1990).


ques:81 K. B. Quest and F. V. Coroniti, J. Geophys. Res. 86, 3289 (1981).


ques:81b K. B. Quest and F. V. Coroniti, J. Geophys. Res. 86, 3299 (1981). 


ques:85a K. B. Quest and F. V. Coroniti, J. Geophys. Res. 90, 1458 (1985).


harr:62 E. G. Harris, Nuovo Cimento. 23, 115 (1962).


daug:11 W. Daughton, V. Roytershteyn, H. Karimabadi, L. Yin, B. J. Albright, B. Bergen and K. J. Bowers, Nat. Phys. 7, 539 (2011). 


fuji:12 K. Fujimoto and R. D. Sydora, Phys. Rev. Lett. 109, 265004 (2012). 


kari:13 H. Karimabadi, and A. Lazarian, Phys. Plasmas 20, 112102 (2013). 


drak:06 J. F. Drake, M. Swisdak, H. Che and M. A. Shay, Nature 443, 553 (2006).


chen:08 L.-J. Chen, A. Bhattacharjee, P. A. Puhl-Quinn, H. Yang, N. Bessho, S. Imada, S. Mühlbachler, P. W. Daly, B. Lefebvre, Y. Khotyaintsev, A. Vaivads, A. Fazakerley and E. Georgescu, Nat. Phys. 4, 19 (2008).


liu:13 Y.-H. Liu, W. Daughton, H. Karimabadi, H. Li, and V. Roytershteyn, Phys. Rev. Lett. 110, 265004 (2013).


haze:75 R. D. Hazeltine, D. Dobrott and T. S. Wang, Phys. Fluids 18, 1778 (1975).


haze:78a R. D. Hazeltine and H. R. Strauss, Phys. Fluids 21, 1007 (1978).


haze:78b R. D. Hazeltine and D. W. Ross, Phys. Fluids 21, 1140 (1978).


maha:78 S. M. Mahajan, R. D. Hazeltine, H. R. Strauss and D. W. Ross, Phys. Rev. Lett. 41, 1375 (1978).


drak:78 J. F. Drake, Phys. Fluids 21, 1777 (1978).


lee:80 X. S. Lee, S. M. Mahajan, and R. D. Hazeltine, Phys. Fluids 23, 599 (1980). 


buss:78 M. N. Bussac, D. Edery, R. Pellat, and J. L. Soule, Phys. Rev. Lett. 40, 1500 (1978).


copp:79 B. Coppi, J. W.-K. Mark, L. Sugiyama, and G. Bertin, Phys. Rev. Lett. 42, 1058 (1979). 


hosh:87 M. Hoshino, J. Geophys. Res. 92, 7368 (1987).


maha:79 S. M. Mahajan, R. D. Hazeltine, H. R. Strauss and D. W. Ross, Phys. Fluids 22, 2147 (1979). 


daug:99 W. Daughton, Phys. Plasmas 6, 1329 (1999).


daug:03 W. Daughton, Phys. Plasmas 10, 3103 (2003). 


baal:12 S. D. Baalrud, A. Bhattacharjee and Y.-M. Huang, Phys. Plasmas, 19, 022101 (2012).


akca:16 C. Akcay, W. Daughton, V. S. Lukin, and Y.-H. Liu, Phys. Plasmas 23, 012112 (2016). 


cowl:86 S. C. Cowley, R. M. Kulsrud and T. S. Hahm, Phys. Fluids 29, 3230 (1986).


frie:61 B. D. Fried and S. C. Conte, The Plasma Dispersion Function (Academic Press, New York, 1961).


yin:08 L. Yin, W. Daughton, H. Karimabadi, B. J. Albright, K. J. Bowers and J. Margulies, Phys. Rev. Lett. 101, 125001 (2008). 


kuzn:85 M. M. Kuznetsova and L. M. Zeleny, Plasma Phys. Controlled Fusion 27, 363 (1985).


glad:90 N. T. Gladd, J. Geophys. Res. 95, 20,889 (1990)


note1:gale Ref. gale:86 takes b̂^s to be the ∥ direction in defining δÂ_∥. However, the k_∥ that is cited must be parallel to  B_o since k·b̂^s = 0. 



foot:dpg The factor of 2 appearing in front of Eq. (<ref>) is missing in Ref. gale:86 (see equation 14), but is required to get the θ = 0 results correct.


wang:93 X. Wang and A. Bhattacharjee, J. Geophys. Res. 98, 19419 (1993).
























































































