The magnetic reconnection experiment has recently seen short wavelength ( k ρ e ∼ 1) lower-hybrid waves near the electron diffusion region in strong guide field reconnection. Based on plasma parameters from the experiment, we perform a three-dimensional fully kinetic simulation in order to investigate the generation of the lower-hybrid waves and their effects on the reconnection process. We find that the low-beta regions around the reconnection site are unstable to the lower-hybrid drift instability propagating in the outflow direction, driven by the difference between the electron and ion outflows. The waves modify the electron distributions, leading to periodic opening and closing of gaps in electron velocity space, and provide a small contribution to the anomalous resistivity. Finally, the simulation results are discussed in the context of space observations and laboratory experiments.

## I. INTRODUCTION

The lower-hybrid drift instability (LHDI) is a plasma instability driven by a diamagnetic current drifting across the magnetic field.^{1–3} During magnetic reconnection, the LHDI has been invoked as a source of the numerous lower-hybrid waves found around the reconnection region in both observations^{4–8} and laboratory experiments.^{9–12} In reconnection regions, the expectation is that short wavelength (electron-gyroradius scale) lower-hybrid waves with a primarily electrostatic character are found at the edges of current sheets, while longer wavelength electromagnetic waves can be found in the current sheet.^{13} However, recent work has shown that when there is a guide field, the shorter wavelength modes are seen inside the current sheet as well.^{6}

The role of these waves in reconnection has been studied in both simulations and observations. At both the magnetopause and in the magnetotail, simulations and observations have shown that the LHDI can cause the reconnection region to become turbulent, leading to enhanced mixing and transport of particles.^{14–17} In regions where the amplitude of the lower-hybrid waves is large, non-gyrotropic perpendicular electron acceleration is seen.^{6} Lower-hybrid waves are also seen in the reconnection outflow regions,^{18–20} which can lead to parallel electron heating.^{21}

Whether the waves contribute to the electron momentum balance and the reconnection electric field is a continuing area of research. It has been suggested that the lower-hybrid waves contribute to “anomalous” resistivity,^{22} but simulations with magnetopause-like parameters have shown that the contributions to the reconnection electric field are small.^{14,15} Observations at the magnetopause slightly away from the electron diffusion region have also shown that the anomalous resistive and viscous terms due to the waves can cancel, leading to a small contribution to the electric field, though particle transport is still present.^{17}

In laboratory experiments, however, recent measurements at the Magnetic Reconnection Experiment (MRX)^{9} have shown that under strong guide fields, the plasma conditions near the electron diffusion region are favorable for the excitation of lower-hybrid waves. Motivated by these experiments, we perform a three-dimensional particle-in-cell simulation of asymmetric reconnection in order to understand how and where the lower-hybrid waves are generated in the system. We show that in the low-beta outflow regions, the plasma is unstable to the LHDI, leading to the generation of lower-hybrid waves. While the effect of the waves on momentum balance is small, we show that the waves cause modification of the electron distribution function, leading to periodic opening and closing of gaps in the electron distribution in the parallel velocity direction.

## II. SIMULATION

^{23,24}The initial setup is an asymmetric current sheet, with the plasma components having Harris-sheet populations in addition to an asymmetric background. The reconnecting magnetic field and temperature are determined by

^{14}

*Q*

_{0},

*Q*

_{1}are the asymptotic upstream values. To satisfy force balance, the initial density is given by

*y*-direction $ B g = 1.8 B 0$ is present, and the electron beta on the high-density side is $ \beta e 0 = 2 \mu 0 n 0 T e 0 / ( B 0 2 + B g 2 ) = 0.3$. The mass ratio is $ m i / m e = 100$, and the ratio of the electron plasma frequency to the cyclotron frequency is $ \omega p e / \omega c e = 2$. A reduced mass ratio (compared to the proton mass) is chosen to provide separation between electron and ion scales while respecting computational limitations with a three-dimensional domain. Other possible effects of the mass ratio are discussed in Sec. IV.

