Kinetic simulations underestimate the effects of waves during magnetic reconnection

In collisionless plasma systems, such as those found in space and astrophysical environments, electron and ion distribution functions show quite some deviation from the Maxwell-Boltzmann distribution (e.g. Graham et al. (2021); Shi et al. (2022); Burch and Phan (2016); Raptis et al. (2022); Maksimovic et al. (2005); Schekochihin et al. (2005)). The use of the Vlasov equation, with three dimensions in velocity space, and three in position space, is necessary to understand the complex dynamics of such physical systems. Kinetic simulations that evolve the particle distributions or approximate the distributions using macroparticles (e.g. Bowers et al. (2008); Markidis et al. (2010); Fonseca et al. (2002); Fujimoto (2018); Juno et al. (2018)) are a powerful tool used to solve the Vlasov equation.

In spite of the advances in modern computing capabilities, fully kinetic simulations, which treat ions and electrons as particles, are still limited. Many physical systems require multiple spatial and temporal scales to be resolved, leading to compromises in the physical parameters used in simulations. Using Earth’s magnetosphere as an example, studies of localised magnetic reconnection regions involve dimensions of tens to hundreds of ion inertial lengths in two- or three-dimensions (e.g. Egedal et al. (2012); Le et al. (2017); Price et al. (2020); Lapenta et al. (2018) ), while studies of the collisionless shock can go from hundreds to up to over a thousand ion inertial lengths (e.g. Lu et al. (2024); Savoini and Lembege (2015); Bessho et al. (2022)). Other simulations studying kinetic phenomena in the magnetotail use an artificially smaller system, yet still require scales of tens to hundreds of ion inertial lengths Totorica and Bhattacharjee (2023); Liu et al. (2014). To simulate these scales while still resolving electron kinetic physics, a reduced electron-ion mass ratio is generally employed, reducing the separation of scales, and the ratio between electron plasma and cyclotron frequencies is also reduced, reducing the ratio between the speed of light and the electron Alfvén speed $c/V_{Ae}$.

Kinetic simulations have been successful in the study of magnetic reconnection, and the use of reduced parameters in such simulations is supported by the result that during collisionless reconnection at ion spatial scales, the reconnection rate is insensitive to the exact electron physics Comisso and Bhattacharjee (2016); Birn et al. (2001); Cassak et al. (2017). However, there are instances where the numerical parameters have affected the results qualitatively. It was shown that large mass ratios allowed the development of a new regime of reconnection with embedded exhaust current layers Le et al. (2013), while Ref. Jara-Almonte et al. (2014) showed that increasing the frequency ratio leads to Debye scale turbulence developing in the reconnecting layer. Artificially low mass ratios also allow the development of drift-kink instabilities that disrupt current sheets Daughton (1998). The use of reduced parameters also significantly affects results in other fields such as shock physics. It is known that observations of electric field fluctuations in shock crossings are much stronger than those found in simulations because of the use of reduced parameters Wilson III et al. (2021); Bohdan et al. (2024).

In this work we use the specific example of lower-hybrid drift waves during magnetic reconnection to quantitatively study the effects of numerical parameters and their importance when coupling electron-scale waves to the ion-scale reconnection region and beyond. Lower-hybrid drift waves are driven by the diamagnetic drift Davidson et al. (1977); Huba and Wu (1976); Krall and Liewer (1971); Davidson and Gladd (1975) and are often found in magnetic reconnection regions. Their importance in the reconnection region is an area of active research Graham et al. (2022); Mozer et al. (2011); Yoo et al. (2024); Chen et al. (2020), where they can contribute to electron momentum balance through correlated density and electric field fluctuations known as anomalous drag, or cause electron heating.

There are numerous simulation, experimental and observational studies of the role of the lower-hybrid drift instability (LHDI) during magnetic reconnection, in which it is shown that they cause particle transport and mixing (e.g. Le et al. (2017, 2018); Price et al. (2020)). Although it is known that generating the waves requires proper scale separation between electrons and ions Daughton (2003); Roytershteyn et al. (2012); Price et al. (2016), existing studies do not discuss saturation and how numerical parameters affect calculations of the waves’ contribution to momentum balance in the larger reconnection region Le et al. (2017, 2018); Price et al. (2016); Roytershteyn et al. (2012). Theoretical studies of the waves themselves Winske and Liewer (1978) consider their saturation, but do not account for their coupling to magnetic reconnection. These limitations mean that quantitative comparisons between simulations and observations or experiments are not adequate when considering the broader effects of the waves.

