Bifurcated Current Sheet Observed on the Boundary of Kelvin-Helmholtz Vortices

On May 5, 2017 MMS observed a bifurcated current sheet at the boundary of Kelvin-Helmholtz vortices (KHVs) developed on the dawnside tailward magnetopause. We use the event to enhance our understanding of the formation and structure of asymmetric current sheets in the presence of density asymmetry, flow shear, and guide field, which have been rarely studied. The entire current layer comprises three separate current sheets, each corresponding to magnetosphere-side sunward separatrix region, central near-X-line region, and magnetosheath-side tailward separatrix region. Two off-center structures are identified as slow-mode discontinuities. All three current sheets have a thickness of ∼0.2 ion inertial length, demonstrating the sub-ion-scale current layer, where electrons mainly carry the current. We find that both the diamagnetic and electron anisotropy currents substantially support the bifurcated currents in the presence of density asymmetry and weak velocity shear. The combined effects of strong guide field, low density asymmetry, and weak flow shear appear to lead to asymmetries in the streamlines and the current-layer structure of the quadrupolar reconnection geometry. We also investigate intense electrostatics waves observed on the magnetosheath side of the KHV boundary. These waves may pre-heat a magnetosheath population that is to participate into the reconnection process, leading to two-step energization of the magnetosheath plasma entering into the magnetosphere via KHV-driven reconnection.


INTRODUCTION
About sixty years ago Dungey (1961) and Axford and Hines (1961) proposed two different models of the solar wind-magnetosphere interaction. The former was based on the concept of magnetic reconnection under large magnetic field shear. The latter was on quasi-viscous interaction in the boundary layer powered by large flow velocity shear. Since then two most important physical processes that lead and regulate the solar wind-Earth's magnetosphere coupling are thought to be magnetic reconnection and the Kelvin-Helmholtz instability (KHI).
Both processes exhibit multi-scale features and either compete or enhance each other. Magnetic reconnection is initiated on the electron-scale size, i.e., in the electron diffusion region (EDR) and then entails dynamics in the ion diffusion region (IDR), and ultimately propagates its effect to the macroscopic region where magnetohydro-dynamics (MHD) governs. On the other hand, Kelvin-Helmholtz waves (KHWs) occurring on the Earth's magnetopause are often generated on the macroscopic scales, i.e., ∼1 R E (Earth radii) (Hasegawa et al., 2004) and then involve kinetic processes occurring on the ion and electron scales as the waves nonlinearly grow into Kelvin-Helmholtz vortices (KHVs). The vortex motion facilitates the formation of thin current sheets between stretched magnetosheath and magnetospheric field lines at the edge of KHVs where magnetic reconnection can occur Otto, 2001, 2004). Observations (Fairfield et al., 2000;Nykyri et al., 2006;Eriksson et al., 2009Eriksson et al., , 2016Hasegawa et al., 2009;Li et al., 2016;Hwang et al., 2020;Kieokaew et al., 2020) have reported ongoing reconnection in such a thin current sheet developed along the boundary of KHVs or a magnetic island as predicted by simulations (Nakamura et al., 2013(Nakamura et al., , 2017. Magnetic reconnection at the edge of KHVs and/or a wave packet inside KHVs (Moore et al., 2016) result in cross-scale energy transport. Cassak and Otto (2011) showed that a flow shear across the symmetric reconnection current sheet decreases the efficiency of the reconnected field line to drive the outflow, similarly to the suppression of reconnection by the diamagnetic effect (Swisdak et al., 2003(Swisdak et al., , 2010. They used the full particle simulation to derive that the reconnection-cutoff velocity shear is the upstream Alfvén speed. Doss et al. (2015) studied the effect of the flow shear in asymmetric reconnection, analytically and numerically predicting that the asymmetric effect allows reconnection to continue even for super-Alfvénic upstream velocity shear. Tanaka et al. (2010) further considered the effect of a guide field as well as a flow shear in asymmetric reconnection, reporting that both an initial upstream flow and the Lorentz force acting inflowing plasmas in the presence of a guide field produce a slanted inflow to the current sheet. Resulting asymmetries in the quadrupolar reconnection current layer is qualitatively similar to the MHD simulation by La Belle-Hamer et al. (1995).
The 2-D MHD simulation for the current sheet across which substantial velocity shear and density jump exist (La Belle-Hamer et al., 1995) indicates that depending on either the competition (occurring on the tailward exhaust region) or the enhancement (sunward exhaust) of the two velocity-shear and densty-asymmetry effects, the structure of the current sheet is often different from the simple 1-D Harris model, showing double off-center peaks in current, i.e., current bifurcation. The bifurcated current sheet has been understood as the Petschek-type reconnection layer, where the reconnection outflow jets ejected from the X-line are bounded by two rotational discontinuities or slow mode shock structures, which split a single reconnecting current sheet.
The formation and structure of asymmetric current sheets in the presence of flow shear, density asymmetry, and guide field have been much less studied. In particular, despite the prediction by La Belle-Hamer et al. (1995), evidence of bifurcated current sheets in KHW/KHV-induced reconnecting layers is rarely reported to this day. The only observation by Cluster  showed that each of the two current sheets constituting the bifurcated layer had a thickness less than the ion inertial length and that the current was likely supported by electrons.
To enhance our understanding of the properties of the magnetic reconnection layer under the combined sheared plasma flow, guide field, and density asymmetry, i.e., typically occurring on the flank-side magnetopause, we use the data from MMS on May 5, 2017. In this paper, we present the observation of a bifurcated current sheet identified on the boundary of KHVs. The following paragraph briefly describes the MMS instruments and data analysis techniques used for the present study (Methods Section). We then investigate plasma and field properties associated with the bifurcated and central current sheets and show that the electrons drifting under both the diamagnetic effect and the magnetic curvature with large temperature anisotropy significantly contribute to the current (The Structure of Current Sheets Section). We also investigate intense electrostatics waves that are predominantly observed on the magnetosheath side of the central current layer (Wave Observation and Analysis Section). Discussion of the formation and structure of the observed current sheet in the presence of flow shear, density asymmetry, and guide field, and the implied cause and effect of the enhanced waves follow in Discussion Section.