The initial width of the current sheet is *L* = *d _{i}*, where

*d*is the ion inertial length on the high-density side. The simulation domain is $ L x \xd7 L y \xd7 L z = 40 d i \xd7 10 d i \xd7 20 d i$ covered by $ 3072 \xd7 768 \xd7 1536$ cells, initialized with 150 particles per species per cell. The simulation uses periodic boundary conditions in the

_{i}*x*and

*y*directions and conducting walls in the

*z*direction. Reconnection is initiated by a small perturbation, which causes an x-line to develop in the center of the domain. In this simulation, the high density plasma is initially in the

*z*< 100 region. Unless otherwise specified, we use simulation units with $ d e = 1 , c = 1$, and $ \omega p e = 1$. Simulation results are shown at $ t \omega c i = 49$, when reconnection is steady-state. Regarding the choice of parameters, the density ratio $ n 1 / n 0$, magnetic field ratio $ B 1 / B 0$, and guide field

*B*were chosen according to experimental parameters, while $ \omega p e / \omega c e$,

_{g}*L*=

*d*were chosen due to computational limitations.

_{i}## III. RESULTS

An overview of the system in a selected *x*–*z* plane is shown in Fig. 1. At this time, there is a well-developed current sheet with a single x-line as seen in the top panel. The strong guide field stabilizes lower-hybrid instabilities near the x-line,^{25,26} causing less visually pronounced structures in the *y* direction than those seen in other asymmetric simulations (e.g., Ref. 14).

The lower-hybrid drift instability that we will discuss in this paper can be seen in the central panel, which shows alternating *u _{ez}* signatures in the upper-right exhaust. The structures resemble those seen in Ref. 27, which discussed the LHDI in the exhaust of a symmetric reconnection region. The lower panel shows the electron beta, with

*β*being smallest in the upper-right quadrant, where the velocity signatures are seen. This is consistent with the LHDI, where reducing

_{e}*β*makes conditions more favorable for the instability.

_{e}^{1,5}

A more detailed illustration of the instability is shown in Fig. 2. Here, the electric field in the *x*–*y* plane at *z* = 101.8 is plotted. The top panel shows *E _{x}*, while the middle panel shows the parallel electric field $ E \u2225 = E \u2192\xb7 b \u2192 \u0302$. The in-plane magnetic field is shown as streamlines in the

*E*plot. The electric field signatures show that there are actually two wave signatures, one associated with the region

_{x}*x*> 215 and one in the region $ 200 < x < 215$. In the

*x*> 215 region, there is a coherent

*E*structure with propagation almost perpendicular to the magnetic field. The structures in the $ 200 < x < 215$ region have a more pronounced $ E \u2225$ component and are oblique modes. There are also strong $ E \u2225$ signatures in the

_{x}*x*> 240 regions, which may be associated with streaming instabilities evolving into phase space holes (e.g., Refs. 28 and 29). Because of the guide field, the electrons are strongly magnetized in this region so that the perpendicular velocity is primarily due to $ E \u2192 \xd7 B \u2192$ drift. As such,

*E*and

_{x}*u*show similar structures in the

_{ez}*x*> 215 region. We note that Ref. 30 shows electron vortex-like structures attributed to the Kelvin–Helmholtz instability, but there is no electron velocity shear in the unstable region in this simulation, meaning that the mechanism for generating these waves is different.

The waves in the $ 200 < x < 215$ region have wavenumber $ k \rho e = 0.65$, and frequency $ \omega = 0.16 \omega p e = 2.17 \omega l h$, where *ω* is measured in the ion frame and *ω _{pe}* and

*ω*are calculated using local plasma and field parameters. These waves propagate approximately $ 100 \xb0$ from the magnetic field. In the

_{lh}*x*> 215 region, the waves have $ \omega = 0.019 \omega p e = 0.24 \omega l h$ in the ion frame with close to perpendicular propagation.