Recently, it has been shown in experiments and observations that during guide-field reconnection (where there is a magnetic field parallel to the initial current sheet), the LHDI is excited by the electron flow in the reconnection exhaust Chen et al. (2020); Yoo et al. (2024). While waves are seen in simulations of these events, there are still discrepancies in both particle acceleration and electron momentum balance related to the wave amplitudes Ng et al. (2023, 2020). We perform simulations using the guide-field reconnection configuration to study the effects of numerical parameters on wave amplitudes and momentum balance. Due to computational limitations, we then use an alternative setup to study the waves driven by the LHDI in isolation and evaluate their scaling with numerical parameters systematically. Our results show that the normalized wave amplitudes and anomalous drag increase as parameters become more realistic.

We first give a brief introduction to the theory of LHDI saturation and how it relates to the reconnection electric field. The LHDI is driven by a diamagnetic current drifting across the magnetic field Davidson et al. (1977); Huba and Wu (1976); Krall and Liewer (1971). The simplest estimate of the saturation amplitude of the LHDI is given in Winske and Liewer (1978) as

where $\mathcal{E}$ is the electric field energy density, $n$ is the local density, $V_{d}$ is the relative velocity between electrons and ions, $m_{e}$ is the electron mass, and $\omega_{pe}/\omega_{ce}$ is the ratio between electron plasma and cyclotron frequencies. This equates the free energy to the wave and particle energy in the electrostatic limit. Other saturation mechanisms such as ion trapping exist and will be discussed later Brackbill et al. (1984); Winske and Liewer (1978); Lavorenti et al. (2021).

To compare the wave fields to magnetic reconnection, we note that in the collisionless limit, the reconnection electric field is generally around $0.1B_{0}v_{A0}$, where $B_{0}$ is the upstream magnetic field and $v_{A0}$ is the ion Alfvén speed calculated using upstream quantities Comisso and Bhattacharjee (2016); Birn et al. (2001); Cassak et al. (2017). In the guide-field configuration we study, $V_{d}$ is given by the difference between the electron and ion outflow speeds, and is proportional to the electron Alfvén speed $v_{Ae0}=B_{0}/\sqrt{\mu_{0}n_{0}m_{e}}$, where $n_{0}$ is the upstream electron density Karimabadi et al. (2007). This can be substituted into the expression for the energy density, and after some manipulation, we find

The first term $\omega_{ce}/\omega_{ce0}$ is the ratio of the local to upstream electron cyclotron frequency and can be treated as a scaling factor in this study, and other quantities on the right-hand side are measured locally. This expression explicitly scales with the mass ratio and the frequency ratio while the normalized reconnection electric field is not sensitive to the mass and frequency ratios.

We perform two types of simulations to study the variation of the lower-hybrid waves with numerical parameters – “reconnection” simulations and “wave” simulations. The two-dimensional reconnection simulations use the same physical parameters as Ng et al. (2023) where the initial conditions consist of a Harris sheet superposed on an asymmetric background, while wave simulations capture the growth of lower-hybrid waves in a current layer. The reconnection simulations illustrate the coupling of lower-hybrid waves to the reconnection process, while the wave simulations allow a study using more realistic parameters.

We first study the development of lower-hybrid waves during asymmetric guide-field reconnection, the initial conditions for the base case are exactly the same as Ng et al. (2023), which was based on MRX experiments. The asymmetric layer is defined by $B_{L}/B_{H}=1.25$, $T_{eL}=T_{eH}=T_{iH}$, $n_{L}/n_{H}=0.5$ and $T_{iL}/T_{eH}=1.23$ where the subscripts represent the asymptotic values on either side of the current sheet. Were this magnetopause reconnection, $L$ would correspond to the magnetosphere and $H$ the magnetosheath. A constant guide field $B_{g}=1.8B_{H}$ is present, and the electron beta on the high density side is $\beta_{eH}=2\mu_{0}n_{H}T_{eH}/(B_{H}^{2}+B_{g}^{2})=0.3$. Instead of performing 3D simulations, we rotate the $x$-$z$ plane in our two-dimensional simulations by $5.7^{\circ}$ so that the magnetic field is almost perpendicular to the simulation plane in the unstable region found in Ng et al. (2023), which is favorable for the excitation of lower-hybrid waves.

