2. Physical models and numerics

A plasma is accurately described by distribution functions $f_s(\boldsymbol {x},\boldsymbol {v},t)$ for each species $s$. The number of particles located between $\boldsymbol {x}$ and $\boldsymbol {x} + \textrm {d}\boldsymbol {x}$ with velocities between $\boldsymbol {v}$ and $\boldsymbol {v} + \textrm {d}\boldsymbol {v}$ is then given by $f(\boldsymbol {x},\boldsymbol {v},t)\,\textrm {d}\boldsymbol {x}\,\textrm {d}\boldsymbol {v}$. For collisionless plasmas, the Vlasov equation determines the evolution of the distribution functions:

(2.1)\begin{equation} \frac{\partial f_{s}}{\partial t} +\boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\nabla} f_{s} + \frac{q_{s}}{m_{s}} (\boldsymbol{E} + \boldsymbol{v} \times \boldsymbol{B}) \boldsymbol{\cdot} \boldsymbol{\nabla}_{v} \,f_{s} = 0. \end{equation}

Various physical quantities can be obtained from the distribution function by taking moments. Multiplying by powers of $\boldsymbol {v}$ and integrating over velocity space leads to expressions for particle density $n_s(\boldsymbol {x}, t) = \int f_s(\boldsymbol {x}, \boldsymbol {v}, t) \,\text {d}\boldsymbol {v}$, mean velocity $\boldsymbol {u}_s(\boldsymbol {x},t) = ({1}/{n_s(\boldsymbol {x}, t)}) \int \boldsymbol {v} f_s(\boldsymbol {x}, \boldsymbol {v}, t) \,\text {d}\boldsymbol {v}$, pressure $\boldsymbol{P}_{\boldsymbol{s}} = m_{s} \int \boldsymbol {v}' \otimes \boldsymbol {v}' f_{s} \,\text {d}\boldsymbol {v}$ and heat flux $\boldsymbol{Q}_{\boldsymbol{s}} = m_{s} \int \boldsymbol {v}' \otimes \boldsymbol {v}' \otimes \boldsymbol {v}' f_{s} \,\text {d}\boldsymbol {v}$. Here, $\otimes$ is the tensor product and $\boldsymbol {v}' = \boldsymbol {v}- \boldsymbol {u}$.

The plasma quantities evolve according to the so called multifluid equations. These can be derived by taking moments of the Vlasov equation and making a physical assumption in order to close the resulting hierarchy. If moments up to the pressure tensor are considered, the 10-moment equations follow (so named because 10 equations need to be solved):

(2.2)\begin{gather} \frac{\partial n_{s}}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} (n_{s} \boldsymbol{u}_{s}) = 0 , \end{gather}

(2.3)\begin{gather}m_{s} \frac{\partial (n_{s} \boldsymbol{u}_{s}) }{\partial t} = n_{s} q_{s} (\boldsymbol{E} + \boldsymbol{u}_{s} \times \boldsymbol{B}) - \boldsymbol{\nabla} \boldsymbol{\cdot} \mathcal{P}_{s} , \end{gather}

(2.4)\begin{gather}\frac{\partial \mathcal{P}^{s}_{ij}}{\partial t} - q_{s} (n_{s} u^{s}_{[i} E_{j]} + \frac{1}{m_{s}} \epsilon_{[ikl} \mathcal{P}^{s}_{kj]} B_{l}) ={-} (\boldsymbol{\nabla} \boldsymbol{\cdot} \mathcal{Q}_{s})_{ij} . \end{gather}

Here $\mathcal {P}_{s} = m_{s} \int \boldsymbol {v} \otimes \boldsymbol {v} f_{s} \,\text {d}\boldsymbol {v}$ and $\mathcal {Q}_{s} = m_{s} \int \boldsymbol {v} \otimes \boldsymbol {v} \otimes \boldsymbol {v} f_{s} \,\text {d}\boldsymbol {v}$ are the second and third moment of the distribution function (multiplied by mass), $\epsilon _{ikl}$ is the Levi–Civita symbol and the square brackets denote the sum over as many permutations of indices as needed to make the tensors symmetric, for example

(2.5)\begin{equation} u_{[i} E_{j]} = u_i E_j + u_j E_i. \end{equation}

This set of equations needs a closure approximation for the divergence of the heat flux which is discussed in § 3. Heat flux $\boldsymbol{Q}$ is related to the distribution function's moments according to ${Q}_{ijk} = \mathcal {Q}_{ijk} - u_{[i} \mathcal {P}_{jk]} + 2 m n u_i u_j u_k$.

The 10-moment equations reduce to five-moment equations, if the temperature (and pressure, respectively) is isotropic. Then (2.3) becomes

(2.6)\begin{equation} m_{s} \frac{\partial (n_{s} \boldsymbol{u}_{s}) }{\partial t} = n_{s} q_{s} (\boldsymbol{E} + \boldsymbol{u}_{s} \times \boldsymbol{B}) - \frac{1}{3} \boldsymbol{\nabla} (2 \mathcal{E}_{s} - m_s n_s u_s^2) - \boldsymbol{\nabla} \boldsymbol{\cdot} (m_s n_s \boldsymbol{u}_s \otimes \boldsymbol{u}_s) , \end{equation}

with scalar energy density $\mathcal {E}_s = ({m_{s}}/{2}) \int v^2 f_{s} \,\text {d}\boldsymbol {v}$. Additionally, assuming zero heat flux in this isotropic limit, (2.4) becomes

(2.7)\begin{equation} \frac{\partial \mathcal{E}_{s}}{\partial t} + \frac{1}{3} \boldsymbol{\nabla} \boldsymbol{\cdot} (\boldsymbol{u}_s (5 \mathcal{E}_s - m_s n_s u_s^2)) - q_s n_s \boldsymbol{u}_s \boldsymbol{\cdot} \boldsymbol{E} = 0. \end{equation}