The three-dimensional structure of these waves is shown in the bottom panel of Fig. 2. To emphasize their structure, regions of $ | u e z | < 0.16$ are transparent, and the rendering shows the filamentary nature of the velocity perturbations superimposed on the current density in the *x*–*z* plane. Three-dimensional magnetic field lines are also shown, once again illustrating the almost perpendicular propagation.

In the rest of this work, we focus on the waves in the *x* > 215 region. We first study the dispersion relation of the waves using the model of Ref. 5, the wave characteristics are compared to the LHDI dispersion relation obtained from linear theory. This model includes electromagnetic effects and parallel electron heat flux. The results are shown in Fig. 3. As can be seen, the plasma is unstable in a narrow range of angles around the $ \theta = 90 \xb0$, which represents propagation perpendicular to the magnetic field, and the most unstable mode is around $ k \rho e = 0.7$. We note that for the waves in the *x* < 215 region, a dispersion analysis using a solver for homogeneous plasmas^{31} shows that when a beam population is present, oblique modes can be excited. Further work on the importance of electron beam populations on the LHDI dispersion where there is a density gradient is ongoing. An eigenmode analysis may also be necessary close to the current sheet where the width is ∼3 times the electron gyroradius.

The waves studied in the simulation have a longer wavelength than waves with the maximum growth rate in the model ( $ k \rho e = 0.48$ vs 0.7), but remain in the unstable region. The real frequency at $ k \rho e = 0.48$ is $ 0.25 \omega l h$, consistent with the measured frequency of the waves. This suggests that the longer wavelength mode may saturate at a larger amplitude. We note that LHDI theory with some finite-beta electromagnetic effects^{1} predicts that the most unstable waves have $ k \rho e \u2248 1$ for these parameters. The importance of the full inclusion of electromagnetic effects on the range of unstable *k* is discussed in Ref. 9.

*y*direction and fluctuating components as follows:

^{15,17}

*y*-averaged quantities and

*δ*represents fluctuating quantities.

*e*is the unit charge. Here, the term proportional to $ \u27e8 \delta n \delta E \u2192 \u27e9$ is usually referred to as anomalous resistivity, while the combination of the $ e \u27e8 \delta ( n u \u2192 ) \xd7 \delta B \u2192 \u27e9$ and $ m \u2207 \xb7 \u27e8 \delta ( n u \u2192 ) \delta u \u2192 \u27e9$ terms is referred to as the anomalous viscosity. Even though the wave propagation is mainly in the

*x*direction, averaging over the

*y*direction gives reliable results compared to averaging over a single wavelength in the direction of propagation. For example, with two simple plane wave quantities and nonzero

*k*,

_{x}*k*, $ f ( x , y ) = A \u2009 sin \u2009 ( k x x + k y y ) , \u2009 g ( x , y ) = B \u2009 sin \u2009 ( k x x + k y y + \varphi )$, we find $ \u27e8 f g \u27e9 y = k y / ( 2 \pi ) \u222b 0 2 \pi / k y f ( x , y ) g ( x , y ) d y$ = $ A B \u2009 cos \u2009 ( \varphi ) / 2$. This is equivalent to integrating in the direction of propagation $ k / ( 2 \pi ) \u222b 0 2 \pi / k A B \u2009 sin \u2009 ( k r ) \u2009 sin \u2009 ( k r + \varphi ) d r = A B \u2009 cos \u2009 ( \varphi ) / 2$.

_{y}As can be seen in Fig. 1, the structure of the current layer around the x-line is laminar, unlike what is seen in simulation studies using magnetopause-like parameters.^{15,16} This is likely due to the lower density ratio reducing the density gradients in the system, which drive the instabilities. The contribution of the anomalous terms to the reconnection electric field is negligible close to the x-line ( $ < 0.01 %$). Instead, since the waves are primarily directed in the *x* direction, we focus on the electron momentum equation in the *x* direction using the diagnostic in Eq. (3). This is shown in Fig. 4. In the region where the lower-hybrid waves are found, the non-ideal electric field is mainly balanced by the pressure tensor divergence, while the anomalous resistivity opposes the $ \u27e8 n \u27e9 e ( \u27e8 E \u2192 \u27e9 + \u27e8 u \u2192 \u27e9 \xd7 \u27e8 B \u2192 \u27e9 )$ term. However, its value in this region is approximately 10% of the non-ideal term, indicating that while the waves facilitate momentum transfer by reducing the electron outflow, the effect in the simulation is small. The relative contribution of the anomalous resistivity to the *y*-directed electric field has a similar scale of approximately 10 to 20% in the region where waves are found.