We show the results from two simulations with $m_{i}/m_{e}=25$ and $m_{i}/m_{e}=400$, both having $\omega_{peH}/\omega_{ceH}=2$ (the realistic value is approximately 125 in the MRX experiment). Both simulations have dimensions $L_{x}\times L_{z}=40d_{i}\times 20d_{i}$, using $1536\times 768$ and $12288\times 6144$ grid cells respectively. The asymmetric conditions cause a density gradient across the outflow region, and the initial configuration has the plasma beta lowest in the upper-right quadrant. Lower-hybrid drift waves develop, similar to the original 3D simulation with $m_{i}/m_{e}=100$ Ng et al. (2023). The fluctuating electric fields of the lower-hybrid waves are illustrated in Figure 1. The normalized wave amplitudes are clearly higher in the mass ratio 400 simulation. The main reason for the larger normalized amplitude is because the normalized outflow $u_{ex}/v_{A0}$ increases as mass ratio increases as implied by Equation (2).

Although these results suggest that the normalized wave electric field increases as realistic parameters are approached, it is challenging to perform a systematic scan because these waves only appear after reconnection has developed, and there are differences in evolution between simulations (such as plasmoid formation that changes the overall structure of the reconnection region). As such, we perform a set of simulations focused only on the generation and evolution of the waves.

The initial setup for a wave simulation is a Vlasov confinement equilibrium with a narrow current layer in the center of the $x$ domain Hewett et al. (1976). This has previously been used for studies of the contribution of the LHDI to anomalous resistivity Winske and Liewer (1978). The ion density is given by

where

The ion temperature is uniform and isotropic, and the electron temperature components $T_{xx}$ and $T_{zz}$ are uniform, while the initial magnetic field, electric field and other electron quantities are determined self-consistently using the recursive relations derived in Ref. Hewett et al. (1976). We set $\alpha=5$, $a=20$ and $n_{\infty}=0.4$, where $\alpha$ and $a$ control the thickness of the current layer and $n_{\infty}$ is the asymptotic density on the positive $x$ side of the simulation. The electron distributions are expressed as a sum of Hermite polynomials. We cut off the infinite sum at 4 as adding higher-order terms does not change the distribution significantly. Similar configurations have been used to study the LHDI in other contexts Lavorenti et al. (2021). Note that the recursive relations here do not depend on the mass ratio.

The computational domain is $L_{x}\times L_{y}=40d_{e}\times 50d_{e}$ covered by $320\times 450$ cells, and the plasma parameters for the baseline case defined at $x=0$ are $\omega_{pe0}/\omega_{ce0}=1$, $T_{i0}/T_{e0}=1.25$ and $\beta_{e}=0.16$. The ion-to-electron mass ratio is $m_{i}/m_{e}=100$ and there are 400 particles per species per cell. These parameters are chosen so the plasma in the current layer is similar to the unstable region in Ng et al. (2023). These simulations are much smaller than the reconnection simulations, allowing us to perform a wider parameter scan by varying $m_{i}/m_{e}$ from 100 to 1836, and $\omega_{pe}/\omega_{ce}$ from $0.5$ to $8$. The resolution is increased appropriately as simulation parameters change so that the electron Debye length is resolved. The number of particles per cell is increased to 3200 for the $\omega_{pe0}/\omega_{ce0}=8$ simulations. The initial conditions are shown in Figure 2 for $m_{i}/m_{e}=1836$, $\omega_{pe0}/\omega_{ce0}=2$.