Physical units in the simulations are normalised as follows. Length is normalised over ion inertial length $d_{i,0}$ based on density $n_0$, velocity over ion Alfvén velocity $v_{A,0}$ based on the magnetic field $B_0$, time over the inverse of the ion cyclotron frequency $\varOmega _{i,0}^{-1}$ and mass over ion mass $m_i$. Further normalisations are vacuum permeability $\mu _0 = 1$ and Boltzmann constant $k_B = 1$.

The Vlasov simulations in this paper use the positive and flux-conservative (known as PFC) method by Filbet, Sonnendrcker & Bertrand (Reference Filbet, Sonnendrücker and Bertrand2001), which is a semi-Lagrangian method, combined with back-substitution (Schmitz & Grauer Reference Schmitz and Grauer2006) for the velocity updates. For the fluid equations we employ a centrally weighted essentially non-oscillating (known as CWENO) method (Kurganov & Levy Reference Kurganov and Levy2000) and the third-order Runge–Kutta scheme by Shu & Osher (Reference Shu and Osher1988). Maxwell's equations are solved by means of the finite-difference time-domain (FDTD) method.

All of the simulations were performed using the graphics processing unit (known as GPU) accelerated ‘muphy2’ multiphysics simulation code developed at the Institute for Theoretical Physics I, Ruhr University Bochum. The code is designed for running kinetic, fluid and hybrid schemes either individually or spatially coupled to each other.

3. Gradient-driven heat flux closure

A heat flux term that models kinetic dissipation mechanisms needs to introduce damping into the 10-moment equations. This can be achieved, for example, with an expression that takes the form of Fourier's law,

(3.1)\begin{equation} \boldsymbol{q} ={-} nD\boldsymbol{\nabla} T, \end{equation}

with heat flux vector $\boldsymbol {q}$ and thermal conductivity $D$. Heat flux approximations based on temperature gradients have been used for plasmas in the classic Braginskii equations already (Braginskii Reference Braginskii1965) and can also be the basis of Landau fluid closures. Hammett & Perkins (Reference Hammett and Perkins1990) and Hammett et al. (Reference Hammett, Dorland and Perkins1992) used phase mixing theory to obtain in 1-D Fourier space $\tilde {q} = - n_{0} \chi _{1} \sqrt {2} ({v_{t}}/{|k|}) \textrm {i} k \tilde {T}$ where tildes denote perturbed quantities and the conductivity $D = \chi _{1} \sqrt {2} ({v_{t}}/{|k|})$ is dependent on the wavenumber $k$ and thermal velocity $v_t$. This closure is suitable for approximating kinetic Landau damping.

Fourier's law is trivially generalised to three dimensions by replacing the heat flux vector with the third-order heat flux tensor ${Q}$ and the scalar temperature with the second-order temperature tensor ${T}$. Then the gradient of the temperature tensor yields a third-order tensor. The conductivity from the Landau fluid closure can also easily be transferred as it contains only the absolute value of the wavenumber. To avoid calculations in Fourier space, Sharma et al. (Reference Sharma, Hammett, Quataert and Stone2006) and Wang et al. (Reference Wang, Hakim, Bhattacharjee and Germaschewski2015) introduced a simplification where the continuous wavenumber field is replaced by a single typical wavenumber $\boldsymbol {k}_0$, and we follow this approach. Landau damping is the kinetic damping of plasma oscillations, thus $1/|\boldsymbol {k}_{s,0}| = d_{s,0}$ is a natural choice, $d_{s,0}$ being the species inertial length based on $n_0$. We choose an initial ansatz $\boldsymbol {q} = - D\boldsymbol {\nabla } P$ slightly different from Fourier's law which originates from comparisons with kinetic simulations. Since Landau damping is modelled, damping is supposed to act on fast changing gradients due to wave-like plasma activities and not on gradients from slowly changing background quantities. Therefore, we subtract an equilibrium pressure which is comparable to the consideration of perturbed quantities as in the Hammett–Perkins closure. For species $s$, the resulting heat flux then has the form

(3.2)\begin{equation} \boldsymbol{Q}_{\boldsymbol{s}} ={-}\frac{\chi}{k_{s,0}}v_{t,s}\boldsymbol{\nabla}(\boldsymbol{P}_{\boldsymbol{s}} - \boldsymbol{P}_{\boldsymbol{s,0}}) \end{equation}

with the dimensionless parameter $\chi$ and an equilibrium pressure $\boldsymbol{P}_{\boldsymbol{s,0}}$. A Maxwell distribution is an equilibrium distribution at which no Landau damping occurs and the related pressure would then be isotropic. This motivation is in accordance with studies where a closure term that drives pressure towards isotropy led to good results (Hesse et al. Reference Hesse, Winske and Kuznetsova1995; Johnson & Rossmanith Reference Johnson and Rossmanith2010; Wang et al. Reference Wang, Hakim, Bhattacharjee and Germaschewski2015). The closure produces wave damping, but due to the approximation of wavenumber and pressure perturbation it does not always give the correct Landau damping rate. It is not obvious which kinetic effects are captured and to what degree, for example Hesse et al. (Reference Hesse, Winske and Kuznetsova1995) included an isotropisation term to account for phase space instabilities.

We need an approximation of $\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol{Q}$ as a closure to the 10-moment equations so we take the divergence. The derivative of $v_{t,s}$ is assumed to be small and is neglected, and we obtain the Laplacian of the pressure tensor which is simply the Laplacian applied to each of its components. Then the final expression (hereafter called the gradient closure) is

(3.3)\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{Q}_{\boldsymbol{s}} ={-} \frac{\chi}{k_{s,0}}v_{t,s}\nabla^{2}(\boldsymbol{P}_{\boldsymbol{s}} - p_s\, {\boldsymbol{\mathsf{id}}}) \end{equation}