Finally, we proceed to study the effect of the LHDI on particle distributions. In Fig. 5, the electron distributions have been reduced to $ f ( x , v \u2225 )$, while the ion distributions have been reduced to $ f ( x , v x )$. Similar to other simulations and observations of asymmetric reconnection, the electrons show two populations: a core population with small positive $ v \u2225$ and a beam population with strong negative $ v \u2225$. The beam population is also anisotropic, with $ T \u22a5 > T \u2225$, and the velocity distribution functions are gyrotropic (not shown). The ion distributions show periodic temperature fluctuations with the same wavelength as the electric field due to the lower-hybrid waves, while the electron distribution shows a periodic opening (pink regions at $ v \u2225 \u2248 \u2212 0.25$) and closing (white shading at the same parallel velocity) of a phase-space hole in $ v \u2225$ correlated with the electric field.

*x*direction and the density gradient is in the

*y*direction. Here, $ \u03f5 = d n 0 / d y$. To calculate the density and

*p*response, velocity moments of $ ( f 1 + \delta f )$ are taken.

_{xx}*p*response, keeping only the terms proportional to

_{xx}*E*is

_{x}*A*is a positive constant, indicating that the temperature approximately is $ \u2212 \pi / 2$ out of phase with the electric field, consistent with what is seen in the simulation results, where the peaks of the electric field show an offset from the regions where the ion distribution is wider. In this calculation, we have used the measured real frequency in the simulation and the growth rate at $ k \rho e = 0.48$ from the linear theory. For simplicity, we have used the

*x*components of the quantities in the calculations as the waves are propagating only $ 9 \xb0$ from the horizontal.

The periodic phase space gaps in the $ f ( x , v \u2225 )$ distribution require a different explanation. Based on test particle tracing shown in Fig. 6, we find that beam particles (blue) are primarily streaming along the magnetic field in the negative *x* and *y* directions, while the core particles (red) are drifting in the positive *x* direction. We introduce a simple model to explain the effect on the distribution function.

We consider static electric and magnetic fields, with a wave propagating in the *x* direction, with $ E \u2192 = ( E 0 \u2009 sin \u2009 k x , E y , 0 )$ and $ ( B \u2192 = B 0 \u2009 sin \u2009 \theta , 0 , B 0 \u2009 cos \u2009 \theta )$, where *θ* is close to zero. This configuration represents an almost-perpendicular electrostatic wave in the frame moving at the phase velocity. For strongly magnetized particles, as in the simulation, the guiding center equations can be used. The lowest order equations of motion are $ v \u2192 \u22a5 = v \u2192 E \xd7 B , \u2009 m d v \u2225 / d t = e E \u2225 \u2212 \mu \u2207 B \u2212 m b \u2192 \xb7 d v \u2192 E \xd7 B / d t$. In this setup, only the $ E \u2225$ term is nonzero in the parallel equation of motion because $ B \u2192$ and *E _{y}* are constant. The only nonuniform contribution to $ E \u2192 \xd7 B \u2192$ is in the

*y*direction, which vanishes when the dot product with $ b \u2192$ is taken.

*ϵ*is an expansion parameter. $ v \u2225 , 0$ is then a constant, and the first-order equation is written as $ m d v \u2225 , 1 / d t = e E 0 \u2009 sin \u2009 ( k x ) \u2009 sin \u2009 \theta $. Using $ v x = v \u2225 , 0 \u2009 sin \u2009 \theta + v \u22a5 , x + \u03f5 v \u2225 , 1 \u2009 sin \u2009 \theta $, the first-order equation is expanded to give

