Drift instabilities in thin current sheets using a two-fluid model with   pressure tensor effects

J. Ng; A. Hakim; J. Juno; A. Bhattacharjee

arXiv:1903.09618·physics.space-ph·June 26, 2019

Drift instabilities in thin current sheets using a two-fluid model with pressure tensor effects

J. Ng, A. Hakim, J. Juno, A. Bhattacharjee

PDF

TL;DR

This paper investigates drift instabilities in thin current sheets using a two-fluid ten moment model that includes pressure tensor effects, comparing results with kinetic simulations to improve global magnetospheric modeling.

Contribution

It demonstrates the effectiveness of a two-fluid ten moment model with pressure tensor effects in capturing drift instabilities, aligning with kinetic results for certain modes.

Findings

01

Good agreement with kinetic calculations for kink mode

02

Electromagnetic and electrostatic lower hybrid instabilities require Landau fluid closure

03

Implications for global magnetospheric simulations discussed

Abstract

The integration of kinetic effects in fluid models is important for global simulations of the Earth's magnetosphere. We use a two-fluid ten moment model, which includes the pressure tensor and has been used to study reconnection, to study the drift kink and lower hybrid drift instabilities. Using a nonlocal linear eigenmode analysis, we find that for the kink mode, the ten moment model shows good agreement with kinetic calculations with the same closure model used in reconnection simulations, while the electromagnetic and electrostatic lower hybrid instabilities require modeling the effects of the ion resonance using a Landau fluid closure. Comparisons with kinetic simulations and the implications of the results for global magnetospheric simulations are discussed.

Equations42

\frac{\partial n}{\partial t} + \frac{\partial}{\partial x _{j}} (n u_{j}) m \frac{\partial}{\partial t} (n u_{i}) + \frac{\partial P _{ij}}{\partial x _{j}} = 0 = n q (E_{j} + ϵ_{ij k} u_{j} B_{k}) .

\frac{\partial n}{\partial t} + \frac{\partial}{\partial x _{j}} (n u_{j}) m \frac{\partial}{\partial t} (n u_{i}) + \frac{\partial P _{ij}}{\partial x _{j}} = 0 = n q (E_{j} + ϵ_{ij k} u_{j} B_{k}) .

P_{ij} \equiv m \int v_{i} v_{j} f d^{3} v

P_{ij} \equiv m \int v_{i} v_{j} f d^{3} v

\frac{\partial E}{\partial t} + \frac{\partial}{\partial x _{j}} (u_{j} (p + E)) = n q u_{i} E_{i} .

\frac{\partial E}{\partial t} + \frac{\partial}{\partial x _{j}} (u_{j} (p + E)) = n q u_{i} E_{i} .

\frac{\partial P _{ij}}{\partial t} + \frac{\partial Q _{ij k}}{\partial x _{k}} = n q u_{[i} E_{j]} + \frac{q}{m} ϵ_{[ik l} P_{k j]} B_{l},

\frac{\partial P _{ij}}{\partial t} + \frac{\partial Q _{ij k}}{\partial x _{k}} = n q u_{[i} E_{j]} + \frac{q}{m} ϵ_{[ik l} P_{k j]} B_{l},

Q_{ij k} \equiv m \int v_{i} v_{j} v_{k} f d^{3} v,

Q_{ij k} \equiv m \int v_{i} v_{j} v_{k} f d^{3} v,

Q_{ij k} = Q_{ij k} + u_{[i} P_{j k]} - 2 mn u_{i} u_{j} u_{k} .

Q_{ij k} = Q_{ij k} + u_{[i} P_{j k]} - 2 mn u_{i} u_{j} u_{k} .

q_{ij k} (x) = n (x) \overset{q}{^}_{ij k} (x)

q_{ij k} (x) = n (x) \overset{q}{^}_{ij k} (x)

\tilde{q}_{ij k} = - i \frac{v _{t}}{∣ k ∣} χ k_{[i} \tilde{T}_{j k]} .

\tilde{q}_{ij k} = - i \frac{v _{t}}{∣ k ∣} χ k_{[i} \tilde{T}_{j k]} .

B_{x} (z)

B_{x} (z)

n (z)

u_{y, s} = \frac{2 \remove c T _{s}}{q _{s} B _{0} L}

u_{y, s} = \frac{2 \remove c T _{s}}{q _{s} B _{0} L}