with $\chi = 3$ and $k_{s,0} = 1/d_{s,0}$. The isotropic pressure is given by $p_s = ({P}_{s,xx} + {P}_{s,yy} + {P}_{s,zz})/3$ and ${\boldsymbol{\mathsf{id}}}$ denotes the identity matrix. The thermal velocity is defined as $v_{t,s} = \sqrt {k_B T_s / m_s}$. We find that $\chi = 3$ is a reasonable value in a broad range of plasma configurations, but modifications to this dimensionless parameter can be made for further improving the agreement with kinetic calculations in specialised scenarios. Equation (3.3) reduces to the isotropisation closure in Wang et al. (Reference Wang, Hakim, Bhattacharjee and Germaschewski2015) with the approximation $\nabla ^2 \approx k_0^2$. The gradient closure was first used in Allmann-Rahn et al. (Reference Allmann-Rahn, Trost and Grauer2018).

For numerical stability and in order to enforce positive temperatures, we set a floor for the diagonal elements of the pressure tensor. The Laplacian is evaluated using second-order central finite differences. This explicit method may need smaller time steps than the solution of the fluid equations themselves. Therefore, the computation of the closure is subcycled with a boundary exchange of the pressure tensor in-between each subcycle. The number of subcycles necessary depends on the cell size and the plasma configuration. In the 3-D reconnection simulation in this paper, for example, eight subcycles were used.

4. 2-D reconnection

4.1. Harris sheet

Magnetic reconnection can develop from a perturbation of the Harris sheet (Harris Reference Harris1962) equilibrium. As initial conditions for the simulation we use parameters similar to the GEM reconnection set-up (Birn et al. Reference Birn, Drake, Shay, Rogers, Denton, Hesse, Kuznetsova, Ma, Bhattacharjee and Otto2001) but with a larger domain of size $L_{x} \times L_{y} = (32{\rm \pi} \times 8{\rm \pi} )d_{i,0} \approx (100 \times 25)d_{i,0}$, resolved by $2400 \times 600$ cells. Further differences are in the ion-to-electron mass ratio which is set to $m_i / m_e = 100$ and in the speed of light, $c = 30v_{A,0}$. The domain is periodic in the $x$-direction, has conducting walls for fields and reflecting walls for particles in the $y$-direction and is translationally symmetric in the $z$-direction. The initial configuration of the magnetic field is $B_{x}(y) = B_{0} \tanh (y/\lambda )$ and the particle density is $n(y) = n_{0} {\rm sech}^{2}(y/\lambda ) + n_{b}$ with $\lambda = 0.5d_{i,0}$ and background density $n_{b} = 0.2n_0$. Temperature is defined by $n_{0} k_{B} (T_{e}+T_{i}) = B_{0}^{2} / (2 \mu _{0})$, $T_{i}/T_{e} = 5$. A perturbation of the magnetic field is added to break the equilibrium which takes the form $\boldsymbol {B} = \hat {\boldsymbol {z}} \times \boldsymbol {\nabla } \psi$ where the perturbation in the magnetic flux is given by $\psi (x,y) = 0.1 \cos (2 {\rm \pi} x / L_{x}) \cos ({\rm \pi} y / L_{y}) B_0 d_{i,0}$. In the Harris equilibrium, the current resulting from the magnetic field configuration is distributed among electrons and ions according to $u_{z,i}/u_{z,e} = T_i/T_e$.

The 10-moment gradient simulation shown in figure 1 has an $x$-point profile comparable to existing PIC simulations with large system sizes. The electron outflow jet near the $x$-point is visible in figure 1(b) in dark blue. Its structure matches the kinetic results shown in figure 2 of Nakamura et al. (Reference Nakamura, Nakamura, Varsani, Genestreti, Baumjohann and Liu2018b), including the inflow and outflow regions along the separatrix boundary. Nakamura et al. (Reference Nakamura, Genestreti, Liu, Nakamura, Teh, Hasegawa, Daughton, Hesse, Torbert and Burch2018a) modelled a reconnection event detected by MMS where the outflow jet is similar. A difference between kinetic and 10-moment fluid models is that the onset of reconnection can take significantly longer in fluid simulations depending on the heat flux closure (Wang et al. Reference Wang, Hakim, Bhattacharjee and Germaschewski2015; Allmann-Rahn et al. Reference Allmann-Rahn, Trost and Grauer2018). This is also the case in the present simulation, and although onset duration is not a particularly important aspect in 2-D Harris sheet reconnection (as current sheet formation is ignored anyway), it is important for the interaction with current sheet instabilities in 3-D simulations.

Figure 1.

Current and velocity profiles in 2-D reconnection. Out-of-plane current $j_{z} / (n_0 v_{A,0})$ (a) and electron outflow velocity $u_{x,e} / v_{A,0}$ (b) at $t = 65{\varOmega _{i,0}^{-1}}$.

After reconnection, a secondary island evolves at the $x$-point location. If perturbations are applied to break symmetry, this island is ejected, but in the perfectly symmetric case shown here, the island stays in the middle of the domain and two new $x$-points are created. When the current sheet enlarges in the $x$-direction, at later times secondary islands also evolve in kinetic simulations with zero guide field (Karimabadi, Daughton & Scudder Reference Karimabadi, Daughton and Scudder2007; Klimas, Hesse & Zenitani Reference Klimas, Hesse and Zenitani2008), but more often when a guide field is present (Drake et al. Reference Drake, Swisdak, Schoeffler, Rogers and Kobayashi2006). Ion-scale islands, also referred to as flux transfer events, have been observed in the magnetosphere (Hasegawa et al. Reference Hasegawa, Kitamura, Saito, Nagai, Shinohara, Yokota, Pollock, Giles, Dorelli and Gershman2016; Hwang et al. Reference Hwang, Sibeck, Giles, Pollock, Gershman, Avanov, Paterson, Dorelli, Ergun and Russell2016). Nevertheless, the 10-moment gradient model seems to be more susceptible to formation of secondary islands.

4.2. Island coalescence