*X*is the particle position at

*t*= 0. The constant of integration

*C*is chosen to make $ v \u2225 , 1 = 0$ at

*t*= 0 so that the initial parallel velocity is $ v \u2225 , 0$.

Physically, this gives the pattern seen in Fig. 5, where the maxima and minima of the separation in electron phase space are approximately $ \pi / 2$ out of phase with the extrema of *E _{x}*, particularly in the

*x*< 230 region. The physical explanation is that particles coming to a point in space from opposite sides of the wave front experience different phases of $ E \u2225$ and, hence, different accelerations. The two electron populations start to become out of phase for

*x*> 230, where the gaps appear more diagonal in $ v \u2225$–

*x*space. A possible explanation is that the coherent wave structure starts to break down for larger

*x*as shown in Fig. 2, particularly affecting the beam particles originating on the right. Note that the model breaks down at $ v 0 , x = 0$, which corresponds to particles moving exactly at the wave phase velocity.

## IV. DISCUSSION AND SUMMARY

The role of the lower-hybrid drift instability and lower-hybrid waves in reconnection regions has been the topic of numerous simulation and observational studies. Recent observations^{5,6,17} and laboratory experiments^{9} have shown the short-wavelength LHDI in and close to the electron diffusion region and have attempted to quantify the contribution of the lower-hybrid waves to the electron momentum equation (including the reconnection electric field) and electron acceleration. Magnetospheric Multiscale mission observations have shown that while the lower-hybrid waves contribute to particle transport and relaxing the density gradient across the magnetopause, their net contribution to the reconnection electric field is small.^{17} In cases with moderate guide fields, the lower-hybrid waves lead to the generation of vortices that demagnetize part of the electron population, causing non-gyrotropic electron acceleration to be seen.^{6} Laboratory experiments^{9} have shown that the conditions in the electron diffusion region are favorable for the excitation of the LHDI, and studies are still ongoing.

Three-dimensional simulations of the lower-hybrid drift instability have shown that for magnetopause conditions, the LHDI causes the development of a turbulent reconnection layer.^{14–16} This causes particle mixing and transport, though the contribution to the reconnection electric field by the lower-hybrid waves has been shown to be small. While anomalous resistivity and viscosity have been seen, they are due to longer wavelength modes causing kinking of the current sheet rather than the lower-hybrid waves.^{15} There are also previous simulations of antiparallel reconnection in MRX^{32} in which longer wavelength $ k \rho e \rho i \u223c 1$ fluctuations are found. In that case, the electron-scale modes are stabilized due to the high beta near the center of the current sheet since there is no guide field.

In the context of these results, we note that the inclusion of a strong guide field also means that the LHDI modes are stabilized near the x-line and primarily propagate in the *x* direction in the lowest *β* outflow region. These waves do not contribute to the reconnection electric field close to the x-line, though they do affect momentum balance in both *x* and *y* directions in the unstable region, acting to reduce the electron outflow velocity. Similar to the earlier simulations, the effect of the anomalous terms is small.

The lower-hybrid waves have an effect on the electron distributions in the reconnection region. As seen in Fig. 5, there are two electron populations, as expected from an asymmetric reconnection configuration.^{14,15,33} Because the waves are not purely perpendicular to the magnetic field, there is a small $ E \u2225$ component that causes periodic opening and closing of a phase space gap between the populations in the $ f ( x , v \u2225 )$ distribution.

One major difference between the simulations and observations or laboratory experiments is the lack of electrostatic lower-hybrid waves close to the x-line even in the strong guide field regime. Earlier simulations show strong LHDI at the separatrix,^{14,16} or in the outflow regions,^{27} but only the longer wavelength electromagnetic modes are seen inside the current sheet.^{13,14} In this simulation, lower-hybrid waves are seen throughout the outflow region, unlike in Ref. 27, where they are confined to regions with enhanced outflow velocity. It is possible that the numerical parameters, such as the reduced mass ratio, are affecting the measurement of the waves in simulations. In observations and experiments, the perpendicular bulk electron velocities regularly exceed the local ion Alfvén speed *v _{A}*, while in the simulations, they are smaller than