The evolution of the system is shown in Figure 3 for the simulation with $m_{i}/m_{e}=1836$ and $\omega_{pe0}/\omega_{ce0}=2$. In this case, the lower-hybrid drift waves develop in less than one ion-cyclotron time, shown by the $E_{y}$ signatures in the center of the domain. The peak fluctuation amplitudes are measured as the current layer starts to break up, corresponding to $t\Omega_{ci}=0.375$ in the Figure. This behavior is similar in all simulations, though the instability growth rate (is proportional to $\omega_{LH}$ Davidson et al. (1977)) relative to the ion cyclotron frequency increases with the mass ratio as $\omega_{LH}/\Omega_{ci}\propto\sqrt{m_{i}/m_{e}}$.

For each simulation, we calculate the root mean square of the fluctuating electric field $\delta E_{y,rms}$, and the anomalous drag term $\langle\delta n_{e}\delta E_{y}\rangle/\langle n_{e}\rangle$ by integrating along the $y$ direction. The maximum values are then evaluated along $x$. For comparison between different simulations, these quantities are normalized by $B_{0}v_{A0}$ evaluated at $x=0,t=0$ in each simulation.

The maximum values of the normalized $\delta E_{y,rms}$ during the simulation are displayed in Figure 4. For the electric field fluctuations, there is a clear increase in the normalized values with mass ratio, in addition to a weaker variation with $\omega_{pe}/\omega_{ce}$.

Fitting the variation with mass ratio gives a scaling of $(m_{i}/m_{e})^{0.56}$ and a correlation coefficient $r=0.98$, showing good agreement with the $\sqrt{m_{i}/m_{e}}$ expression in Equation (2), indicating the importance of the mass ratio to the normalized wave amplitudes in simulations.

The variation of the electric field amplitude does not show exact agreement with the $\omega_{pe}/\omega_{ce}/\sqrt{1+(\omega_{pe}/\omega_{ce})^{2}}$ relationship, but the general trend of a steeper increase at small $\omega_{pe}/\omega_{ce}$ and smaller increase at larger, more realistic conditions holds. Equation (2) does not include electron temperature or electromagnetic effects, which partially account for the larger discrepancies at $\omega_{pe}/\omega_{ce}\leq 1$.

The largest initial ratio between the drift velocity and ion thermal speed $V_{d}/v_{thi}\approx 3.5$ is found in the $m_{i}/m_{e}=1836$ simulations, which is within the empirical limits suggested by Ref. Winske and Liewer (1978) for the validity of Equation (1), and close to the $V_{d}/v_{thi}\approx 3$ boundary between Equation (1) and the ion trapping saturation mechanism found by Brackbill et al. (1984). While we see accelerated ions (not shown), there is surprisingly good agreement with Equation (1) at the larger $\omega_{pe}/\omega_{ce}$ relevant to laboratory and magnetospheric conditions. The electron resonance is unlikely to be important due to the low beta Huba et al. (1978).

The variation of the normalized anomalous drag term $\langle\delta n_{e}\delta E_{y}\rangle/\langle n_{e}\rangle$ is shown in Table 1. Similar to $\delta E_{rms}$, there is an increase of this term with the mass ratio. However, the scaling is stronger, with an order of magnitude increase going from $m_{i}/m_{e}=100$ to $1836$. There is no clear trend in the variation with $\omega_{pe}/\omega_{ce}$, with an increase with $\omega_{pe}/\omega_{ce}$ at mass ratio $400$ and a decrease at realistic mass ratio. While the electric field fluctuations scale like $\sqrt{m_{i}/m_{e}}$ as shown earlier, the density fluctuations in the simulations do not change with the mass ratio.

The additional increase of the drag term with mass ratio can be explained by the phase difference between $\delta n_{e}$ and $\delta E_{y}$ Yoo et al. (2024). We examine the quasi-linear expression for the anomalous collision frequency from Davidson and Gladd (1975),