Reconnection studies are often initialised with a perturbed current sheet as described in the previous section. This approach results in reconnection rates and an $x$-point structure comparable to reconnection events in the magnetosphere. The procedure of reconnection onset, however, differs from that in nature since the formation of the current sheet is not considered. A model that can be applied to the magnetosphere should also be able to capture the physics of current sheet formation appropriately. This can be tested using the island coalescence reconnection set-up where kinetic simulations showed lower reconnection rates for larger system sizes (Stanier et al. Reference Stanier, Daughton, Chacón, Karimabadi, Ng, Huang, Hakim and Bhattacharjee2015), and the 10-moment fluid model with the gradient closure can reproduce the reconnection rate scaling (Allmann-Rahn et al. Reference Allmann-Rahn, Trost and Grauer2018). The correct time development of reconnection is particularly important for the interplay with current sheet instabilities such as the LHDI. In order to see if the kinetic results can be reproduced using the parameters chosen here for Harris sheet reconnection and the LHDI, we perform island coalescence simulations for different system sizes with $\chi = 3$ in (3.3).

The initial geometry is defined by density $n = n_{0} (1 - \epsilon ^{2}) / (\cosh (y/\lambda ) + \epsilon \cos (x/\lambda ))^{2} + n_{b}$ and magnetic potential $A_{z} = - \lambda B_{0} \ln ( \cosh (y/\lambda ) + \epsilon \cos (x/\lambda )$ with a system size of $L_{x} \times L_{y} = (2 {\rm \pi} \lambda \times 4 {\rm \pi} \lambda )\ d_{i,0}$. For a more detailed description of the set-up see Allmann-Rahn et al. (Reference Allmann-Rahn, Trost and Grauer2018) (also Ng et al. (Reference Ng, Huang, Hakim, Bhattacharjee, Stanier, Daughton, Wang and Germaschewski2015) and Stanier et al. (Reference Stanier, Daughton, Chacón, Karimabadi, Ng, Huang, Hakim and Bhattacharjee2015)). Reconnection rate is evaluated using the magnetic flux, which is the integral over $B_x$ from the $O$- to the $X$-point $\varPsi = \int ^{{X}}_{{O}}B_{x}\,{\textrm {d} y}$ and is normalised over $B'$, the maximum of the magnetic field's absolute value at $x=0$ and time $t = 0$. With $v'_{A} = B'/\sqrt {\mu _{0} n_{0} m_{i}}$, the normalised reconnection rate is $E_{R} = ({\partial \varPsi }/{\partial t}) / (B' v'_{A})$.

In figure 2 the average reconnection rate is shown depending on the parameter $\lambda$ which determines the system size. The dotted line represents the scaling of $\langle E_R \rangle$ with $\lambda$ found by Stanier et al. (Reference Stanier, Daughton, Chacón, Karimabadi, Ng, Huang, Hakim and Bhattacharjee2015) in kinetic PIC simulations up to $\lambda = 25d_{i,0}$, and the 10-moment gradient model comes close to $\langle E_R \rangle \propto (\lambda / d_{i,0})^{-0.8}$. Slightly better agreement can be achieved by tuning $\chi$ and $k_{s,0}$, but it is important that the closure works fine with the same parameters in different configurations – and that seems to be the case for $\chi = 3$ and $k_{s,0} = 1/d_{s,0}$. In addition to the 10-moment gradient results, we show performance of the five-moment model, which is able to capture the LHDI (cf. § 5) and was not considered in previous studies of island coalescence. The drop of reconnection rate is much smaller and the scaling with system size is only somewhat stronger than that of Hall MHD in Stanier et al. (Reference Stanier, Daughton, Chacón, Karimabadi, Ng, Huang, Hakim and Bhattacharjee2015). The 10-moment model using the isotropisation closure or a non-local closure outperforms the five-moment model, although they do not reach the scaling of the kinetic and 10-moment gradient models (Ng et al. Reference Ng, Hakim, Bhattacharjee, Stanier and Daughton2017; Allmann-Rahn et al. Reference Allmann-Rahn, Trost and Grauer2018).

Figure 2.

Scaling of average reconnection rate with system size in island coalescence reconnection: 10-moment gradient and five-moment simulations.

5. LHDI

When magnetic reconnection is simulated in two dimensions, the assumed approximation is that the plasma is homogeneous along the out-of-plane dimension. Thin current sheets, as they are present in reconnection, are often susceptible to instabilities though, and these can cause highly turbulent features in the out-of-plane dimension. The initial set-up for a 2-D analysis of various current sheet instabilities can again be a Harris sheet equilibrium. The connection to the 2-D reconnection set-up is that instead of the $x$–$y$-plane, now the $y$–$z$-plane is simulated for one fixed point in the $x$-direction. So the initial configuration is again $n(y) = n_{0} {\rm sech}^{2}(y/\lambda ) + n_{b}$, $B_{x}(y) = B_{0} \tanh (y/\lambda )$, $n_{0} k_{B} (T_{e}+T_{i}) = B_{0}^{2} / (2 \mu _{0})$ and $c = 30v_{A,0}$. For the velocity profiles it follows that $u_{z,s} = 2 T_{s,0} B_0 {\rm sech}^{2}(y/\lambda ) / (q_s n_s \lambda \mu _0)$. The instabilities are initiated by adding perturbations which can either be random noise or explicit modes that are supposed to be examined. In PIC simulations, the particle noise also excites instabilities, and when single precision calculations are used (especially combined with the FDTD method) even the noise from machine precision drives instabilities. Therefore we utilise double precision calculations in all cases in order to avoid unphysical results. There are four important aspects to consider when studying the LHDI: the growth rate; the duration of onset; the saturation; and the potential excitation of secondary instabilities. These characteristics are most influenced by current sheet thickness, background density and ion–electron temperature ratio. Generally, the growth rate is higher for thin current sheets, low background density and high $T_i/T_e$.

5.1. Electrostatic branch

The electrostatic LHDI is often the first current sheet instability to arise and generates fluctuations located at the edges of the current sheet. For studying the electrostatic LHDI branch we choose parameters as in Ng et al. (Reference Ng, Hakim, Juno and Bhattacharjee2019) so that there is the possibility of comparing the results. The parameters are: half-thickness of the current sheet $\lambda = \rho _{i,0} \approx 0.67 d_{i,0}$,; background density $n_b = 0.001n_0$; and temperature ratio $T_{i,0}/T_{e,0} = 10$, where $\rho _{i,0} = \sqrt {k_B T_{i,0} / m_i} / \varOmega _{i,0}$ is the ion gyroradius. These are suitable conditions for a fast development of the LHDI. Ng et al. (Reference Ng, Hakim, Juno and Bhattacharjee2019) used a mass ratio of $m_i/m_e = 36$ which is interesting for magnetosphere modelling where reduced mass ratios are common. Nevertheless, it is necessary to consider more realistic mass ratios, both for obtaining physically relevant results and for comparing numerical models. In order to analyse the influence of mass ratio, we perform simulations using different electron masses. Perturbations that exist in realistic plasmas can be modelled by adding random noise to the initial set-up. That way, the contribution and interaction of all modes as well as the occurrence of other instabilities than the LHDI can be observed. For the simulations in this subsection, linearly distributed random noise is added to the initial $B_{x}$ with a magnitude of $10^{-4} B_0$. The simulated domain is of size $L_y \times L_z = (12.8 \times 6.4) \rho _{i,0} \approx (8.6 \times 4.3) d_{i,0}$ with a spatial resolution of $512 \times 256$ for $m_i/m_e \leqslant 100$ and $1024 \times 512$ for higher mass ratios. In the Vlasov simulations, the extent of velocity space is ${\pm } 8 v_{t,s}$ resolved by $2.1$ cells per $v_{t,e}$. In the $m_i/m_e = 250$ case, for example, electron velocity space goes from $-27 v_{A,0}$ to $27 v_{A,0}$ and ion velocity space from $-5 v_{A,0}$ to $5 v_{A,0}$ with a total resolution of $1024\times 512\times 114^2$ cells.

Growth rate, onset duration and the mode number that develops are shown in figure 3 for the Vlasov model, the 10-moment model using the gradient closure, the five-moment model and the 10-moment model using the isotropisation closure in Wang et al. (Reference Wang, Hakim, Bhattacharjee and Germaschewski2015). Growth rate $\gamma$ is determined via a fit in the phase where $||E_z||_{\infty } \propto \exp (\gamma t)$ and we define the instability's onset duration as the time until $||E_z||_{\infty } \geqslant 0.02B_0 v_{A,0}$ which was the value where exponential growth started in the simulations.

Figure 3.

Scaling of LHDI characteristics with mass ratio for different models. (a) Growth rate of the instability, (b) onset duration and (c) number of modes.

The growth rate agrees reasonably well between the kinetic Vlasov model and the 10-moment gradient model. For the $m_i/m_e = 100$ case, time development of the electric field is shown in figure 4 where the slope of the curves correspond to the growth rates. The growth rates are $2.2\varOmega _{i,0}$ and $2.4\varOmega _{i,0}$ for the Vlasov and 10-moment gradient models, respectively, whereas growth rate is higher in the five-moment model ($3.2\varOmega _{i,0}$) and significantly lower in the 10-moment isotropisation model ($0.67\varOmega _{i,0}$). The latter approximates the divergence of heat flux as $\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol{Q}_{\boldsymbol{s}} = k_{s,0}v_{t,s}\ (\boldsymbol{P}_{\boldsymbol{s}} - p_s\,{\boldsymbol{\mathsf{id}}})$ and we choose $k_{s,0} = 1/d_{s,0}$ for the free parameter which was used in reconnection studies. When $k_{s,0} \rightarrow \infty$, the isotropisation closure approaches the five-moment limit (isotropy) and the growth rate increases accordingly. In this case, however, reconnection cannot be modelled well. Vlasov simulations with high mass ratios are not feasible, but we performed 10-moment gradient simulations with realistic mass ratio, and the growth rate goes up to $4\varOmega _{i,0}$. This should be considered when interpreting magnetosphere simulations where severely reduced mass ratios such as $m_i/m_e = 25$ are employed.

Figure 4.

Time development of the supremum norm of $E_{z}$ in LHDI simulations ($m_i/m_e = 100$).

The onset duration is, for Harris sheet reconnection simulations, as relevant as the growth rate since the influence of the LHDI at a certain time results from the combination of both. It is important to note that the onset duration depends heavily on the magnitude of the initial perturbation. But when the perturbation is fixed, the duration until exponential growth starts can be examined for the various models. In figure 5 the structure and magnitude of the electric field at a certain time is shown, which gives an impression of the combined effect of onset duration and growth rate. The difference between the 10-moment gradient model and the Vlasov model is small and the slightly higher $E$-field value in the Vlasov case results from the earlier onset. The five-moment model is already in a state where the fastest growing mode has saturated and secondary modes are excited. This rapid development can be attributed more to the very fast onset than to the higher growth rate. The five-moment results are in contrast to the findings of Ng et al. (Reference Ng, Hakim, Juno and Bhattacharjee2019) who reported a much too slow LHDI growth using this model. They agree, however, with a study by TenBarge et al. (Reference TenBarge, Ng, Juno, Wang, Hakim and Bhattacharjee2019) who modelled an MMS reconnection event and found the five-moment model to be not well-suited due to explosive growth of the LHDI. Our findings are also consistent with the $k_{s,0} \rightarrow \infty$ limit of the 10-moment isotropisation closure. The onset time in the 10-moment isotropisation model using normal parameters is very late, decreases with mass ratio, and for $m_i/m_e=250$ no instability grows at all during the simulated time span.

Figure 5.

Here $E_z/(B_0 v_{A,0})$ in simulations of the LHDI at $t = 4.75\varOmega _{i,0}^{-1}$ for the Vlasov model, the 10-moment gradient model and the five-moment model ($m_i/m_e = 100$).

The mode number shown in figure 3 is defined as the $M=k_z L_z/(2 {\rm \pi} )$ that develops out of random noise in the respective simulation and thus corresponds to the fastest growing mode (or one of the fastest, if growth rates are very close for several modes). The mode numbers do not agree between the Vlasov and fluid models for low mass ratios. While the number of modes itself is of little relevance for reconnection, it is crucial to always look at the respective fastest growing modes when comparing different models. Ng et al. (Reference Ng, Hakim, Juno and Bhattacharjee2019) explicitly initialised the $M=8$ mode in the $m_i/m_e=36$ case which is the fastest mode for the fluid models but a rather slow mode in the Vlasov model (we obtain a growth rate of $1.1\varOmega _{i,0}$ for the $M=8$ mode in agreement with Ng et al. (Reference Ng, Hakim, Juno and Bhattacharjee2019)). Therefore, only looking at a single mode is not an appropriate method of measuring the capability of fluid models to reproduce the kinetic LHDI. At mass ratio $m_i/m_e=250$ the Vlasov and the 10-moment gradient model results become closer and mode numbers keep increasing towards realistic mass ratio in the fluid simulations.

It is interesting to compare the simulation results with linear kinetic theory, especially at high mass ratios where the Vlasov simulations are not available for validating the fluid model. An expression for the local dispersion relation of the LHDI according to linear kinetic theory is given in (52) of Davidson et al. (Reference Davidson, Gladd, Wu and Huba1977). The dispersion relation can be solved numerically to give the dependence of growth rate on the wavenumber. The local density and magnetic field values were taken at the $y$-location of the fastest modes in the simulations. The results from linear theory are shown in figure 3 with dotted lines. The growth rate agrees well between theory and the nonlinear Vlasov simulation at $m_i/m_e = 36$ but the theory estimates slightly higher growth rates than the kinetic simulation at more realistic mass ratios. This suggests that nonlinear electron effects play an important role. In the theoretic prediction of the mode number the discrepancy is more evident; there is excellent agreement between kinetic theory and the Vlasov simulations at low mass ratios, but at $m_i/m_e = 250$ the theory does not reproduce the higher mode number from the simulation. When the Vlasov simulation results are not available, the gradient fluid model closely matches the growth rates predicted by linear kinetic theory. Concerning the mode number of the fastest growing mode, the theory yields better results at low mass ratios, but the fluid simulation can capture nonlinear effects and the interaction between various modes leading to more reliable results when realistic mass ratio is approached.

The mass ratio dependence of LHDI onset and growth rate simulations supports our choice of $k_{s,0} = 1/d_{s,0}$ in the gradient closure. More heat flux due to the closure leads to increased LHDI growth rates and earlier onset. Since $d_{s,0} \propto \sqrt {m_s}$, the higher $k_{s,0}$ compensates the increased heat flux caused by $v_{t,s}$ in the closure at small electron masses and prevents the electrons from being thermalised. If $k_{s,0}$ is chosen to be independent of mass, the LHDI grows too fast at high mass ratios, whereas with $k_{s,0} = 1/d_{s,0}$ there is good agreement between kinetic and 10-moment gradient results.

5.2. Electromagnetic branch

The electromagnetic LHDI has lower growth rates compared with the electrostatic branch, but it penetrates into the centre of the current sheet and thus can potentially have more significant impact on (for example) reconnection. Therefore, we test whether the gradient fluid model can also reproduce these more complex electromagnetic modes. A detailed analysis of the electromagnetic branch using kinetic theory and simulations can be found in Daughton (Reference Daughton2003). Since we are interested in the effect of the electromagnetic LHDI within reconnection, we choose parameters and resolutions as in the fluid simulation of 3-D reconnection that is discussed in the next section. The parameters are $n_b = 0.2n_0$, $T_i/T_e = 5$, $\lambda = 0.5d_{i,0}$, $m_i/m_e=100$ and initial random noise added to electron density with a magnitude of $\delta n_e = 10^{-8}n_0$. It should be noted that an electromagnetic LHDI branch is also present in Vlasov and fluid simulations when the set-up from the previous section is used. The domain is now of size $L_y \times L_z = (6 {\rm \pi} \times 2 {\rm \pi} ) d_{i,0}$ resolved by $450\times 150$ cells.

In figure 6 the fluctuations in the magnetic field caused by the electromagnetic LHDI are shown for kinetic Vlasov and 10-moment fluid simulations next to the current sheet density to demonstrate the influence of the mode. Onset duration of both the electrostatic and the electromagnetic LHDI modes is longer using the fluid model compared with the kinetic simulation with these plasma parameters. Although an accurate representation of onset times in all parameter regimes would of course be desirable, this behaviour is consistent with the longer reconnection onset of the fluid model that was mentioned in § 4. Growth rate of the electromagnetic modes that are shown in figure 6 (concerning the amplitude of magnetic field fluctuations) is $\gamma \approx 0.41\varOmega _{i,0}$ in the Vlasov simulation and $\gamma \approx 0.21\varOmega _{i,0}$ in the fluid simulation. In order to still compare the influence of the mode in both models, we choose times in the plot when the magnetic fluctuations have a similar magnitude. The kinetic Vlasov solution has the typical spatial structure of the electromagnetic mode at $t = 26\varOmega _{i,0}^{-1}$ which is not perfectly reproduced by the 10 moment gradient model ($t = 54\varOmega _{i,0}^{-1}$), but in the latter the mode does also arise in the centre of the current sheet and the wavelength is similar. At this stage the electromagnetic LHDI already causes significant kinking of the current sheet and the magnitude of the magnetic field fluctuations is high. In contrast to the findings in Daughton (Reference Daughton2003) the saturation is not at a moderate level. Instead, the fluctuations in the magnetic field keep growing so that $\tilde {B}_x > 1B_0$, causing significant modifications to the current sheet and magnetic field profile as can be seen in the lower panels of figure 6 (times $t = 28\varOmega _{i,0}^{-1}$ and $t = 57\varOmega _{i,0}^{-1}$, respectively). The reason for the high saturation level is unclear, but of course the results cannot be directly compared with Daughton (Reference Daughton2003) because the plasma parameters are different. In Harris sheet reconnection simulations, the reconnection process is started by an initial perturbation to the magnetic field. If strong fluctuations in the magnetic field due to the electromagnetic LHDI arise before reconnection, they can also initiate reconnection leading to multiple $x$-lines and turbulence. The high saturation level of the LHDI is surprising, considering that many studies suggest that very strong kinking of magnetosphere current sheets is primarily caused by slower instabilities such as the ion–ion kink instability or Kelvin–Helmholtz type instabilities (Lapenta & Brackbill Reference Lapenta and Brackbill2002; Lapenta et al. Reference Lapenta, Brackbill and Daughton2003; Baumjohann et al. Reference Baumjohann, Roux, Le Contel, Nakamura, Birn, Hoshino, Lui, Owen, Sauvaud and Vaivads2007). A detailed analysis of the saturation of the electromagnetic LHDI branch and its effect on current sheet stability using kinetic simulations in various parameter regimes goes beyond the scope of this paper but is a subject of future research.

Figure 6.

Electromagnetic mode of the LHDI in Vlasov and 10-moment gradient simulations. Shown are the deviation of the magnetic field from the initialised profile $\tilde {B}_x = B_x - \tanh (y/\lambda )$ and the out-of-plane current density $j_z$.

The fluid model cannot exactly represent the complex electromagnetic LHDI modes, but the key properties and the effect on the current sheet are in agreement with the Vlasov simulation including high saturation levels of the magnetic field fluctuations.

6. 3-D reconnection

In order to examine the effect of the LHDI on reconnection using the fluid model, we study the interplay of the two in a 3-D set-up that is identical to the one of the Harris sheet reconnection simulation in § 4.1 – except for the additional dimension and a tiny initial perturbation of electron density. The domain now has the size $L_{x} \times L_{y} \times L_{z} = (32{\rm \pi} \times 8{\rm \pi} \times 2{\rm \pi} )d_{i,0} \approx (100 \times 25 \times 6)d_{i,0}$ with a resolution of $2400 \times 600 \times 150$ cells and periodic boundary conditions in the $z$-direction. Linearly distributed random noise is added to the initial electron density with a maximum magnitude of $\delta n_e = 10^{-8} n_0$ so that instabilities like the LHDI can develop. The perturbation's magnitude influences the LHDI's onset duration and therefore has great impact on the effect of the LHDI, as will be discussed later. Using a small random perturbation of electron density, we give current sheet instabilities maximum freedom in their development.

Contour renderings of the current density $j_{z}$ are given in figure 7 after the reconnection current sheet has formed and at a later stage when turbulence dominates. Figure 7(a) shows the same time as figure 1 for the 2-D case and similar features are visible with a slightly smaller $x$-line current sheet in the 3-D simulation. Fluctuations in the third dimension are moderate in the electron diffusion region but strong in the outflow where the electromagnetic LHDI has impaired the current sheet. The electrostatic LHDI mode also affects the current density and can be seen in the contours at $t=65\varOmega _{i,0}^{-1}$ left and right of the $x$-line. Around $t=68\varOmega _{i,0}^{-1}$ a secondary island (plasmoid) forms and is ejected in the negative the $x$-direction. Later, the electron diffusion region becomes more turbulent and the $x$-line current sheet broadens, as is visible in figure 7(b).

Figure 7.

Contour plots of $j_{z} / (n_0 v_{A,0})$ from the 3-D reconnection simulation at times (a) $t = 65\varOmega _{i,0}^{-1}$ and (b) $t = 95\varOmega _{i,0}^{-1}$.

The reconnection rate is shown in figure 8(a) for both two and three dimensions, evaluated by averaging the value of $E_z$ at the $x$-line over $z$. We define the $x$-line as the location in the $x$–$y$ plane where $B_y$ averaged over $z$ goes through zero to the left (i.e. the negative $x$-direction) of the maximum of $B_y$. When two $x$-lines are present due to a magnetic island, this corresponds to the one at larger $x$-position as the island is ejected in the negative $x$-direction in the 3-D simulation. The island formation is visible in the form of a bump around $t=68\varOmega _{i,0}^{-1}$ and is associated with an increase of reconnection rate. The island is not ejected in the 2-D simulation with perfect symmetry which may cause the slightly higher reconnection rate at later times in the 2-D run. The earlier start of fast reconnection in three dimensions can likely be attributed to current sheet thinning and electron heating caused by the LHDI. After fast reconnection has started, the increased electron temperature and weaker density gradient in the electron diffusion region leads to a decay of the LHDI which may be the reason for the temporary drop of the reconnection rate at $t=48\varOmega _{i,0}^{-1}$. The development of instabilities at the $x$-line is demonstrated in figure 9 at three exemplary points in time. The LHDI (a) is weaker and arises later at the $x$-line than in the rest of the current sheet and when fast reconnection has started it decays and leaves moderate turbulence behind (b) caused by secondary modes. In the later stage, a different instability emerges (c) which has common features with the CSSI shown in figure 3 of Fujimoto & Sydora (Reference Fujimoto and Sydora2017) concerning wavelength, $y$-extent and magnitude. To confirm whether the instability in our simulation is indeed the CSSI will be left for future work as it is not the primary subject of this study. However, the structure of the electric field at the $x$-line does agree with the kinetic PIC simulation of 3-D reconnection in Fujimoto & Sydora (Reference Fujimoto and Sydora2017), although there the turbulence is much stronger around the $x$-line. In contrast, Nakamura et al. (Reference Nakamura, Genestreti, Liu, Nakamura, Teh, Hasegawa, Daughton, Hesse, Torbert and Burch2018a) found negligible three-dimensionality in the electron diffusion region. This is not necessarily a contradiction and is likely caused by the different background densities which were $0.044 n_0$ in the former and $0.3 n_0$ in the latter study, attended by fast or slow growth of the LHDI, respectively. The low temperature ratio of $T_i/T_e = 3$ in Nakamura et al. (Reference Nakamura, Genestreti, Liu, Nakamura, Teh, Hasegawa, Daughton, Hesse, Torbert and Burch2018a) may have further slowed down the LHDI development. Our set-up with $n_b=0.2 n_0$ and $T_i/T_e = 5$ lies between the two which is also evident in the amount of turbulence present. Recent studies on the role turbulence and the LHDI in reconnection have discussed a possible increase in reconnection rate due to modifications in Ohm's law caused by turbulence around the $x$-line (so called anomalous dissipation). The focus was on the electromagnetic LHDI in the reconnecting current sheet and there are findings that support the importance of turbulence (Price et al. Reference Price, Swisdak, Drake, Cassak, Dahlin and Ergun2016, Reference Price, Swisdak, Drake, Burch, Cassak and Ergun2017) and some that estimate a rather low influence of anomalous dissipation (Le et al. Reference Le, Daughton, Chen and Egedal2017, Reference Le, Daughton, Ohia, Chen, Liu, Wang, Nystrom and Bird2018). The small electromagnetic fluctuations at the $x$-line in our 3-D fluid simulation at earlier time (figure 9b) do have a similar structure as those shown in figure 7(b) of Price et al. (Reference Price, Swisdak, Drake, Burch, Cassak and Ergun2017), but in our case contributions of the electromagnetic LHDI at the $x$-line decay fast because of the low density gradient. Therefore, the LHDI cannot make long-term contributions apart from inducing the later CSSI-like instability. It should be noted that this could potentially be different in the reconnection scenario that aforementioned studies used where density and temperature asymmetries are present. The possible influence of the CSSI-like instability on reconnection rate will be an interesting topic of future research.

Figure 8.

(a) Reconnection rate, (b) supremum norm of $E_z$ in 3-D reconnection compared with 2-D LHDI and (c) LHDI onset duration depending on initial perturbation.

Figure 9.

Instabilities at the $x$-line in 3-D reconnection. Shown are $y$–$z$-slices of $E_z/(B_0 v_{A,0})$ at times (a) $t=42\varOmega _{i,0}^{-1}$, (b) $t=54\varOmega _{i,0}^{-1}$ and (c) $t=98.5\varOmega _{i,0}^{-1}$.

For agreement of the fluid simulations of 3-D reconnection with published kinetic PIC simulations, the relation between instability onset duration and reconnection onset duration is a central issue since both depend on the closure approximation. If reconnection onset takes too long with respect to LHDI onset, the electromagnetic LHDI can cause turbulent reconnection at multiple $x$-lines before reconnection due to the initially applied $B$-field perturbation takes place. The turbulence in the electron diffusion region is also largely determined by the time that the LHDI has been present before the reconnection process leads to a decay of the LHDI around the $x$-line. We compared the LHDI development in 2-D Vlasov and 10-moment gradient simulations with the one within 3-D reconnection using the same initial conditions for the current sheet. The electric field's evolution is given in figure 8(b). Comparing the 10-moment simulation of pure LHDI with the reconnection simulation, the development is very similar indicating that the LHDI is the main source of fluctuations in $E_z$, as expected. Its onset is slightly later within reconnection because of the broader current sheet that develops in the course of the reconnection process. Looking at the Vlasov simulation, the development is again comparable but onset time is much earlier and growth rate is also higher ($1.08\varOmega _{i,0}$ compared with $0.61\varOmega _{i,0}$ in the fluid model). The faster onset of the LHDI using the kinetic model in this set-up is consistent with the faster onset of reconnection so that in the end the effect of the LHDI in 10-moment gradient simulations resembles the effect in kinetic simulations. It should, however, be noted that the initially applied perturbation of $\delta n = 10^{-8} n_0$ is small compared with intrinsic PIC noise present in kinetic simulations of 3-D reconnection. A low initial noise level is necessary for agreement with previous PIC studies, which is interesting because reconnection onset and LHDI onset are delayed to similar extent compared with Vlasov simulations. The reason is that the electromagnetic branch of the LHDI has a high saturation level so that the magnetic field fluctuations can initiate turbulent reconnection at multiple $x$-lines.

The onset duration of the LHDI is important because it influences the the time when strong fluctuations in the magnetic field arise and largely determines the impact that the LHDI can have around the $x$-line. In order to examine the dependence of LHDI onset on the noise level, we apply initial perturbations of different magnitudes $\delta n$ to electron and ion density in the current sheet set-up used for the Harris sheet reconnection simulations, and measure the onset duration (time until $||E_z||_{\infty } \geqslant 0.04B_0 v_{A,0}$). The results are shown in figure 8(c). The onset duration has an approximately logarithmic dependence $t_{onset} \sim \log (\delta n / n_0)$ in this set-up using a resolution as in the 3-D simulation. There is a wide range between $t_{onset} = 30\varOmega _{i,0}^{-1}$ when $\delta n = 10^{-10}n_0$ and $t_{onset} = 4\varOmega _{i,0}^{-1}$ when $\delta n = 10^{-2}n_0$ and the time development is substantially different even when comparing, for example, an initialisation using $\delta n = 10^{-6}n_0$ and one using $\delta n = 10^{-4}n_0$. There is also competition between the LHDI onset and the reconnection onset, which is not significantly influenced by the level of initial random noise but can be controlled by the magnitude of perturbation in the magnetic flux $\psi$. Thus, depending on the initial perturbation one can get highly turbulent reconnection with $x$-lines due to current sheet instabilities or weak effect of instabilities and a reconnection that is essentially 2-D. In fluid simulations there is the additional difficulty that onset times are affected by the closure approximation. Past studies of 3-D reconnection that used kinetic PIC simulations typically relied on the discrete particle noise to trigger instabilities. When comparing such simulations it is important to consider that instability onset depends on the, respectively used, numbers of particles, resolutions and noise reduction techniques. Discrepancies in past publications concerning the influence of instabilities and turbulence may partly be attributed to varying onset times caused by different noise levels. Generally, a prediction of the influence of instabilities in magnetospheric reconnection events using simulations with a preformed current sheet is a difficult task because accurate information on perturbation levels at the beginning of the simulated time span would be needed. Three-dimensional Harris sheet simulations provide valuable insight into the possible scenarios concerning turbulence in magnetospheric reconnection though.