*v*. From the dispersion, the lower normalized perpendicular velocity leads to a lower growth rate of the LHDI due to the reduced amount of free energy present.

_{A}^{5}The amplitude of the fluctuations is also affected, with $ \Delta B / B \u2248 1 %$ in this simulation, and $ < 5 %$ in Ref. 27, while they are comparable to the guide field magnitude in Ref. 6. We do note that the parameters in Ref. 27 are more similar to the event in Ref. 6, where the guide field was smaller.

Because these simulations require three-dimensions to capture the LHDI, reduced parameters are necessary due to computational limitations, artificially increasing the Alfvén speed and thermal velocities. It is possible that this is affecting where the lower-hybrid waves are found in simulations and their effect on the plasma. Additionally, the artificial $ \omega p e / \omega c e$ ratio effectively increases the temperature of the plasma species, which may also affect the generation of waves. Preliminary results can be found in the Appendix, and further study on this issue will be required.

In conclusion, we have performed a three-dimensional simulation of asymmetric reconnection with a strong guide field based on MRX experiments. We show that the electrostatic LHDI is found in the low-beta outflow region as predicted by theory, though the contribution to momentum balance and the reconnection electric field is small. However, the lower-hybrid waves do modify the particle distribution functions, with the most striking result being a periodic opening and closing of a phase space hole in the electron distribution. The lack of lower-hybrid waves close to the x-line in simulations compared to experiments and observations raises questions about the parameters used to study the LHDI in reconnection simulations.

This work was supported by NASA Grant Nos. 80NSSC21K1462, NNH20ZDA001N, and 80NSSC21K1795.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Jonathan Ng:** Conceptualization (equal); Formal analysis (equal); Methodology (equal); Software (equal); Visualization (equal); Writing – original draft (equal). **Jongsoo Yoo:** Conceptualization (equal); Funding acquisition (lead); Writing – review & editing (equal). **Li-Jen Chen:** Supervision (equal); Writing – review & editing (equal). **Naoki Bessho:** Writing – review & editing (equal). **Hantao Ji:** Project administration (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are openly available in NASA HECC, Ref. 34.

### APPENDIX: SCALING OF GROWTH RATE WITH PARAMETERS

Here, we present the results of changing the parameters $ u / v A , m i / m e$, and $ \omega p e / \omega c e$ in the LHDI calculation using the dispersion relation in Ref. 5. These results are shown in Fig. 7. When increasing $ u e x / v A$, keeping all other parameters unchanged, the growth rate increases.

Changing $ m i / m e$ affects the value of *v _{A}*. As such, we provide results where the unnormalized value of

*u*is kept constant and where $ u / v A$ is kept constant by changing the electron velocity. When $ u / v A$ is kept constant, the growth rate remains approximately constant. Otherwise, $ u / v A$ increases as the ion mass is increased, leading to an increased growth rate.

In these calculations, changing $ \omega p e / \omega c e$ effectively scales the value of $ | B |$. This affects other derived quantities, such as $ v A , \rho e$, and *β*. In this case, we keep *β* constant by scaling the electron and ion temperatures. When $ u / v A$ is constant, the growth rate decreases as larger (more realistic) values of $ \omega p e / \omega c e$ are reached. When $ u / v A$ is not kept constant, the growth rate increases with larger $ \omega p e / \omega c e$.

In a reconnecting system, we do not expect $ u / v A$ to remain constant as the electron outflows can reach fractions of the electron Alfvén velocity. However, the outflow can be affected by the pressure anisotropy and guide field; therefore, a conclusive answer of how these parameters affect the LHDI in such a system cannot be given here.