METHODS
The four MMS spacecraft (Burch et al., 2016a) fly in lowinclination and highly elliptical orbits. We used the magnetic field data with a time resolution of 10-ms in burst mode, the electric field data with a 0.122-ms time resolution in burst mode, and ion and electron data in burst mode with a 150-ms and 30-ms time resolution, respectively, a 11.25°angular resolution, and an energy range of ∼10 eV-26 keV.
We determined boundary normal coordinates (LMN) by performing minimum directional derivative (MDD) analysis (Shi et al., 2005): three eigenvectors corresponding to the medium, minimum, and maximum eigenvalues (λ) of the matrix, (∇B)(∇B) T constitute the l, m, and n axes, respectively, in the LMN coordinates. To determine the propagation velocity of the Frontiers in Astronomy and Space Sciences | www.frontiersin.org November 2021 | Volume 8 | Article 782924 current layer, we performed spatio-temporal difference (STD) analysis (Shi et al., 2006). The current-sheet-normal propagation velocity is consistent with the value calculated from a fourspacecraft timing analysis (Paschmann et al., 1998). To investigate the wave propagation of the electrostatic waves, we used the maximum variance analysis (Sonnerup and Scheible, 1998;Siscoe and Suey 1972) of the electric field. To further investigate the wave mode, frequency, and growth rate, we performed the linear kinetic instability analysis using BO code (Xie, 2019) with input parameters obtained by fitting the observed ion and electron distributions functions to a sum of multi-component Maxwellian distributions.