(ω - k_{y} u_{s}) n_{1, s} - k_{y} n_{0} v_{1 y, s} + i (n_{0}^{'} v_{1 z, s} + n_{0, s} v_{1 z, s}^{'}) = 0 (ω - k_{y} u_{s}) v_{1 y, s} - \frac{i q _{s}}{m _{s}} E_{y} - i Ω_{s} v_{1 z, s} - \frac{k _{y} P _{1, s}}{n _{0} m _{s}} = 0 (ω - k_{y} u_{s}) v_{1 z, s} - \frac{i q _{s}}{m _{s}} E_{z} + i Ω_{s} v_{1 y, s} + i \frac{q _{s} v _{s}}{m _{s}} B_{1 x} + i Ω_{s} \frac{n _{1, s}}{n _{0}} u_{s} + i \frac{P _{1, s}^{'}}{m _{s} n _{0}} = 0 (ω - k_{y} u_{s}) P_{1, s} + i v_{1 z, s} P_{0, s}^{'} - γ k_{y} P_{0} v_{1 y, s} + iγ P_{0, s} v_{1 z, s}^{'} = 0 ω B_{1 x} - (k_{y} E_{z} + i E_{y}^{'}) = 0 ω E_{y} - i B_{1 x}^{'} + i s \sum q_{s} (n_{0} v_{1 y, s} + n_{1 s} v_{s}) = 0 ω E_{z} - k_{y} B_{1 x} + i s \sum q_{s} n_{0} v_{1 z, s} = 0,

(ω - k_{y} u_{s}) n_{1, s} - k_{y} n_{0} v_{1 y, s} + i (n_{0}^{'} v_{1 z, s} + n_{0, s} v_{1 z, s}^{'}) = 0 (ω - k_{y} u_{s}) v_{1 y, s} - \frac{i q _{s}}{m _{s}} E_{y} - i Ω_{s} v_{1 z, s} - \frac{k _{y} P _{1, s}}{n _{0} m _{s}} = 0 (ω - k_{y} u_{s}) v_{1 z, s} - \frac{i q _{s}}{m _{s}} E_{z} + i Ω_{s} v_{1 y, s} + i \frac{q _{s} v _{s}}{m _{s}} B_{1 x} + i Ω_{s} \frac{n _{1, s}}{n _{0}} u_{s} + i \frac{P _{1, s}^{'}}{m _{s} n _{0}} = 0 (ω - k_{y} u_{s}) P_{1, s} + i v_{1 z, s} P_{0, s}^{'} - γ k_{y} P_{0} v_{1 y, s} + iγ P_{0, s} v_{1 z, s}^{'} = 0 ω B_{1 x} - (k_{y} E_{z} + i E_{y}^{'}) = 0 ω E_{y} - i B_{1 x}^{'} + i s \sum q_{s} (n_{0} v_{1 y, s} + n_{1 s} v_{s}) = 0 ω E_{z} - k_{y} B_{1 x} + i s \sum q_{s} n_{0} v_{1 z, s} = 0,

(ω - k_{y} u_{s}) n_{1, s} - k_{y} n_{0} v_{1 y, s} + i (n_{0}^{'} v_{1 z, s} + n_{0, s} v_{1 z, s}^{'}) = 0 (ω - k_{y} u_{s}) v_{1 y, s} - \frac{i q _{s}}{m _{s}} E_{y} - i Ω_{s} v_{1 z, s} - \frac{k _{y} P _{1 y y, s}}{n _{0} m _{s}} + i \frac{P _{1 y z, s}^{'}}{n _{0} m _{s}} = 0 (ω - k_{y} u_{s}) v_{1 z, s} - \frac{i q _{s}}{m _{s}} E_{z} + i Ω_{s} v_{1 y, s} + i \frac{q _{s} u _{s}}{m _{s}} B_{1 x} + i Ω_{s} \frac{n _{1, s}}{n _{0}} u_{s} - k_{y} \frac{P _{1 y z, s}}{n _{0} m _{s}} + i \frac{P _{1 z z, s}^{'}}{m _{s} n _{0}} = 0 (ω - k_{y} u_{s}) P_{1 xx, s} - k_{y} P_{0, s} v_{1 y, s} + i P_{0, s} v_{1 z, s}^{'} + i P_{0, s}^{'} v_{1 z, s} + i ∣ k_{0, s} ∣ v_{t, s} (P_{1 xx, s} - \frac{P _{1, s}}{3}) = 0 (ω - k_{y} u_{s}) P_{1 y y, s} + i P_{0, s} (i k_{y} v_{1 y, s} + v_{1 z, s}^{'}) - 2 k_{y} P_{0, s} v_{1 y, s} + i v_{1 z, s} P_{0, s}^{'} - 2 i Ω_{j} P_{1 y z, j} + i ∣ k_{0, s} ∣ v_{t, s} (P_{1 y y, s} - \frac{P _{1, s}}{3}) = 0 (ω - k_{y} u_{s}) P_{1 y z, s} + i P_{0, s} v_{1 y, s}^{'} - k_{y} P_{0, s} v_{1 z, s} - i Ω_{s} (P_{1 z z, s} - P_{1 y y, s}) + i ∣ k_{0, s} ∣ v_{t, s} P_{1 y z, s} = 0 (ω - k_{y} u_{s}) P_{1 z z, s} + i P_{0, s} (i k_{y} v_{1 y, s} + v_{1 z, s}^{'}) + 2 i P_{0, s} v_{1 z, s}^{'} + i v_{1 z, s} P_{0, s}^{'} + 2 i Ω_{s} P_{1 y z, s} + i ∣ k_{0, s} ∣ v_{t, s} (P_{1 z z, s} - \frac{P _{1, s}}{3}) = 0,

(ω - k_{y} u_{s}) n_{1, s} - k_{y} n_{0} v_{1 y, s} + i (n_{0}^{'} v_{1 z, s} + n_{0, s} v_{1 z, s}^{'}) = 0 (ω - k_{y} u_{s}) v_{1 y, s} - \frac{i q _{s}}{m _{s}} E_{y} - i Ω_{s} v_{1 z, s} - \frac{k _{y} P _{1 y y, s}}{n _{0} m _{s}} + i \frac{P _{1 y z, s}^{'}}{n _{0} m _{s}} = 0 (ω - k_{y} u_{s}) v_{1 z, s} - \frac{i q _{s}}{m _{s}} E_{z} + i Ω_{s} v_{1 y, s} + i \frac{q _{s} u _{s}}{m _{s}} B_{1 x} + i Ω_{s} \frac{n _{1, s}}{n _{0}} u_{s} - k_{y} \frac{P _{1 y z, s}}{n _{0} m _{s}} + i \frac{P _{1 z z, s}^{'}}{m _{s} n _{0}} = 0 (ω - k_{y} u_{s}) P_{1 xx, s} - k_{y} P_{0, s} v_{1 y, s} + i P_{0, s} v_{1 z, s}^{'} + i P_{0, s}^{'} v_{1 z, s} + i ∣ k_{0, s} ∣ v_{t, s} (P_{1 xx, s} - \frac{P _{1, s}}{3}) = 0 (ω - k_{y} u_{s}) P_{1 y y, s} + i P_{0, s} (i k_{y} v_{1 y, s} + v_{1 z, s}^{'}) - 2 k_{y} P_{0, s} v_{1 y, s} + i v_{1 z, s} P_{0, s}^{'} - 2 i Ω_{j} P_{1 y z, j} + i ∣ k_{0, s} ∣ v_{t, s} (P_{1 y y, s} - \frac{P _{1, s}}{3}) = 0 (ω - k_{y} u_{s}) P_{1 y z, s} + i P_{0, s} v_{1 y, s}^{'} - k_{y} P_{0, s} v_{1 z, s} - i Ω_{s} (P_{1 z z, s} - P_{1 y y, s}) + i ∣ k_{0, s} ∣ v_{t, s} P_{1 y z, s} = 0 (ω - k_{y} u_{s}) P_{1 z z, s} + i P_{0, s} (i k_{y} v_{1 y, s} + v_{1 z, s}^{'}) + 2 i P_{0, s} v_{1 z, s}^{'} + i v_{1 z, s} P_{0, s}^{'} + 2 i Ω_{s} P_{1 y z, s} + i ∣ k_{0, s} ∣ v_{t, s} (P_{1 z z, s} - \frac{P _{1, s}}{3}) = 0,

i (\nabla \cdot q)_{xx} i (\nabla \cdot q)_{y y} i (\nabla \cdot q)_{y z} i (\nabla \cdot q)_{z z} = i \frac{4}{9 π} k_{y} v_{t} (P_{1, xx} - n_{1} T_{0}) = i (\frac{4}{π} k_{y} v_{t} (P_{1, y y} - n_{1} T_{0})) + \frac{2}{3} \frac{4}{π} v_{t} \frac{\partial P _{1, y z}}{\partial z} = i (\frac{2}{3} \frac{4}{π} k_{y} v_{t} P_{1, y z}) + \frac{1}{3} \frac{4}{π} v_{t} \frac{\partial ( P _{1, z z} - n _{1} T _{0} )}{\partial z} = i (\frac{4}{9 π} k_{y} v_{t} (P_{1, z z} - n_{1} T_{0}))

i (\nabla \cdot q)_{xx} i (\nabla \cdot q)_{y y} i (\nabla \cdot q)_{y z} i (\nabla \cdot q)_{z z} = i \frac{4}{9 π} k_{y} v_{t} (P_{1, xx} - n_{1} T_{0}) = i (\frac{4}{π} k_{y} v_{t} (P_{1, y y} - n_{1} T_{0})) + \frac{2}{3} \frac{4}{π} v_{t} \frac{\partial P _{1, y z}}{\partial z} = i (\frac{2}{3} \frac{4}{π} k_{y} v_{t} P_{1, y z}) + \frac{1}{3} \frac{4}{π} v_{t} \frac{\partial ( P _{1, z z} - n _{1} T _{0} )}{\partial z} = i (\frac{4}{9 π} k_{y} v_{t} (P_{1, z z} - n_{1} T_{0}))

ω F + k_{y} A_{y} F + A_{z} \frac{\partial F}{\partial z} + S F = 0,

ω F + k_{y} A_{y} F + A_{z} \frac{\partial F}{\partial z} + S F = 0,

ω_{r}

ω_{r}

γ

1 + \frac{ω _{p e}^{2}}{ω _{ce}^{2}} (1 + \frac{ω _{p e}^{2}}{c ^{2} k _{y}^{2}}) - \frac{2 ω _{p i}^{2}}{k _{y}^{2} v _{t i}^{2}} (1 + β_{i} /2) \frac{k _{y} V _{d i}}{ω + k _{y} V _{d i}} + \frac{2 ω _{p i}^{2}}{k _{y}^{2} v _{t i}^{2}} [1 + ζ_{i} Z (ζ_{i})] = 0

1 + \frac{ω _{p e}^{2}}{ω _{ce}^{2}} (1 + \frac{ω _{p e}^{2}}{c ^{2} k _{y}^{2}}) - \frac{2 ω _{p i}^{2}}{k _{y}^{2} v _{t i}^{2}} (1 + β_{i} /2) \frac{k _{y} V _{d i}}{ω + k _{y} V _{d i}} + \frac{2 ω _{p i}^{2}}{k _{y}^{2} v _{t i}^{2}} [1 + ζ_{i} Z (ζ_{i})] = 0

n_{1} = - n_{0} \frac{e ϕ _{1}}{T _{0}} R (ζ)

n_{1} = - n_{0} \frac{e ϕ _{1}}{T _{0}} R (ζ)

R_{5 m o m} (ζ) = \frac{1}{γ - 2 ζ ^{2}}

R_{5 m o m} (ζ) = \frac{1}{γ - 2 ζ ^{2}}

R_{n o n l oc a l} (ζ) = \frac{χ _{1} - i ζ}{χ _{1} - 3 i ζ - 2 χ _{1} ζ ^{2} + 2 i ζ ^{3}},

R_{n o n l oc a l} (ζ) = \frac{χ _{1} - i ζ}{χ _{1} - 3 i ζ - 2 χ _{1} ζ ^{2} + 2 i ζ ^{3}},

R_{l oc a l} (ζ) = \frac{α - i ζ}{5 α /3 - 3 i ζ - 2 α ζ ^{2} + 2 i ζ ^{3}}

R_{l oc a l} (ζ) = \frac{α - i ζ}{5 α /3 - 3 i ζ - 2 α ζ ^{2} + 2 i ζ ^{3}}

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Full text

Drift instabilities in thin current sheets using a two-fluid model with pressure tensor effects

Abstract

The integration of kinetic effects in fluid models is important for global simulations of the Earth’s magnetosphere. We use a two-fluid ten moment model, which includes the pressure tensor and has been used to study reconnection, to study the drift kink and lower hybrid drift instabilities. Using a nonlocal linear eigenmode analysis, we find that for the kink mode, the ten moment model shows good agreement with kinetic calculations with the same closure model used in reconnection simulations, while the electromagnetic and electrostatic lower hybrid instabilities require modeling the effects of the ion resonance using a Landau fluid closure. Comparisons with kinetic simulations and the implications of the results for global magnetospheric simulations are discussed.

\draftfalse\journalname

JGR: Space Physics

Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08543, USA Princeton Plasma Physics Laboratory, Princeton, NJ 08544, USA Institute for Research in Electronics and Applied Physics, University of Maryland, College Park, MD 20742, USA

{keypoints}

The drift kink and lower hybrid drift instabilities are studied using a \changeten moment modelten-moment fluid model.

Inclusion of the non-gyrotropic pressure tensor improves agreement with kinetic results for the kink mode

Ion physics of the lower hybrid drift instability can be reproduced using a nonlocal heat flux closure.

1 Introduction

Thin current sheets are often found in the Earth’s magnetosphere, and are unstable to a variety of modes, including the tearing mode, drift-kink mode and lower hybrid drift instability (LHDI).

The drift-kink mode is an ion scale mode ( $k\rho_{i}\sim 1$ ) driven by the streaming of ions and electron, and was once thought to be a possible mechanism for substorm onset, becoming the subject of theoretical and numerical studies using fluid and kinetic theory (Daughton, \APACyear1999\APACexlab\BCnt1, \APACyear1998; Pritchett \BOthers., \APACyear1996; Zhu \BBA Winglee, \APACyear1996; Yoon \BOthers., \APACyear1998; Ozaki \BOthers., \APACyear1996)\changecitation to Ozaki added. However, it was shown (Daughton, \APACyear1999\APACexlab\BCnt2) that the electron-ion drift-kink instability has a strongly reduced growth rate at the physical mass ratio. More recently, there has been work on ion-ion kink instabilities driven by the velocity difference between background and current carrying ions (Karimabadi, Daughton\BCBL \BOthers., \APACyear2003; Karimabadi, Pritchett\BCBL \BOthers., \APACyear2003).

Compared to the drift kink instability, the electrostatic lower hybrid drift instability (LHDI) has shorter wavelength, with a broad range of wavenumbers $(m_{e}/m_{i})^{1/4}<k\rho_{e}<1$ with frequency $\omega\approx\Omega_{lh}\sim\sqrt{\Omega_{ce}\Omega_{ci}}$ (Daughton, \APACyear1999\APACexlab\BCnt2; Davidson \BOthers., \APACyear1977). These fluctuations are located at the edge of the current sheet, where the density gradient is strongest, and have been observed in space, experiments and simulations (Bale \BOthers., \APACyear2002; Carter \BOthers., \APACyear2002; Lapenta \BOthers., \APACyear2003). While the electrostatic LHDI does not always enhance \removealways reconnection by itself due to its location away from the centre of the current sheet, it can alter the structure of the current layer due to its comparatively faster growth rate and drive secondary instabilities such as the drift-kink or Kelvin-Helmholtz instabilities (Price \BOthers., \APACyear2016; Lapenta \BOthers., \APACyear2003; Daughton, \APACyear2003). Additionally, there is also a longer wavelength $k\sqrt{\rho_{e}\rho_{i}}\sim 1$ electromagnetic lower hybrid mode which has a lower growth rate (Daughton, \APACyear2003). This instability can be observed at the centre of the current sheet and can influence the reconnection process (Roytershteyn \BOthers., \APACyear2012). Within the magnetosphere, there have been observations of the LHDI at both the magnetopause (Graham \BOthers., \APACyear2017) and magnetotail (Zhou \BOthers., \APACyear2009).

The Earth’s magnetosphere is comprised mostly of a collisionless plasma. Global simulations of the magnetosphere have relied so far mostly on single-fluid MHD, which is inadequate for collisionless plasmas. In recent years, we have attempted to extend fluid models to incorporate more kinetic effects in magnetospheric systems using higher moment models (Wang \BOthers., \APACyear2018). While these models have been successful in simulating large reconnecting systems (Wang \BOthers., \APACyear2015, \APACyear2018; Ng \BOthers., \APACyear2015, \APACyear2017; Allmann-Rahn \BOthers., \APACyear2018), there have not been detailed studies on how well the drift instabilities are represented by the models. Though there is some work on these instabilities in field-reversed configurations (Hakim \BBA Shumlak, \APACyear2007), existing fluid theory for current sheets (Daughton, \APACyear1999\APACexlab\BCnt1; Yoon \BOthers., \APACyear2002, \APACyear1998; Pritchett \BOthers., \APACyear1996) shows some discrepancies with kinetic theory (Davidson \BOthers., \APACyear1977; Daughton, \APACyear1999\APACexlab\BCnt2, \APACyear2003), and does not include the pressure tensor which is evolved by the ten moment model.

In the light of these attempts, it is important to understand if the extended fluid equations can model these instabilities, and if inclusion of the pressure tensor and associated closure improves the agreement between kinetic and fluid models. It is also necessary to determine if the same closures which give good agreement with kinetic studies of reconnection can simultaneously describe the instabilities. One area of interest is the growth rate of the kink and sausage modes, where the fluid calculations can have faster growth rates at shorter wavelengths (Daughton, \APACyear1999\APACexlab\BCnt1; Yoon \BOthers., \APACyear2002; Pritchett \BOthers., \APACyear1996), while in kinetic theory, the fastest growing kink mode is around $kL\sim 1$ , $L\sim\rho_{i}$ being the length scale of gradients in the equilibrium, and the sausage mode is stable (Daughton, \APACyear1999\APACexlab\BCnt2).

This paper is focused on linear eigenmode calculations of the drift instabilities in Harris sheets (Harris, \APACyear1962) using the five and ten moment models. The five moment model is a standard two fluid model with isotropic pressure, and reduces to Hall MHD in the limit of $m_{e}\to 0$ , $c\to\infty$ and $n_{i}=n_{e}$ , while the ten moment model includes the effects of an anisotropic pressure tensor and a heat flux closure. Our results show that the ten moment model is able to model the drift kink instability and magnetic reconnection simultaneously, while a proper treatment of the lower hybrid instabilities requires capturing the ion kinetic response using a Landau fluid closure for the heat flux, though the instability still appears when using a simple fluid model. The remainder of the paper is organised as follows: Section 2 describes the moment models and closures used in the calculations, and Section 3 describes the linear eigenmode calculations. The results of the kink and LHDI calculations are shown in Sections 4 and 5, with some discussion of the appropriate closure to use for the LHDI in Section 5. Finally, comparisons between fluid and a fully kinetic Vlasov-Maxwell simulation are presented in Section 5.2, and we conclude in Section 6.

2 Moment equations

For each species, the fluid equations are obtained by taking velocity moments of the Vlasov equations. This leads to

[TABLE]

where $\mathcal{P}_{ij}$ is the second moment of the distribution function

[TABLE]

In the five moment model, the pressure is assumed to be isotropic, and we evolve the energy equation in addition to the continuity and momentum equations.

[TABLE]

Here $\mathcal{E}=\tfrac{1}{2}nmu^{2}+nm\epsilon$ , where $\epsilon=P/[(\gamma-1)nm]$ is the internal energy per unit mass. For this paper we use $\gamma=5/3$ .

The ten moment model evolves the full pressure tensor according to

[TABLE]

where $\mathcal{Q}_{ijk}$ is the third moment of the distribution function

[TABLE]

and the square brackets denote a sum over permutations of the indices (e.g. $u_{[i}E_{j]}=u_{i}E_{j}+u_{j}E_{i}$ ). Following (Wang \BOthers., \APACyear2015) one can write $\mathcal{Q}_{ijk}$ in terms of the heat flux tensor $Q_{ijk}\equiv m\int(v_{i}-u_{i})(v_{j}-u_{j})(v_{k}-u_{k})fd^{3}v$

[TABLE]

For collisionless plasmas in the unmagnetised limit, we use a three-dimensional extension of the Hammett-Perkins closure, which can be expressed as follows for both electrons and ions (Hammett \BBA Perkins, \APACyear1990):

[TABLE]

where $\hat{q}_{ijk}$ in Fourier space is $\tilde{q}_{ijk}$ and is calculated as

[TABLE]

Here $\tilde{T}_{jk}$ is the Fourier transform of the deviation of the local temperature tensor from the mean. The $1/|k|$ scaling makes this a non-local closure when expressed in real space (Hammett \BOthers., \APACyear1992; Snyder \BBA Hammett, \APACyear2001) and provides a 1 to 3 pole Padé approximation of various components of the dielectric tensor. The coefficient $\chi=\sqrt{4/9\pi}$ is the best fit value for the diagonal $q_{iii}$ component and reduces to the closure in (Hammett \BBA Perkins, \APACyear1990; Hammett \BOthers., \APACyear1992) in the 1-D limit. This closure has been used to study reconnection in the context of magnetic island coalescence, and gives better agreement with kinetic results than Hall MHD (Ng \BOthers., \APACyear2017).

Due to the computational costs involved in calculating the nonlocal heat flux, relaxation of the pressure tensor to local isotropy is a more common approximation, and has been used successfully in large scale studies of reconnection and magnetospheres (Wang \BOthers., \APACyear2015; Ng \BOthers., \APACyear2015; Wang \BOthers., \APACyear2018). With this model the heat flux divergence term is replaced by $\partial_{i}Q_{ijk}=v_{t}|k_{0}|(P_{ij}-P\delta_{ij})$ (Wang \BOthers., \APACyear2015; Ng \BOthers., \APACyear2015; Hesse \BOthers., \APACyear1995; Yin \BOthers., \APACyear2001), where $v_{t}=\sqrt{2T/m}$ is the thermal velocity of the associated species and $k_{0}$ is a free parameter for each species.

As this work is focused on understanding if the drift instabilities exist within the ten-moment model and whether they will be present in global simulations and interact with reconnection, we study both local relaxation and the nonlocal closure over a variety of parameter regimes.

3 Eigenmode calculations

To study the current sheet instabilities, we begin with the exact Harris equilibrium (Harris, \APACyear1962). The magnetic field and density are described by

[TABLE]

with species drift velocities

[TABLE]

The temperature is determined by the equilibrium condition \change $\beta=1$$\beta_{e}+\beta_{i}=1$ . Here $\beta_{s}$ is the species plasma beta defined as $2\mu_{0}n_{0}T_{s}/B_{0}^{2}$ .

We consider perturbations about the equilibrium in the form $f(y,z,t)=f_{1}(z)\exp\left(i\left(k_{y}y-\omega t\right)\right)$ , with no variation in the $x$ direction (parallel to the equilibrium magnetic field). This is orthogonal to the usual $2$ -D plane used in reconnection studies. For the modes we are studying, the perturbed quantities $B_{y},B_{z},E_{x}$ and $v_{x}$ are identically zero (Pritchett \BOthers., \APACyear1996), and in the ten-moment model, the pressure tensor components $P_{xy}$ and $P_{xz}$ are also zero. This leads to reduced systems of $11$ and $17$ equations for the five and ten-moment models respectively (for two species).

For the five moment equations, they are (normalised to $c=1$ , $d_{i0}=1$ , $\omega_{pi0}=1$ \addin simulation units):

[TABLE]

where the primes represent $z$ derivatives and $\Omega_{s}=q_{s}B_{0}(z)/m_{s}$ .

The linear ten-moment equations are as follows:

[TABLE]

where $P_{1}=P_{1xx}+P_{1yy}+P_{1zz}$ is the perturbed trace of the pressure tensor and Maxwell’s equations remainin the same. The modifications are the additional equations for the pressure tensor components and the replacement of the pressure gradient by the divergence of the pressure tensor in the momentum equations. The terms proportional to $|k_{0,s}|$ in the pressure tensor evolution represent the local isotropisation discussed in Section 2.

When using the nonlocal closure, we replace the relaxation terms in Eq. (LABEL:eq:tenmom) with the following expressions for the nonlocal heat flux

[TABLE]

Here we have only kept the $k_{y}$ terms in $q_{ijk}\propto k_{[i}T_{jk]}/|k|$ .

Although it is possible to reduce the five-moment system to a single second-order differential equation which is amenable for analysis (Yoon \BOthers., \APACyear2002), the additional equations in the ten-moment system make it somewhat difficult to use the same method. Instead, we note that the equations can be written as

[TABLE]

where $\mathcal{A}_{y},\mathcal{A}_{z}$ and $\mathcal{S}$ are coefficient matrices. The instabilities of the system can then be found directly by discretizing the equations and solving for the eigenvalues of the resulting matrix. In this work we used 6th order central differences to calculate the $z$ derivatives. The equations are solved from $z=-12.8L$ to $z=12.8L$ . The resolution of kink modes and longer wavelength lower hybrid modes typically requires fewer than 250 grid points. For shorter wavelength lower hybrid modes, which are more localised and can have finer structure, we use a smaller domain $z=-6.8L$ to $z=6.8L$ , and 250 points. Once specific eigenvalues are found, convergence is tested by increasing resolution by a factor of four and using a sparse solver to find the closest solutions to the selected eigenvalue.

One feature of this method compared to the search methods employed by (Daughton, \APACyear1999\APACexlab\BCnt2; Yoon \BOthers., \APACyear2002) is that we find all the modes of the system (limited by the resolution and numerical method), and post-processing is necessary to identify the unstable modes of interest.

4 Drift-kink instability

The solution of Eq. (15) for both systems leads to a spectrum of eigenmodes over a range of $k_{y}$ . In this – and the following – section we compare the five and ten moment solutions for the drift-kink and lower hybrid instabilities. Where possible, we use similar parameters to the kinetic calculations in the literature (Daughton, \APACyear1999\APACexlab\BCnt2, \APACyear2003).

We begin by studying the case of an electron-positron plasma, $m_{i}/m_{e}=1,T_{i}/T_{e}=1,\rho_{i}/L=0.5,v_{t,e}=0.25c$ . This particular set of parameters has been studied in earlier work (Daughton, \APACyear1999\APACexlab\BCnt1; Pritchett \BOthers., \APACyear1996) and is a useful basis for direct comparison. Figure 1 shows the fastest growing kink modes for this configuration. \changeThe five moment results are comparable to those ofThe variation of the growth rate with $k_{y}$ shown in Fig. 1 \addshows good agreement with the results of (Pritchett \BOthers., \APACyear1996), with a maximum growth rate of $\gamma/\Omega_{ci}=0.22$ , while the ten-moment result shows a maximum of $\gamma/\Omega_{ci}=0.17$ at a longer wavelength with $k_{y}L\approx 1$ . This is in better agreement with the linear Vlasov results in (Daughton, \APACyear1998, \APACyear1999\APACexlab\BCnt1). With the ten-moment model, there is a plateau for $k_{y}L>1.5$ , which is sensitive to the value of $k_{0,i}$ used. For this set of results we used \addlocal relaxation with $k_{0,e}=1/d_{e}$ , $k_{0,i}=1/d_{i}$ , a choice similar to that used in earlier reconnection studies (Ng \BOthers., \APACyear2015; Wang \BOthers., \APACyear2015).

At long wavelengths, both models approach the dashed lines, which show the incompressible solution (Daughton, \APACyear1999\APACexlab\BCnt1)

[TABLE]

The equation systems (LABEL:eq:fivemom) and (LABEL:eq:tenmom) support a spectrum of eigenmodes. In Fig. 2 we show the mode structure of unstable odd and even harmonics for both five and ten moment models at a fixed wavenumber $k_{y}L=0.5$ . Other physical parameters are $m_{i}/m_{e}=64,\rho_{i}/L=0.7$ and $T_{i}/T_{e}=1$ . In the left column, the five moment eigenfunctions are shown, with both odd and even (kink and sausage) modes supported by the system. The real frequencies are consistent with the ion diamagnetic frequency, with $k_{y}u_{i}/\Omega_{ci0}=0.245$ . In the right column, the ten moment eigenmodes are shown. We were only able to find a single kink mode growing at a similar growth rate to the five moment solutions, with the sausage mode growth rate more than a factor of three smaller.

The scaling of the models with physical parameters is shown in Figs. 3 and 4. In Fig. 3 the scaling of the growth rates and frequencies of the kink and sausage modes with the ratio of ion and electron temperatures is shown. Here we are comparing modes with structure similar to those shown in Fig. 2. In both models, the growth rates increase as $T_{i}/T_{e}$ decreases, in agreement with kinetic calculations (Daughton, \APACyear1999\APACexlab\BCnt2). The differences between the models are evident in the sausage mode growth rates, where the five moment model shows a sausage mode growing at almost the same rate as the kink mode, similar to (Pritchett \BOthers., \APACyear1996), while the sausage mode in the ten-moment model grows 3 to 4 times more slowly than the kink mode across the range of temperatures.

In large scale simulations, the use of a reduced ion/electron mass ratio is common in order to reduce computational costs. It is thus important to understand how the instabilities scale with $m_{i}/m_{e}$ to ensure that the reduced models do not excite unrealistic instabilities, particularly since the kink instability can potentially disrupt current sheets. In kinetic theory, it is known that the drift-kink instability growth rate is greatly reduced at higher $m_{i}/m_{e}$ (Daughton, \APACyear1999\APACexlab\BCnt2, \APACyear1998). Fig. 4 shows how the two models scale with mass ratio. The five moment model shows an increase in growth rate with mass ratio, with the maximum growth rate being found at shorter wavelengths. In contrast, the kink mode in the ten-moment model shows a decrease in the growth rate as $m_{i}/m_{e}$ increases, with the fastest growth rate occurring at $k_{y}\rho_{i}\sim 1$ , in agreement with kinetic results. The differences between these models show the importance of keeping the non-isotropic pressure tensor in modeling the kink instability.

4.1 Scaling with relaxation parameters

The local ten-moment model we use has free parameters, the relaxation constants $k_{0,s}$ for the different species. In previous studies (Wang \BOthers., \APACyear2015; Ng \BOthers., \APACyear2015; Wang \BOthers., \APACyear2018), it was found that setting $k_{0,s}\sim 1/d_{s}$ was suitable for modeling magnetic reconnection. It is thus important to understand how the kink instabilities are affected by different $k_{0}$ and if the values used for reconnection are suitable for studying these instabilities.

We perform two scaling studies, one in which we hold the ion relaxation parameter constant at $1/d_{i}$ , and one in which we hold the electron parameter constant at $1/d_{e}$ . \addThe mass ratio remains $m_{i}/m_{e}=64$ . The results are shown in Fig. 5. For the kink instability, the variation of $k_{0,i}$ has a greater effect on the maximum growth rate. As $k_{0,i}$ is increased from $1/d_{i}$ to $100/d_{i}$ , the ions are isotropised and the fluid model for ions more closely resembles the five moment model, with an increase in growth rate and a shift of the fastest growing mode to longer wavelength. Decreasing the value of $k_{0,i}$ has a small impact on the growth rate, with an increase of $<0.01\Omega_{ci0}$ over two orders of magnitude. The effect of the electron relaxation parameter is comparatively small, with an increase in growth rate at smaller $k_{0,e}$ \change,. \changeindicatingThese results indicate that retaining the additional ion physics is sufficient for describing the kink mode.

5 Lower hybrid drift instability

The equation systems (LABEL:eq:fivemom) and (LABEL:eq:tenmom) also support the lower hybrid drift instability (LHDI) (Daughton, \APACyear2003; Yoon \BOthers., \APACyear2002). These instabilities can be found at either $k_{y}(\rho_{i}\rho_{e})^{1/2}\sim 1$ or $k_{y}\rho_{e}\sim 1$ . The shorter wavelength modes ( $k_{y}\rho_{e}\sim 1$ ) have frequency on order of $\omega_{lh}\approx\sqrt{\Omega_{ci}\Omega_{ce}}$ (Davidson \BOthers., \APACyear1977; Daughton, \APACyear2003) and are localised around the edge of the current sheet, while the longer wavelength modes have a lower frequency and can penetrate to the centre of the current sheet (Daughton, \APACyear2003).

We first review the local kinetic theory of the LHDI in order to guide our understanding of how to approximate the LHDI using fluid models. In the cold electron limit, for $k_{y}\rho_{e}\gg 1$ modes, the local dispersion relation of the LHDI can be written as (Davidson \BOthers., \APACyear1977)

[TABLE]

where $Z(\zeta_{i})$ is the ion plasma dispersion function, $\zeta_{i}=\omega/k_{y}v_{ti}$ and $V_{di}$ is the ion diamagnetic drift velocity. Note that this is the dispersion relation in the ion rest frame, so any comparisons with our Eqs (LABEL:eq:fivemom) and (LABEL:eq:tenmom) should be Doppler shifted. In this limit the ion kinetic response is assumed to be unmagnetised, and gradients of perturbed quantities in the $z$ direction are neglected.

The dispersion relation has known unstable solutions in the \changefluidadiabatic ( $\zeta_{i}\gg 1$ ) and kinetic ( $\zeta_{i}\ll 1$ ) \noteadded subscript limits (Hirose \BBA Alexeff, \APACyear1972), through coupling of the drift and lower hybrid wave or the ion resonance. Because the ions can be treated as unmagnetised, the nonlocal closure of (Hammett \BBA Perkins, \APACyear1990) or 3-d generalisation of (Ng \BOthers., \APACyear2017) would be the best fluid model for capturing the kinetic ion physics. A discussion of how well the fluid models approximate $1+\zeta_{i}Z(\zeta_{i})$ is in the appendix.

For the electrons, since $k_{y}\rho_{e}\sim 1$ , we do not expect the nonlocal closure to be applicable in this regime for the cross-field heat flux. In a more general situation, a gyrofluid model with finite Larmor radius effects would be appropriate (Snyder \BBA Hammett, \APACyear2001; Tassi \BOthers., \APACyear2018). However, in a reconnecting current sheet with no guide field, such models would be inapplicable close to the centre of the sheet as the magnetic field becomes close to zero. In this case, particularly in the cold electron limit where the electron dynamics affect the rate but not the instability threshold (Davidson \BOthers., \APACyear1977), a ten-moment model with reduced isotropisation may be a better description. The role of the electron closure is discussed in Section 5.1.3.

5.1 Results

In this section we study the electrostatic and electromagnetic LHDI using the five moment model, the ten moment local model with the closure used in reconnection models ( $k_{0,s}=1/d_{s}$ ), the local model with $k_{0,i}=k_{y}$ and no electron relaxation, and the ten moment model with a nonlocal closure for ions and no electron relaxation. We use two current sheets, a thicker sheet with $\rho_{i}=L$ and hotter ions $T_{i}/T_{e}=5$ , where Eq. (18) would be most applicable, and a current sheet with the parameters in (Daughton, \APACyear2003) for direct comparison with kinetic work.

5.1.1 Electrostatic LHDI

We first present the results of calculations for the electrostatic LHDI. The parameters used here are $m_{i}/m_{e}=256$ , $T_{i}/T_{e}=5$ , $\omega_{pe}/\Omega_{ce}=5$ and $\rho_{i}/L=1$ . Modes are calculated for $k_{y}\rho_{e}=0.5$ . \addThe current sheets support a spectrum of unstable modes (“harmonics” in (Daughton, \APACyear2003)) \addand we show four for each model and configuration we study in the figures.

The eigenmodes are shown in Fig. 6, where we use the \removefive moment model, the local ten moment model with $k_{0,e}=0,k_{0,i}=k_{y}$ and the ten moment model with nonlocal ions and $k_{0,e}=0$ . For this set of parameters, the \addfive-moment model and ten\add-moment model with $k_{0,s}=1/d_{s}$ \changewasare stable to lower hybrid instabilities.

\change

In this instance, the five moment model is comparatively stable, with the growth rate being $5$ to $10$ times smaller than calculated using the nonlocal ten moment model. ThisIn this case, the stability of the five-moment model is likely due to its inability to model the ion response correctly, which \changeare beis important in thicker sheets with a smaller equilibrium drift velocity (Davidson \BOthers., \APACyear1977). Both ten moment calculations show that the LHDI is present, and the model using the nonlocal closure for the ions has a growth rate and structure that is consistent with local kinetic theory, which gives a growth rate of $\gamma/\Omega_{ci0}=2.9$ in the region around $z/L=1.5$ .

Figure 7 shows the eigenmode calculations for the second set of parameters, with $\rho_{i}/L=2$ , $T_{i}=T_{e}$ , $m_{i}/m_{e}=512$ and $\omega_{pe}/\omega_{ce}=5$ , also used in (Daughton, \APACyear2003). This is a comparatively thinner current sheet, and all the models are unstable to the LHDI, with similar mode structures but different growth rates. The five moment model has the fastest growing modes, while the ten moment models with $k_{0,e}=0$ have growth rates and frequencies differing by less than $\Omega_{ci0}$ . For similar real frequency, these modes have a faster growth rate than in the kinetic calculation of (Daughton, \APACyear2003). The model with $k_{0,s}=1/d_{s}$ shows the lowest growth rates, which is a general trend reproduced in the next sections.

5.1.2 Electromagnetic LHDI

For the longer wavelength electromagnetic LHDI, we first present the results of calculations with the thicker $\rho_{i}=L$ current sheet with the same parameters as the previous \changesecondsection, but wavelength $k_{y}\sqrt{\rho_{i}\rho_{e}}=1$ . Selected eigenmodes are shown in Fig. 8, where the local model labeled (1) again refers to the closure with $k_{0,e}=0$ , $k_{0,i}=k_{y}$ , and the model labeled (2) has $k_{0,s}=1/d_{s}$ . Here the local model with $k_{0,s}=1/d_{s}$ shows very weakly unstable modes, with structures reminiscent of the electrostatic modes, while the five moment model and other ten moment models show broader mode structures which extend to the centre of the current sheet, which is expected for these modes (Daughton, \APACyear2003).

With the thinner sheet used in (Daughton, \APACyear2003), we find the modes shown in Fig. 9. Again, the ten moment model with the parameters used in reconnection studies ( $k_{0,s}=1/d_{s}$ ) shows much more stable modes around the drift frequency $k_{y}u_{i}\approx 8.0\Omega_{ci0}$ . The five moment model, ten moment model with nonlocal ions and with $k_{0,e}=0,k_{0,i}=k_{y}$ are able to capture the electromagnetic LHDI, though the growth rates show quantitative differences with the results of (Daughton, \APACyear2003), though the frequencies and growth rates have a similar range of values. We believe the discrepancies are due to the limitations of our electron model, which will be demonstrated in the next section.

5.1.3 Sensitivity to electron model

Although the ten-moment model contains non-gyrotropic pressure effects, it is not clear how this affects the calculations of the lower hybrid instabilities where $k_{y}\rho_{e}\sim 1$ . As we did with the kink instability, we perform a scaling of the electron relaxation parameter and study how the fastest growing lower hybrid mode varies with $k_{0,e}$ . In these calculations we use the nonlocal ten moment ion model as it is the best approximation to the ion kinetic response.

The parameters of the current sheet used in this scan are $m_{i}/m_{e}=256$ , $T_{i}/T_{e}=5$ , $\omega_{pe}/\Omega_{ce}=5$ and $\rho_{i}/L=1$ , and we use $k_{y}\rho_{e}=0.5$ . The results are shown in Fig. 10, with calculations using nonlocal and five moment electrons also plotted for reference. The growth rate of the instability is larger in the limits $k_{0}\to 0$ and $k_{0}\to\infty$ , and has a minimum for intermediate values of $k_{0}$ .

In the limit of $k_{0}\gg k_{y}$ , the electrons are isotropised, so that the ten moment model approaches the five moment limit, while as $k_{0}\to 0$ , the pressure tensor is allowed to evolve freely. The calculation with nonlocal electrons has growth rate close to the case with $k_{0}=k_{y}$ , which is consistent with the $k_{y}$ dependence of the nonlocal heat flux.

\add

Which closure approximation best captures the LHDI remains a question for further study. As mentioned earlier, Equation (18) \addis the local dispersion relation in the cold electron limit, where finite Larmor radius effects can be neglected. For arbitrary $T_{e}$ , however, it is still necessary to incorporate these effects (Davidson \BOthers., \APACyear1977) \add. Outside the current sheet, electrons are strongly magnetised, so a gyrofluid closure with finite Larmor radius approximations would work ( $k_{y}=k_{\perp}$ for this geometry) (Snyder \BBA Hammett, \APACyear2001; Tassi \BOthers., \APACyear2018), \addwhile within the current sheet, where electrons are unmagnetised, the nonlocal closure discussed above would correctly capture the electron response. A transition between the two limits will thus be necessary to capture the instability correctly in the ten-moment model.

\remove

In order to properly describe the electron dynamics, a more complex closure with spatial dependence will be required as the electrons are magnetised outside the current sheet while being unmagnetised within the sheet.

5.2 Simulations of the LHDI

\note

changed to section 5.2

The five and ten moment equation systems have been implemented in the finite-volume version of the Gkeyll code, which uses a high-resolution wave propagation method for the hyperbolic part of the equations and a point implicit method for the source terms (Hakim \BOthers., \APACyear2006; Hakim, \APACyear2008), and has previously been used to study magnetic reconnection (Wang \BOthers., \APACyear2015; Ng \BOthers., \APACyear2015, \APACyear2017). The kinetic simulations use the discontinuous Galerkin \addfinite-element Vlasov-Maxwell solver of Gkeyll 2.0 (Juno \BOthers., \APACyear2018). \addBecause the Vlasov code uses a discontinuous Galerkin method, we require a basis function expansion in each cell, and we choose piecewise quadratic basis functions from the Serendipity Element family. Details on the particulars of the basis expansion can be found in (Arnold \BBA Awanou, \APACyear2011; Juno \BOthers., \APACyear2018).

The simulations presented below use the parameters $\rho_{i}=L$ , $v_{t,e}=0.06$ , $m_{i}/m_{e}=36,T_{i}/T_{e}=10$ , with simulation domain $L_{y}\times L_{z}=6.4L\times 12.8L$ . A background plasma with $n_{b}=0.001n_{0}$ is introduced for numerical stability. In the fluid simulations the grid size was $N_{x}\times N_{y}=256\times 512$ . The kinetic simulations are run in two velocity dimensions (2X2V) as the LHDI (with no guide field and $k=k_{y}$ ) does not depend on the out-of-plane velocity. The configuration space dimensions are the same as the fluid simulations, but use $96\times 192$ cells, with quadratic Serendipity elements (Arnold \BBA Awanou, \APACyear2011)\notechanged cite to citep in each cell. The electron velocity domain ranges from $-8v_{t,e}$ to $8v_{t,e}$ while the ion velocity domain ranges from $-6v_{t,i}$ to $6v_{t,i}$ in each direction. Two cells with quadratic serendipity elements are used per species thermal velocity. Quadratic Serendipity elements give us roughly a factor of 3 in additional sub-cell resolution for a total amount of resolution of 2.3 grid points per $\rho_{e}$ (calculated using the asymptotic $B$ field) and 6 cells per thermal velocity.

These parameters were chosen to balance computational costs while maintaining $T_{i}\gg T_{e}$ and the local approximation where gradients of the perturbed quantities in the $y$ direction are much larger than gradients in the $z$ direction. In the simulations, an initial $m=8$ perturbation is imposed, which corresponds to $k_{y}L\approx 7.9$ or $k_{y}\rho_{e}\approx 0.41$ and is close to the wavelength with the maximum growth rate predicted by local theory (Davidson \BOthers., \APACyear1977). For these parameters, the predicted kinetic growth rate is $\gamma=1.1\Omega_{ci0}$ and the most unstable region is at $z/L\approx 1.6$ . \addIn order to compare the (2X2V) simulations to the fluid simulations, we use an adiabatic index of 2 for the five-moment simulations, and modify the ten-moment model to relax only the in-plane components of the pressure tensor when using the local closure. We do not modify the nonlocal model as the closure does not couple the out-of-plane diagonal component of the pressure tensor to the in-plane components.

Fig. 11 shows a comparison between the structure of $E_{y}$ in kinetic and fluid simulations of the LHDI at $t\Omega_{ci0}=6$ . The local closure used $k_{0,i}=k_{y}$ and $k_{0,e}=0$ . The calculation with the nonlocal closure for ions also used $k_{0,e}=0$ . From the simulations, the measured growth rates of the $m=8$ mode were \change $0.39\Omega_{ci0}$$0.34\Omega_{ci0}$ for the local closure, $0.84\Omega_{ci0}$ for the nonlocal closure and $1.1\Omega_{ci0}$ for the kinetic model. The mode was found to be stable when using the five moment model. As can be seen in the lower panels, the structure of the LHDI in this regime is well described by the nonlocal model, with the growth rate about $24\%$ slower than in the kinetic simulation. When using the local model, in spite of the slower growth rate, the LHDI does eventually develop with a similar structure. \addWhen performing the equivalent fluid simulations with three velocity dimensions – with an adiabatic index of $5/3$ or relaxing all the components of the pressure tensor – which would be used in simulations of physical systems, the results are similar. The five-moment model is stable for these parameters, while the local relaxation shows a small increase in the growth rate to $0.39\Omega_{ci0}$ .

The role of ion kinetic effects is highlighted in Fig. 12, which shows a cut of the ion distribution function $f(v_{y},v_{z}=0)$ at the edge, where the perturbed electric field is confined to, and centre of the current sheet. These points are marked in Fig. 11. At the edge of the current sheet, the ion resonance can be seen in the upper panel, where the phase velocity from the theoretical solution and the initial drift velocity are marked. At the centre of the sheet, the distribution remains close to Maxwellian, consistent with the electrostatic LHDI being confined to the edge of the current sheet.

6 Conclusion

We have performed calculations of the drift-kink and lower hybrid drift instabilities for Harris sheets using the five and ten moment two-fluid models. For the drift-kink instability, the ten-moment model has growth rates and wavenumbers comparable to the results of Vlasov-Maxwell calculations, unlike the five-moment, or standard two-fluid, model, which has faster growing modes at larger mass ratios and wavenumbers. The growth rates are not sensitive to the relaxation parameter in the range $k_{0,s}\lesssim 1/d_{s}$ . Additionally, the sausage moment is damped by the ten-moment model, which is consistent with kinetic studies Pritchett \BBA Coroniti (\APACyear1996); Daughton (\APACyear1999\APACexlab\BCnt2). Although the kink mode has a lower growth rate at high mass ratio, this result does not preclude its excitation as a secondary instability (Lapenta \BBA Brackbill, \APACyear2002), or the growth of ion-ion kink instabilities in the ten moment model (Karimabadi, Daughton\BCBL \BOthers., \APACyear2003; Karimabadi, Pritchett\BCBL \BOthers., \APACyear2003), which could be the topic of future work.

The results are consistent the fluid work of (Pritchett \BOthers., \APACyear1996; Daughton, \APACyear1999\APACexlab\BCnt1) in the long wavelength regime, and the scaling with physical parameters such as temperature ratio and mass ratio are consistent with kinetic models (Daughton, \APACyear1999\APACexlab\BCnt2). There are some differences compared to the fluid results of (Yoon \BOthers., \APACyear2002) due to the treatment of the pressure term in the momentum equation. The kink modes are more sensitive to the ion model used, which may be useful if reduced electron models are used to save on computational costs.

\add

In global simulations, the importance of using a closure that captures the kink mode correctly has been seen in global simulations of Ganymede (Wang \BOthers., \APACyear2018; Ng, \APACyear2019). In Ng (\APACyear2019), \addit was shown that when using a five-moment model to study magnetosphere dynamics, a kink instability was excited in the magnetotail after the formation of the tail current sheet. While the current sheet did not disrupt in this case, the instability caused the formation of a large-scale corrugated structure. In the ten-moment simulations (Wang \BOthers., \APACyear2018; Ng, \APACyear2019), \addwhich are expected to show better agreement with kinetic models for this instability, the growth of the kink instability was much reduced.

The LHDI can be observed in both five and ten moment models with the appropriate choice of closure. Based on the results of kinetic theory (Davidson \BOthers., \APACyear1977; Hirose \BBA Alexeff, \APACyear1972), it is clear that the ions should be modeled using a nonlocal closure, or a relaxation with the $k_{0,i}\approx k_{y}$ so that the ion resonance can be captured. However, there is sensitivity to the electron model used, and the instabilities have the largest growth rates when $k_{0,e}\to 0$ or $\infty$ . For parameters used in reconnection studies ( $k_{0,s}=1/d_{s}$ ) (Wang \BOthers., \APACyear2015; Ng \BOthers., \APACyear2015), both electromangetic and electrostatic instabilities are damped, with the electromagnetic modes being stable in some regimes.

Finally, we have performed comparisons of fluid and Vlasov-Maxwell simulations to show that the LHDI can be observed in our fluid simulations. The distribution function information demonstrates the importance of the ion resonance and applicability of the local kinetic theory for thicker sheets, and illustrates the utility of the Vlasov-Maxwell code in analysing distributions due to the lack of particle noise as compared to particle-in-cell simulations.

In the context of global simulations, where we would want to model magnetic reconnection in addition to these instabilities, it may not be possible to capture the kink, LHDI and reconnection simultaneously in certain regimes. Due to computational constraints, it would be very computationally intensive to use the nonlocal closure for the ions in global studies, and setting $k_{0,i}=k_{y}$ to capture the LHDI would only be appropriate for a small range of wavenumbers, and likely excite the drift-kink unphysically. Additionally, the LHDI has a reduced growth rate when using closure parameters similar to those in reconnection studies, though this may not be an issue for the electrostatic LHDI in sufficiently thin sheets. A compromise may be the use of the five moment model for electrons, though this would miss the electron pressure tensor effects on reconnection. Even with compromises, it is evident that the ten moment model is a significant improvement over MHD models and captures some key physics features of fully kinetic simulations, which cannot be used for global space weather studies.

There is potential for further development of the closure with the use of temperature gradients, which provides some heat flux while remaining computationally tractable (Allmann-Rahn \BOthers., \APACyear2018)\notechanged cite to citep, but more study on how this model affects reconnection and instabilities is required. With the current ten-moment model, the drift-kink instability and reconnection can be studied simultaneously, which avoids the unphysical growth of the kink mode in five moment models that can disrupt the current sheet.

Acknowledgements.

J. Ng, A. Hakim and A. Bhattacharjee are supported by NSF Grant AGS-1338944 and DOE contract DE-AC02-09CH11466. The work of J. Juno was supported by a NASA Earth and Space Science Fellowship, Grant no. 80NSSC17K0428. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. Data for the figures are available online at (Ng \BOthers., \APACyear2018). \addWe are grateful to J. TenBarge for valuable discussions. The calculations in this work made use of the scipy library (Jones \BOthers., \APACyear2001–).

Appendix A Electrostatic response in the various plasma models

How well the fluid models approximate the ion response is determined by the term proportional to $R(\zeta_{i})=1+\zeta_{i}Z(\zeta_{i})$ in (18). This can be calculated by solving the 1-D dispersion relation in the fluid models and finding the perturbed density (Hammett \BBA Perkins, \APACyear1990)

[TABLE]

where $\phi_{1}$ is the perturbed electrostatic potential.

In the five moment model with no thermal conductivity or viscosity, the response can be written as

[TABLE]

where $\gamma$ is the adiabatic constant (we use $5/3$ in this paper).

The nonlocal ten-moment model has (Hammett \BBA Perkins, \APACyear1990)

[TABLE]

with $\chi_{1}=2/\sqrt{\pi}$ , and the relaxation to local isotropy has

[TABLE]

where $\alpha=k_{0}/k$ , with $k_{0}v_{t}$ being the relaxation rate. The dependence of the response function on the relaxation parameter is quite explicit and shows why the choice of $k_{0}$ is important the regime where the ion resonance is important for the LHDI. There is also a subtle difference between relaxing the pressure to local isotropy and relaxing temperature fluctuations to the equilibrium temperature $T_{0}$ , which is responsible for the $5/3$ in the first term of the denominator.

Bibliography42

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1Allmann-Rahn \B Others . ( \APA Cyear 2018) \APA Cinsertmetastar allmann:2018 {APA Crefauthors} Allmann-Rahn, F., Trost, T. \BCBL \BBA Grauer, R. \APA Cref Year Month Day 2018. \BBOQ \APA Crefatitle Temperature gradient driven heat flux closure in fluid simulations of collisionless reconnection Temperature gradient driven heat flux closure in fluid simulations of collisionless reconnection. \BBCQ \APA Cjournal Vol Num Pages Journal of Plasma Physics 843905840307. {APA Cref DOI} 10.101 · doi ↗
2Arnold \BBA Awanou ( \APA Cyear 2011) \APA Cinsertmetastar arnold:2011 {APA Crefauthors} Arnold, D \BPBI N. \BCBT \BBA Awanou, G. \APA Cref Year Month Day 2011. \BBOQ \APA Crefatitle The serendipity family of finite elements The serendipity family of finite elements. \BBCQ \APA Cjournal Vol Num Pages Foundations of Computational Mathematics 113337–344. \Print Back Refs \Current Bib
3Bale \B Others . ( \APA Cyear 2002) \APA Cinsertmetastar bale:2002 {APA Crefauthors} Bale, S \BPBI D., Mozer, F \BPBI S. \BCBL \BBA Phan, T. \APA Cref Year Month Day 2002. \BBOQ \APA Crefatitle Observation of lower hybrid drift instability in the diffusion region at a reconnecting magnetopause Observation of lower hybrid drift instability in the diffusion region at a reconnecting magnetopause. \BBCQ \APA Cjournal Vol Num Pages Geophysical Research Letters 292433-1–33-4. {APA Cref URL} · doi ↗
4Carter \B Others . ( \APA Cyear 2002) \APA Cinsertmetastar carter:2002 {APA Crefauthors} Carter, T \BPBI A., Yamada, M., Ji, H., Kulsrud, R \BPBI M. \BCBL \BBA Trintchouk, F. \APA Cref Year Month Day 2002. \BBOQ \APA Crefatitle Experimental study of lower-hybrid drift turbulence in a reconnecting current sheet Experimental study of lower-hybrid drift turbulence in a reconnecting current sheet. \BBCQ \APA Cjournal Vol Num Pages Physics of Plasmas 983272-3288. {APA Cref URL} http://dx.d · doi ↗
5Daughton ( \APA Cyear 1998) \APA Cinsertmetastar daughton:1998 {APA Crefauthors} Daughton, W. \APA Cref Year Month Day 1998. \BBOQ \APA Crefatitle Kinetic theory of the drift kink instability in a current sheet Kinetic theory of the drift kink instability in a current sheet. \BBCQ \APA Cjournal Vol Num Pages Journal of Geophysical Research: Space Physics 103A 1229429–29443. {APA Cref URL} http://dx.doi.org/10.1029/1998 JA 900028 {APA Cref DOI} 10.1029/1998 JA 900028 \Print Back Refs \C · doi ↗
6Daughton ( \APA Cyear 1999 \APA Cexlab \B Cnt 1) \APA Cinsertmetastar daughton:1999 kink {APA Crefauthors} Daughton, W. \APA Cref Year Month Day 1999 \B Cnt 1. \BBOQ \APA Crefatitle Two-fluid theory of the drift kink instability Two-fluid theory of the drift kink instability. \BBCQ \APA Cjournal Vol Num Pages Journal of Geophysical Research: Space Physics 104A 1228701–28707. {APA Cref URL} http://dx.doi.org/10.1029/1999 JA 900388 {APA Cref DOI} 10.1029/1999 JA 900388 \Print Back Refs \Curr · doi ↗
7Daughton ( \APA Cyear 1999 \APA Cexlab \B Cnt 2) \APA Cinsertmetastar daughton:1999 {APA Crefauthors} Daughton, W. \APA Cref Year Month Day 1999 \B Cnt 2. \BBOQ \APA Crefatitle The unstable eigenmodes of a neutral sheet The unstable eigenmodes of a neutral sheet. \BBCQ \APA Cjournal Vol Num Pages Physics of Plasmas 641329-1343. {APA Cref URL} http://dx.doi.org/10.1063/1.873374 {APA Cref DOI} 10.1063/1.873374 \Print Back Refs \Current Bib · doi ↗
8Daughton ( \APA Cyear 2003) \APA Cinsertmetastar daughton:2003 {APA Crefauthors} Daughton, W. \APA Cref Year Month Day 2003. \BBOQ \APA Crefatitle Electromagnetic properties of the lower-hybrid drift instability in a thin current sheet Electromagnetic properties of the lower-hybrid drift instability in a thin current sheet. \BBCQ \APA Cjournal Vol Num Pages Physics of Plasmas 1083103-3119. {APA Cref URL} http://dx.doi.org/10.1063/1.1594724 {APA Cref DOI} 10.1063/1.1594724 \Print · doi ↗