where $\zeta_{i}=\omega/(kv_{thi})$, $u_{ey}$ is the electron velocity driving the instability, $Z$ is the plasma dispersion function and $k_{M}$ is the wavenumber at maximum growth. Assuming $k_{M}\rho_{e}\sim 1$, after some manipulation and normalization, the anomalous drag term is proportional to $(\omega_{pe}/\omega_{ce})^{2}/(1+(\omega_{pe}/\omega_{ce})^{2})\sqrt{m_{i}/m_{e}}\text{Im}(\zeta_{i}Z(\zeta_{i}))$. We find that $\zeta_{i}<1$ for all but the $\omega_{pe}/\omega_{ce}=8$, $m_{i}/m_{e}=1836$ case, and $\zeta_{i}$ increases with mass ratio, such that $\text{Im}(\zeta_{i}Z(\zeta_{i}))$ increases for the studied parameters (the maximum $\text{Im}(\zeta_{i}Z(\zeta_{i}))$ is at $\zeta_{i}\approx 0.7$).

Although this work focuses on the amplitudes of electrostatic lower-hybrid drift waves in guide-field reconnection, there are wider implications, both for lower-hybrid waves and kinetic simulations in general. Generalizing the setup to 3D allows the excitation of lower-hybrid drift waves with a small $k_{\parallel}$ with parallel electric fields that accelerate electrons Lavorenti et al. (2021); Ng et al. (2023); Marshall et al. (2022); Cairns and McMillan (2005). Variations in wave amplitudes would naturally affect the expected electron energies. We also briefly address the cancellation of the drag term due to anomalous viscosity Graham et al. (2022); Price et al. (2020). In our simulations, the ratio between the magnitude of the drag and viscous terms is on average $\approx 1.3$, such that the cancellation is not perfect. It is likely that this result is due to the shorter wavelengths seen in our simulations ($k\rho_{e}\approx 1$) compared to the magnetopause observations ($k\rho_{e}\approx 0.4$) Graham et al. (2022), meaning that electrons are not as strongly magnetized. Magnetotail observations show agyrotropic distributions where electron gyroradii have similar scales to the wavelength Chen et al. (2020), suggesting that the waves can contribute to momentum balance there.

In antiparallel or weak guide field reconnection simulations, electrostatic lower-hybrid waves are found outside the current sheet Daughton (2003); Mozer et al. (2011), and do not contribute significantly to the reconnection electric field at the x-line. Nonetheless, the conclusions of this work apply to the wave amplitudes. Eigenmode studies Daughton (2003) begin with a Harris equilibrium Harris (1962), where the relative velocity between electrons and ions $V_{d}\propto T/(eB_{0}L)$ where $e$ is the unit charge and $L$ is the width of the current sheet. One may perform a similar analysis for the electric field and find that $\delta E/(B_{0}v_{A0})\sim((\beta_{i}+\beta_{e})/L)(c/\omega_{LH})/\sqrt{4(1+\omega_{pe}^{2}/\omega_{ce}^{2})}$. If $L$ is $d_{i}$ scale, the mass ratio does not affect the normalized electric field, but if $L$ is $d_{e}$ scale, the factor of $\sqrt{m_{i}/m_{e}}$ reappears. In an evolving system, the current layer thins to sub-ion scales during reconnection, and the electron drift velocity increases, suggesting that the relative effects of the wave electric field would increase when using realistic parameters. Systems with electron-scale reconnection also see stronger flows, though the generation of the LHDI may be limited by their spatial extent. While waves at the separatrix or outside the current sheet may not be important for reconnection, observations show that they still affect the electron temperature and mixing Wang et al. (2022); Graham et al. (2022); Le et al. (2018).

To summarize, we have shown that even though the reconnection rate in ion-scale regions is not sensitive to the reduced parameters used in kinetic simulation, the global effects of lower-hybrid drift waves are underestimated, both in terms of their amplitudes relative to the reconnection electric field and the contribution to momentum balance through correlated density and electric field fluctuations. The results of this paper are applicable to simulations of other systems where electron-scale waves and instabilities couple to larger scales. For example, in collisionless shock simulations, fluctuating electric fields are much smaller than the quasi-static electric fields, in contrast to observations where the opposite is true Wilson III et al. (2021); Bohdan et al. (2024). Our results suggest that a careful analysis of how waves develop and saturate in multi-scale simulations is necessary, as current simulations could be underestimating their relative effects.

This work was supported by the MMS Mission, NASA Grants 80NSSC21K1462, 80NSSC21K1795, NNH20ZDA001N, and 80NSSC24K0094, and NSF Grant AGS2010231. Simulations were performed using NASA HECC resources.