THE STRUCTURE OF CURRENT SHEETS
From 1920 to 2320 UT on May 5, 2017, MMS observed quasiperiodic perturbations of the dawnside tailward magnetopause, as reported by Hwang et al. (2020). Figure 1A shows (i) the interplanetary magnetic field obtained from ACE OMNI-HRO 1-min data and (ii) the magnetic field at dawnside tailward magnetopause encountered by MMS4 (x, y, and z components in blue, green, and red with the magnetic strength in black) during 1930-2100 UT in Geocentric Solar Magnetospheric (GSM) coordinates. The ACE HRO data provide the timeshifted IMF (interplanetary magnetic field) at a model bow shock nose location (Russell et al., 1983). The event occurred within a period of mainly northward and slightly sunward/dawnward IMF (Figures 1Ai). MMS4 observed quasi-periodic fluctuations with a period of ∼2.5-6 min (Figures 1Aii). Hwang et al. (2020) showed via boundary-normal analyses that the fluctuations were likely to be attributed to nonlinear KHWs. During this internal, we identified a thin current sheet formed at the boundary of the KHV at ∼2009:49 UT (the vertical blue line in Figure 1A) when MMS traversed the boundary from the magnetospheric side to the magnetosheath side. Figures 1Bi shows the x component of the magnetic field (B) in GSM coordinates. The negative-to-positive reversal of B x observed by MMS1, 2, 3, and 4 (black, red, green, and blue) indicates a current sheet. The B x profiles are, however, different from a Harris-sheet hyperbolic tangent profile, displaying local dip and peak and plateau (MMS2) or gentle slope (MMS134) between them, indicative of a bifurcated current sheet. The medium-to-minimum (maximum-tomedium) eigenvalue ratio is ∼14.6 (6.4), indicating a reliable calculation. The MDD-derived LMN coordinates are consistent with the LMN coordinates derived from minimum variance analysis (MVA) (Sonnerup and Scheible, 1998; Siscoe and Suey 1972) within 12.5°, 10.9°, and 6.7°differences along l, m, and n for all four spacecraft. We use the LMN coordinates throughout following figures and analyses. Figure 2A shows the tetrahedral configuration of the four MMS spacecraft in LMN around its barycenter at (−13.9, −17.9, −4.8) GSM Earth radii (R E ). The notable difference in B x between MMS2 and MMS134 (Figures 1Bi) most likely came from the n-directional separation, as seen in their LN-plane projections. About 179 km separation along n as well as the average spacecraft separation of ∼156 km/s are comparable to the ion inertial length (d i ) based on the magnetosheath values (d i ∼185 km).
Figures 1Bvii shows the tetrahedral-averaged B in LMN. The negative-to-positive B l reversal is denoted by "C2" and the vertical dashed gray line. The tetrahedral-averaged electric current calculated from particle moments perpendicular to B (Figures 1Bviii) shows three overall peaks between "L" and "T". They comprise two larger peaks before and after "C2" and a smaller peak at ∼"C2". Using the STD-driven current-sheet normal velocity ( We assume the overall currentsheet-normal velocity to be 67 km/s along -n. Then, the three current sheets before, at/around, and after "C2" with a duration of ∼0.65, 0.50, and 0.65 s (Figures 1Bviii and Figure 4B) has a thickness of ∼43.6, 33.5, and 43.6 km. Since these values correspond to ∼0.24, 0.18, and 0.24 d i (∼9.7, 7.4, and 9.7 electron inertial length, d e ∼ 4.5 km in this event) similar to Hasegawa et al. (2009), the current in these sub-ion scale current sheets is expected to be supported by electrons.
Due to the large spacecraft separation compared to the current sheet thickness, investigation of the detailed structure of the current layer should be made using an individual spacecraft. We use the data from MMS1 (Figure 3 with all vector parameters in LMN). Figures 3A,B shows the l (blue), m (green), and n (red) components of the magnetic (B) and electric (E) fields. The leading and trailing edges ("L" and "T", magenta and red dashed lines) denote dip and hump in B l , decreases in B m , and increase (ΔE n ∼8.5 mV/m) and decrease (ΔE n ∼−2.0 mV/m) in E n at "L" and "T", respectively. (The latter two signatures correspond to the Hall features as illustrated in Figure 2B to be discussed in the following paragraphs) The reversals in B l and B n are marked by "C" and "C*" (black and gray dashed lines), respectively.
Variations of the ion density and ion/electron temperatures (Figures 3C,D) together with ion/electron energy spectrograms (Figures 3H,I) show that MMS1 crossed the current sheet from the more magnetospheric region (prior to "C") to the more magnetosheath region (after "C"). Intense electric field fluctuations (marked by "TB" and two vertical dashed cyan lines) are seen in the magnetosheath side of the current sheet (to be discussed in Wave Observation and Analysis Section).
The ion velocity between "L" and "T" varies from slower tailward flow (smaller −V i,l ) to faster tailward flow (larger −V i,l ) across "C" (marked by the blue arrow in Figure 3E) around V i,l −154 km/s (the blue dotted line). This indicates the sunward exhaust region (before "C") to the antisunward exhaust region (after "C") of the current sheet, which was convecting antisunward along with the KHV propagation.
Therefore, MMS most likely crossed the overall current sheet from the sunward magnetospheric quadrant to the antisunward magnetosheath quadrant of the reconnection plane with a large guide field, B g (B m ) ∼ 1.5 |B l | (at 2009:46.0 UT; Figure 3A) out of the plane. The trajectory of MMS is denoted by the dashed orange Frontiers in Astronomy and Space Sciences | www.frontiersin.org November 2021 | Volume 8 | Article 782924 arrow in Figure 2B, where "L", "C", "C*", and "T" correspond to those shown at the top of Figure 3A. The aforementioned Hall magnetic and electric field signatures are illustrated in green and red, respectively, around "L" and "T".
The electron velocity shows more complicated patterns than the ion velocity. Beyond the same variation of the slower-to-faster tailward outflow jets across "C" (blue arrow in Figure 3F), the slower tailward flow between "L" and "C" includes a short duration of the sunward flow (+V e,l ) during 2009:47.9-48.1 UT (magenta arrow). This indicates the existence of a narrow (∼13.4 km, ∼3 d e ) electron-current layer embedded in the outflow region (marked by "V e channel" in Figure 2B). Its counterpart may exist in the tailward exhaust region (between "C" and "T") with a faster tailward jet before "T" (magenta arrow in Figure 3F). Figure 2B shows possible electron flow streamlines in dashed blue arrows that may explain the observed flow channels.
Around "L" and "T", V e,l sharply changes its sign. The enhanced tailward flow before/at 'L' and the sunward flow at/ after "T" (black arrows in Figure 3F) are associated with electrons streaming toward an X-line in the separatrix region (see solid blue arrows in Figure 2B) (Egedal et al., 2005(Egedal et al., , 2008Hwang et al., 2017Hwang et al., , 2018. Pitch-angle distributions of the low-(<200 eV) and mid-(200 eV < energy <2 keV) energy electrons support this, showing the enhancement of the parallel flux at/around "L" and "T" (red arrows in Figures 3J,K). These counter-streaming electron flows (±V e,l ) across ∼"L" and ∼"T" sustain the Hall field along the separatrix.
These electron populations carry the electric current (current density, J) around "L" and "T". Figure 4B shows the l, m, and n components of J calculated from both ion and electron moments (solid blue, green, and red profiles). Overplotted are the ion current (dot-dashed light blue, light green, and orange) and the electron current (dotted blue, dark green, and red). The current (in particular, J m ) is mostly carried by electrons. Both J l and J m between "L" and "T" show the three-peak structure with two larger peaks before/after "C" and a smaller peak located at the B n reversal ("C*"), as demonstrated in J || (black arrows in Figure 4C). Thus, we speculate that the two larger peaks correspond to one of each pair of a bifurcated current sheet and the central peak is associated with an X-line ( Figure 2B).
In the sub-ion scale current layers such as this event, ion velocities perpendicular to B can be different from the E × B drift while electrons mostly follow E × B. Figures 4D-F shows ion (red) and electron (blue) velocities perpendicular to B compared with the E × B drift (black). Ion perpendicular velocities relatively agree with the E × B trend, but showing a substantial deviation from E × B around "L" and "T". Electrons show a more notable deviation from E × B during "L"− 'T' as denoted by yellow arrows in Figures 4D-F. This can result from a certain level of electron agyrotropy or other perpendicular drifts such as diamagnetic and/or magnetic curvature drifts (Norgren et al., 2018).
To see the level of electron agyrotropy, we use Q √ that quantifies the level of agyrotropy (Swisdak, 2016) as shown in green in Figure 4G: 0 for gyrotropy and 1 for maximal agyrotropy. In general, the agyrotropy is weak, showing a bit higher level of Q √ in the magnetospheric side than the magnetosheath side, as predicted by a higher temperature and lower density for magnetospheric electrons ( Figures  3C,D). Local peaks around "L", "C" and "T" are insignificant. On the other hand, Joule dissipation in the electron frame, J · (E + V e × B) shown in black, blue, and red profiles representing the total, parallel, and perpendicular components to B, shows fluctuating or positive values between "L" and "T". The dissipation (mostly along B) is enhanced during a later half of "TB", where intense wave activities are found. To understand the electron deviation from E × B and the origin of a pair of off-centered (bifurcated) currents, we plot the l, m, and n components of the measured electron perpendicular current, J e,⊥ (black profiles in Figures 4H-J), compared with those of the electron E × B current (−en e E×B B 2 , blue), the electron diamagnetic current ( B×∇P e,⊥ B 2 , orange), and the electron anisotropy current taking into account the influence of curvature drifts (P e,|| − P e,⊥ ) B×(B·∇)B B 4 , green, where n e is the electron density, and P e,|| and P e,⊥ are the electron pressures parallel and perpendicular to B (Zelenyi et al., 2004). Due to the large spacecraft separation (∼d i ), we cannot calculate gradient terms using four-spacecraft measurements. Instead, we use ∇ ≈ 1/(dt V MMS across structure ) −1/(dt V structure ), where dt is a sampling cadence of the electron data. The normal component of V structure was derived from MDD as described earlier. The tangential (l and m) components of V structure are, however, largely uncertain for the 1-D structure (Figures 1Biii). We use the l and m components of the background ion bulk velocity before "L" (averaged for 2009:46.0-47.5 UT), giving rise to V structure (−173, 61, −67) km/s in LMN. The red profiles in Figures 4H-J are the sum of the E × B, diamagnetic, and anisotropy currents, and show better agreements with the measured J e,⊥ than each of the three contributions. We find that both the diamagnetic and anisotropy currents significantly contribute to the current at/around "L". The diamagnetic current predominantly supports the current at/around "T" or between "C*" and "T". The anisotropy effect is most dominant in the n component of J e,⊥ and between "L" and "C" (Note a large electron temperature anisotropy between "L" and "C" in Figure 3D).

WAVE OBSERVATION AND ANALYSIS
We investigate the intense waves observed intermittently within the current layer (marked by 1, 2, and 3 in Figure 4K and stronger wave activities (throughout 4-7) observed during "TB". Figures 4L,M show the power spectral density (PSD) of E and B. The waves are mostly electrostatic and enhanced near or below the electron cyclotron frequency (f CE ) or the ion plasma frequency (f PI ) and above the lower-hybrid frequency (f LH ). Figure 5i shows the waveform of E decomposed into parallel (red) and perpendicular (blue) components with respect to B for timing 1, 4, and 7 ( Figure 4K). To estimate the propagation direction (k) of the electrostatic waves, we use maximum variance analysis of E for each interval. Results shown in Table 1 (a) demonstrate that these waves propagated parallel or antiparallel to B.
To understand the generation of theses waves, we perform linear instability analysis using the electron and ion distribution functions at timing 1, 4, and 7. Each particle distribution is modeled by a sum of two Maxwellian distributions, and the best fitting parameters (density, thermal speed, and beam drift speed) are listed in Table 1 (b-g). Figure 5ii-v shows MMS observations of 2-D reduced electron (ii) and ion (iii) distributions and their 2-D reduced model distributions (iv, v) Figures 5vi-vii show 1-D reduced distributions in the V || axis for detailed comparisons between model (red) and observation (black). The modeled distributions agree well with the MMS observation for all timing 1, 4, and 7. We note that a cold ion population exists throughout these times and bidirectional electron beams exist in timing 1 and 4, but are flattened at timing 7.
By making use of these modeled-distribution parameters ( Table 1 b-g), we perform the linear kinetic instability analysis using BO code (Xie, 2019). Figures 5viii-ix show real and imaginary parts of the growing mode in ω − k space. The frequency ω and the growth rate c are normalized to the ion plasma frequency ω pi , and the wave number k is normalized to k 0 ω pi /(10 3 km/s). The ω − c plot in Figure 5x shows more clearly in which frequency range waves are generated.
At timing 1, a low-frequency mode is generated in the range of 0.2 − 0.3ω pi as well as a broader spectrum in the range of 0 − 1.2ω pi . The low-frequency mode is the fastest growing mode, which is observed only at/around timing 1 ( Figure 4L). The phase speed of the fastest growing mode is ∼130 km/s at the maximum growth rate, less than the ion acoustic speed (∼250 km/s).
At timing 4, the frequency range of wave generation is much broader than timing 1. Two distinct modes are derived. One locates below ω pi with a peak at ∼ 0.6ω pi . The other locates in the range of 0 − 2.0ω pi and peaks at ∼1.0ω pi . Their phase speeds are ∼300 km/s and ∼420 km/s, respectively. Because the frequency ranges of the two modes are overlapped as well as their growth rates are comparable to each other, the two wave modes may not be distinguished in the observation. The superposition of these waves might explain that the waveforms (Figure 5Bi) slightly deviate from sinusoidal. At timing 7, the modeled distribution produces no growing mode most likely due to the flattened electron distribution. This indicates that the bi-directional electron beams are a major freeenergy source for the generation of the observed electrostatic waves.

DISCUSSION
In this paper, we report a bifurcated current sheet developed on the boundary of KHVs propagating along the flank-side magnetopause, across which both plasma flow shear and density asymmetry exist under a large guide field, B g ∼ 1.0 |B l | (on the magnetosheath side) to 1.5 |B l | (on the magnetospheric side). Via discussion on The Structure of Current Sheets Section, we speculate the trajectory of MMS that followed the dashed orange arrow in Figure 2B across the reconnection plane.
The overall current density profiles show three peaks (Figures 4K,L; green shades in Figure 2B), each observed in the proximity to the magnetospheric-side, sunward separatrix region (around "L"), the central, near-X-line region ("C-C*"), and the magnetosheath-side, tailward separatrix region (around "T"). The slower-tailward to faster-tailward jets across the central current sheet, i.e., reconnection outflows, demonstrate that the two off-centered signatures are corresponding to two rotational discontinuities or slow mode shocks in the Petschek reconnection geometry. 1) Tangential (B l ) and normal (B n ) components of B are non-zero at/around "L" and "T". 2) Decrease in |B l | from upstream (inflow region) to downstream (outflow region) of "L" and "T" (along magenta and red arrows in Figure 3A) indicates that the magnetic field bends toward n.
3) The magnetic field strength or pressure decrease from upstream to downstream (magenta and red arrows in Figures 3A,G). 4) The plasma density and pressure increase (magenta and red arrows in Figures 3C,G) across "L" and "T". All 1-4) features support that the two discontinuities are identified as slow modes.
For the two periods between "L" and ∼"C" and between ∼"C" and "T", we performed a Walén test separately for ions and electrons (Scudder et al., 1999;Hwang et al., 2016). The ion flow in the deHoffmann-Teller frame (V HT ) showed a linear correlation with the ion Alfvén velocity with a correlation coefficient of ∼0.8 for both intervals, but only 0.1−0.2 of the ion Alfvén speed; electrons did neither satisfy the Walén relation 1 | Parameters of modeled distributions used for linear analysis, where n 0 is total density and c is the speed of light. Each distribution is modeled as a sum of two Maxwellian distributions.
(a)k in LMN, angle between (k, B), max-to-mid eigenvalue ratio of MVA nor display a correlation (not shown). This indicates that 1) the reconnection layer has not been fully developed (equivalently, the MMS orbit was too close to the X-line) or 2) other accelerations exerted on the current sheet. The sub-ion scale current sheet (The Structure of Current Sheets Section) and the significant contribution from the diamagnetic drift and/or the electron anisotropy drift (Figures 4H-J) indicate these possibilities 1-2), respectively. Lin and Lee (1994), for asymmetric guide-field reconnection with no velocity shear, predicted the formation of different discontinuities between on the magnetosheath-side separatrix region (a time-dependent intermediate shock and a slow expansion wave, which evolves to a slow shock with time) and the magnetospheric-side separatrix region (a time-dependent intermediate shock and a weak slow shock). Such multiple discontinuities were not observed in this event, where the whole current layer is on a sub-ion scale, i.e., possibly due to the MMS trajectory being too close to the X-line. Further comparison is hindered since MMS did not traverse the entire exhaust region either side of X. The MHD simulation (La Belle-Hamer et al., 1995) for asymmetric no-B g reconnection with a flow shear also predicted the formation of an intermediate shock on the magnetosheath-side, sunward separatrix region (the upper-left quadrant of Figure 2B), which was, again, not traversed by MMS.
We note a short duration of the sunward electron jet between "L" and "C". V e,n and V i,n are more negative and less negative across "C" (red arrows in Figure 3F). Thus, the plasma streamlines between "L"−"C" and "C"−"T" might not be symmetric (dashed blue arrows in Figure 2B).
La Belle-Hamer et al. (1995) and Tanaka et al. (2010), indeed, predicted such asymmetry in the reconnection geometry under the density gradient and velocity shear. In the tailward exhaust region (i.e., between "C*" and "T"), the outflow (−L) is in the same direction as the upstream magnetosheath flow (−L). Thus, a smaller force is required to drive the outflow. On the other hand, the larger density/inertia on the magnetosheath side requires a larger accelerating force to drive the outflow. As a result, the effects of shear flow and density gradient compete with each other, which results in the streamlines less deformed (blue dashed arrows in Figure 2B) and makes the field transition layer broader (note that the structure was significantly 2-D toward "T" in Figures 1Biii) as predicted by Figure 4 of La Belle-Hamer et al. (1995) and Figure 7 of Tanaka et al. (2010).
In the sunward magnetosheath-side exhaust region, the outflow (+L) is opposite to the upstream magnetosheath flow (−L), requiring a larger accelerating force to drive the outflow. The shear-flow and density-gradient effects enhance each other, forming a narrow field reversal region and putting the accelerated flow on the magnetospheric side of the field reversal ( Figure 4 of La Belle-Hamer et al., 1995;Figure 7 of Tanaka et al., 2010). We speculate that the observed narrow sunward electron jet on the sunward, magnetospheric quadrant (between "L" and "C") is consistent with this prediction.
It may be notable that although such asymmetric streamlines are indicated by ion flows in the MHD (La Belle-Hamer et al., 1995) and particle-in-cell (Tanaka et al., 2010) simulations, the electron velocity appears to mostly represent the asymmetry in the present event. This implies that the aforementioned combined effects of the shear flow and density asymmetry are valid for the electron streamlines, in particular, in this sub-ion scale current layer.
We also note that the upstream flow difference across the current sheet is quite weak in this event (∼6% of the parallel Alfvén speed on either side of the current sheet) while B g is strong. Tanaka et al. (2010) showed that the combination of density gradient and guide field led to the similar effect obtained by the combination of density gradient and flow shear. Thus, we conclude that the combined effects of strong guide field, low density asymmetry (ρ sh ρ sp ∼2.2), and weak flow shear appear to derive asymmetries in the streamlines and the current-layer structure of the quadrupolar reconnection geometry, as illustrated in Figure 2B.
We estimate how these asymmetries would modify the reconnection rate, using the formula derived by Doss et al. (2015): where E 0, asym is the asymmetric reconnection rate in the absence of upstream shear flow (Cassak and Shay, 2007), and V A is the hybrid Alfvén speed, where ρ sh and ρ sp (B sh and B sp ) represent the magnetosheath and magnetospheric mass density (magnetic field intensity). We use the magnetospheric and magnetosheath (upstream) values obtained at ∼2009:47.0 UT and at ∼2009:50.2 UT, respectively, which give the velocity difference, ΔV e,l ∼124 km/s across the current sheet and the density ratio, ρ sh ρ sp ∼2.2. Our estimate of the second term in the parenthesis of Eq. 1 is ∼ −0.003. Thus, the effect of the combined velocity shear and density asymmetry will have little influence in the asymmetric reconnection rate in this event. At the magnetotail flanks, the current sheet developed along the KHV boundary will mainly move with the bulk tailward velocity of the KHV. The velocity shear (ΔV e,l ) is expected to be small, significantly reduced from the initial upstream velocity shear, which is the case as shown in Figures 3E,F. Electrons mainly carried the current for the present event, and ion contribution to the currents is limited up to ∼27% of the total current ( Figure 4B), which is expected for the sub-ion scale current sheet. The three current density humps have a thickness of ∼43.6, 33.5, and 43.6 km, i.e., ∼0.24, 0.18, and 0.24 d i (∼9.7, 7.4, and 9.7 d e ), respectively, demonstrating the sub-ion scale current layer.
Numerous theoretical and simulation studies for the magnetotail (i.e., symmetric) environment have been performed to understand the formation of the current sheet bifurcation. Among various mechanisms proposed, one important factor is temperature anisotropy. Sitnov et al. (2004) suggested that the bifurcation is caused by weak ion temperature anisotropy with T i,⊥ 1.1−1.2 T i,|| . Zelenyi et al. (2004) and Jiang and Lu (2021) suggested that the bifurcation can be caused by the electron Frontiers in Astronomy and Space Sciences | www.frontiersin.org November 2021 | Volume 8 | Article 782924 pressure anisotropy (P e,⊥ > P e,|| ), which decreases the current sheet density at the center of the current sheet, via the electron anisotropy drift contribution. Schindler and Hesse (2008) also showed P e,⊥ > P e,|| during the formation of a bifurcated current sheet (with a half-thickness of ∼d i ) embedded in an initially wider (∼5d i ) current sheet under quasisteady compression. In our observation, we note the opposite electron anisotropy, P e,|| > P e,⊥ throughout the current layer and most enhanced in the magnetospheric-side, sunward-exhaust region. This results in a significant contribution of the electron anisotropy current in supporting the bifurcated current along n direction ( Figure 4J). A larger contribution from the diamagnetic current was observed in the magnetosheath-side, tailwardexhaust region (Figures 4H-J). Therefore, both the diamagnetic and electron anisotropy currents substantially support the bifurcated currents in the presence of density asymmetry and velocity shear.
A statistical study of the bifurcated current sheets using Cluster data (Thompson et al., 2006) indicated that the narrower the current sheets are, the more likely they are bifurcated. This, together with our present study, may suggest that the electrons play a major role in the formation of the bifurcated current sheet in both symmetric and asymmetric environment.
We investigated intense electrostatics waves that were predominantly observed on the magnetosheath side of the central current layer using linear kinetic analysis for selected timings, 1, 4, and 7 ( Figure 4K). At timing 1 and 4, the electron distributions contain clear bi-directional beams with growing wave modes produced, while at timing 7 they show a plateau distribution with no growing mode wave. Still, large differences in the wave generation between timing 1 and 4 imply that various types of waves could be generated by the bi-directional beams that are ubiquitous in the KHV-induced reconnection sites Hwang et al., 2020).
Unlike the linear kinetic instability theory, we observed the electrostatic wave at timing 7. The waveforms at timing 7 (Figures 5Ci), however, indicate highly nonlinear waves. They are possibly propagated to the MMS location, after having been generated remotely. We also note that the first wave signature observed near "L" or 1 in Figure 4K corresponds to the location where the lowenergy (cold) magnetosheath ion reaches after penetrating into the magnetospheric side, as indicated by the plasma density (black in Figure 3C) and the red arrow in Figure 3H.
Therefore, we speculate that ion may play an important role in generating different types of waves. According to Omura et al. (1996), ions could change types of generated waves depending on the ion temperature and the ion drift by interacting resonantly with waves generated by electrons or scattering electrons. As a result, various types of waves could be generated such as ion acoustic wave, electron solitary wave, electron hole, and Langmuir wave. Nonlinear wave mode and its evolution cannot be studied by linear analysis. Further study using kinetic simulation is required for understanding how ion dynamics affect nonlinear wave mode.
This observation, however, implies that the electrostatic waves observed predominantly in the magnetospheath side of the KHV boundary may pre-heat the cold magnetosheath population that is to participate into the reconnection process moving toward an X-line via/along the inflow/separatrix region. This may explain the higher-energy (200 eV < energy <2 keV) electrons streaming toward X along the magnetosheath-side separatrix region (red arrow in Figure 3K). Large Joule dissipation during the period of the enhanced wave activity ( Figure 4G) also supports this twostep energization of the magnetosheath plasma entering into the magnetosphere via KHV-driven reconnection.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/ Supplementary Material.