Multi-scale evolution of Kelvin–Helmholtz waves at the Earth's magnetopause during southward IMF periods

At the Earth’s low-latitude magnetopause, the Kelvin–Helmholtz instability (KHI), driven by the velocity shear between the magnetosheath and magnetosphere, has been frequently observed during northward interplanetary magnetic field (IMF) periods. However, the signatures of the KHI have been much less frequently observed during southward IMF periods, and how the KHI develops under southward IMF has been less explored. Here, we performed a series of realistic 2D and 3D fully kinetic simulations of a KH wave event observed by the Magnetospheric Multiscale (MMS) mission at the dusk-flank magnetopause during southward IMF on September 23, 2017. The simulations demonstrate that the primary KHI bends the magnetopause current layer and excites the Rayleigh–Taylor instability (RTI), leading to penetration of high-density arms into the magnetospheric side. This arm penetration disturbs the structures of the vortex layer and produces intermittent and irregular variations of the surface waves which significantly reduces the observational probability of the periodic KH waves. The simulations further demonstrate that in the non-linear growth phase of the primary KHI, the lower-hybrid drift instability (LHDI) is induced near the edge of the primary vortices and contributes to an efficient plasma mixing across the magnetopause. The signatures of the large-scale surface waves by the KHI/RTI and the small-scale fluctuations by the LHDI are reasonably consistent with the MMS observations. These results indicate that the multi-scale evolution of the magnetopause KH waves and the resulting plasma transport and mixing as seen in the simulations may occur during southward IMF. VC 2022 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http:// creativecommons.org/licenses/by/4.0/). https://doi.org/10.1063/5.0067391


I. INTRODUCTION
The Kelvin-Helmholtz instability (KHI) becomes unstable when the plasma shear flow is super-Alfv enic for the magnetic field component parallel to the shear flow (Chandrasekhar, 1961).This unstable condition is easily satisfied when the magnetic field is oriented nearly perpendicular to the direction of the shear flow.At the Earth's low-latitude magnetopause, this magnetic field configuration appears when the interplanetary magnetic field (IMF) is strongly northward or southward.Indeed, clear signatures of surface waves and flow vortices, which could be generated by the KHI, have been frequently observed around the low-latitude magnetopause during periods of strong northward IMF (e.g., Sckopke et al., 1981;Kokubun et al., 1994;Slinker et al., 2003;Kivelson and Chen, 1995;Fairfield et al., 2000;Hasegawa et al., 2004;2006;Foullon et al., 2008;Kavosi and Raeder, 2015;Moore et al., 2016).These magnetopause KH waves and vortices have been believed to cause efficient mass, momentum, and energy transfer across the magnetopause and effectively contribute to forming the Earth's low-latitude boundary layer (LLBL), where plasmas of magnetosheath and magnetospheric origins are mixed, during the northward IMF periods (e.g., Nakamura, 2021 and references therein).
On the other hand, although the magnetopause boundary layer can be unstable for the KHI, even for the southward IMF in the MHD regime, as long as the magnetic field component parallel to the shear flow is weak enough (Chandrasekhar, 1961), the signatures of the magnetopause KH waves and vortices have been much less frequently observed during periods of southward IMF (e.g., Kavosi and Raeder, 2015).Hwang et al. (2011) reported a Cluster observation event of nonlinear KH vortices during a southward IMF period.In this event, observed plasma and field variations were irregular and temporally intermittent, indicating that the structure of the KH vortices was being disturbed.A recent three-dimensional (3D) fully kinetic simulation under pure southward IMF conditions demonstrated that strong evolution of magnetic reconnection across thin current layers formed within the KH waves quickly destroys the wave and vortex structures (Nakamura et al., 2020a).This reconnection-driven decay process of the KH vortex may explain the intermittent and irregular variations of the observed KH vortices as well as the low observational probability of the periodic KH waves/vortices during southward IMF.However, it was difficult to resolve the thin current layers and confirm the occurrence of the vortexinduced reconnection (VIR) in the reported Cluster event because of the insufficient time resolution in plasma measurements.
High-time-resolution fields and plasma data collected by the Magnetospheric Multiscale (MMS) mission (Burch et al., 2016) have been used to resolve small-scale physics within the KH waves and vortices, such as the VIR, as recently shown for some magnetopause KH wave events during the northward IMF (Eriksson et al., 2016a;2016b;Li et al., 2016;Vernisse et al., 2016;Stawarz et al., 2016;Tang et al., 2018;Hasegawa et al., 2020;Hwang et al., 2020;Kieokaew et al., 2020).In these northward IMF events, signatures of the VIR, such as reconnection outflow jets (Eriksson et al., 2016a;Li et al., 2016), energy dissipation within the electron diffusion region (Eriksson et al., 2016b), turbulent spectra caused by turbulent evolution of the VIR (Stawarz et al., 2016;Hasegawa et al., 2020), and the formation of multiple flux ropes and their interactions (Hwang et al., 2020;Kieokaew et al., 2020) were reported.3D fully kinetic simulations of one of these MMS events on September 8, 2015 showed that the simulated VIR signatures are reasonably consistent with many of the above observation signatures (Nakamura et al., 2017a;2017b;Nakamura et al., 2020b;Nakamura, 2021).Although past theoretical and numerical studies suggested that the plasma shear flow can reduce the rate of spontaneous reconnection and weaken the resulting solar wind transport across the magnetopause (e.g., Cassak and Otto, 2011), these realistic simulations of the observed VIR further showed that the VIR, which is a strongly driven reconnection process controlled by the super-Alfv enic vortex flow and whose rate is much higher than that of spontaneous reconnection, leads to an efficient solar wind transport.Based on the consistencies between the observations and simulations, these studies indicated that the KHI and subsequent occurrence of the VIR would indeed contribute to efficient solar wind transport during northward IMF.
Based on these previous studies, our companion paper, Blasl et al. (unpublished), reported the first MMS observations of the KH waves during southward IMF.In this event, MMS observed the intermittent and irregular variations of the surface waves, which can be interpreted as being formed by the KHI during the southward IMF and magnetosheath magnetic field.Although clear VIR signatures as reported in the above MMS observation events for the northward IMF were not found in this event, the high-time-resolution measurements of MMS frequently detected small-scale fluctuations, which can be interpreted as being generated by the lower-hybrid drift instability (LHDI), excited near the edge of the surface waves.To investigate this event in more detail, in the present paper, we perform a series of 2D and 3D fully kinetic simulations with parameters matched to this MMS event.The simulation results are consistent with the observations in terms of both large-scale surface wave signatures and smallscale LHDI fluctuations.The simulations also demonstrate that the primary KH waves induce the secondary Rayleigh-Taylor instability (RTI) at the surface bent by the KHI.The RTI forms high-density arms penetrating into the low-density magnetospheric side.This arm penetration quickly disturbs the primary KH wave structures and produces intermittent and irregular variations of the surface waves, leading to a reduction in the observational probability of the primary KH waves.Interestingly, this RTI-related reduction of the observational probability of the KH waves proceeds faster than the above VIRrelated reduction for the northward IMF, indicating that the secondary RTI may also be a key process that makes it more difficult to detect the KH waves during southward IMF.
This paper is organized as follows.Section II describes the details of the simulation model employed in this paper.Section III presents the overall simulation results focusing on the large-scale variations of the surface waves, while in Sec.IV we focus on the secondary processes induced within the primary waves.In Sec.V, we summarize the results and discuss differences in the evolution of the KH waves between the northward and southward IMF conditions.

II. MODEL A. Simulation settings
We performed a series of 2D (in the x-y plane) and 3D fully kinetic simulations that model a magnetopause crossing event observed by MMS on September 23, 2017 [Blasl et al. (unpublished)], using the high-performance particle-in-cell code VPIC (Bowers et al., 2008;2009).The x, y, and z coordinates in the simulations correspond to the direction along the velocity shear ($the magnetosheath flow in the equatorial plane), the boundary normal ($magnetosheath-tomagnetosphere), and the out-of-the-vortex-plane ($south-to-north), respectively.The initial density, magnetic field, and ion (and electron) bulk velocities across the magnetopause boundary layer are set to the values obtained from the observations near the KH vortex-like interval 15:33-15:35 UT [see Blasl et al. (unpublished) for more details of this interval].Denoting the magnetosheath and magnetospheric sides as 1 and 2, respectively, we first select the density (n 10 , n 20 ), the magnetic field (B x10 , B z10 , B x20 , B z20 ), and the bulk velocities (U x10 , U z10 , U x20 , U z20 ) on the two sides from the observations.To set the initial equilibrium, we here neglect the y-component of the initial magnetic field and velocities, which are negligibly small compared to the x and z components.We then set the initial density profiles by connecting these values across the magnetopause using a tanh(y/L 0 ) function

ARTICLE
scitation.org/journal/php(Nakamura and Daughton, 2014) as , where L 0 ¼ 2.5d i is the initial half thickness of the shear layer and d i ¼ c/x pi is the ion inertial length based on n 0 ¼ n 10 .To set the bulk velocities, particles (ions and electrons) are initialized with drifting Maxwellian velocities whose drift velocities are U x10 and U x20 for the magnetosheath and magnetospheric particles, respectively.The initial magnetic field is setup as B x,z (y) and ion flows and electron density are introduced to satisfy the Harris type variation of B x and B z .To satisfy the force balance, the temperatures for the magnetospheric ion and electron components T i,e20 are set to be higher than the magnetosheath components T i,e10 , where the ion-to-electron temperature ratio is fixed as T i0 /T e0 ¼ 5.0.
The set of values obtained from MMS data in regions 1 and 2 are n 10 /n 20 ¼ 8.0, B z10 ¼ ÀB 0 , B z20 ¼ B 0 , U x10 ¼ V 0 /2, and U x20 ¼ ÀV 0 /2, where n 10 ¼ 8 cm À3 , B 0 ¼ 12 nT, jV 0 j ¼ 290 km/s ¼ 3.0V A (V A : Alfv en speed based on n 0 and B 0 ).Note that the system is set to be in a drifting frame of reference with half the velocity of the magnetosheath flow.Since the in-plane magnetic field components (B x10 and B x20 ) are known to easily change within the vortex layer as a result of the vortex motion (e.g., Fairfield et al., 2007;Nakamura et al., 2008), it is difficult to determine B x10 and B x20 from the observations during the vortex-like fluctuating interval.Nevertheless, since the spacecraft crossed a relatively steady magnetosheath-like interval just before the vortex-like interval, we took B x10 from this interval as 0.17-0.2B0 .Regarding B x20 , since the spacecraft did not cross a clear magnetosphere-like region near the vortex-like interval, we tested various values from 0 to 0.2B 0 .In this paper, we show representative simulation runs with (B x10 , B x20 ) ¼ (0.2B 0 , 0) and (0.17B 0 , 0.17B 0 ) (i.e., a case with no in-plane field on the magnetospheric side and a case with uniform in-plane field) as listed in Table I.The total plasma beta on the magnetosheath side is b 1 ¼ 1.5, the ratio between the electron plasma frequency and the gyrofrequency is x pe /X e ¼ 2.0 and the ion-to-electron mass ratio is m i /m e ¼ 100 for all runs.The system is periodic in x (and in z for 3D), and y boundaries are modeled as perfect conductors for the fields and reflecting for the particles.
In this paper, we will show the results from four runs listed in Table I.We will first show a large-scale 2D run (run-A), in which the system size is L x Â L y ¼ 100d i Â 100d i ¼ 6144 Â 6144 cells with a total of 1.5 Â 10 10 (super)particles.In this run, we added no specific mode of initial perturbations and the KH instability grows from the random particle noises.Then, to investigate the three-dimensional effects, we will show local 3D (run-B) and 2D (run-C) runs, which feature one wavelength KH mode (m x ¼ 1) with the system size L x Â L y Â L z ¼ 15d i Â 30d i Â 10.4d i ¼ 864 Â 1728 Â 600 cells with a total of 3.6 Â 10 11 particles for 3D, and L x Â L y ¼ 15d i Â 30d i ¼ 864 Â 1728 cells with a total of 6.0 Â 10 9 particles for 2D.In these local runs, to initiate the one wavelength KH mode, we added an initial weak flow perturbation as dU iy ¼ dU ey ¼ 0:02V 0 exp À0:5ðy= ½ L 0 Þ 2 sinð2px=L x Þ.In addition, to investigate the effects of the in-plane magnetic field, we will also show a local 2D run (run-D) with the same setting as run-C except that the in-plane field is initially set to be zero.Note that for the 3D run (run-B), which used $2 Â 10 4 cores of MareNostrum at Barcelona Supercomputing Center (BSC) for more than 10 2 h, we setup the largest system size within the computer resource limitation to reproduce the MHD-scale (>d i ) primary KH waves without being affected by non-MHD effects in their initial growth phase, although the size is still about ten times smaller than the estimated wavelength (k KH $ 10 4-5 km $10 2-3 d i ) of the KH waves in the MMS observations [Blasl et al. (unpublished)].

B. Initial equilibrium
In the present simulations, the pressure balance in the boundary normal (y) direction is satisfied in the initial conditions.However, when the magnetic field strength within the current layer is significantly weaker than that in the outside regions like the Harris-type current sheet as employed in the present simulations, the fluid-type equilibrium across the layer is easily disturbed by kinetic effects.One of the main causes for this disturbance is the gyro-motion of ions located near the center of the layer.The orbits of these ions are larger due to the smaller magnetic field within the layer.As the simulation proceeds, the simulated particles start their gyro-motions and the orbits of the ions that initially located near the center but on one side of the layer can enter the other side of the layer because of the finite gyro-radius.When there is initially a large density jump across the layer as in the present simulations, the net number of the crossing ions that originally located on the two sides is not balanced, leading to the disturbance of the equilibrium.Figure 1 shows the time evolutions of the 1D profiles in the y direction of the ion density, the out-of-plane magnetic field component, and the shearing flow component for the large-scale 2D run (run-A).As the simulation proceeds, the highdensity plasma penetrates into the lower-density side as predicted.After a few ion gyro-periods, a new quasi-equilibrium state, in which the profiles are no longer largely changed, is accomplished.Note that the KH instability starts to grow after this new equilibrium is accomplished, as will be shown in Fig. 2 and Sec.III. , we see that 6 KH waves are growing near the center of the boundary layer [Fig.2(a)].After that, a part of the surface waves, which convects toward the low-density side, grows as highdensity arms penetrating into the low-density side [Figs.2(c) and 2(d)].Notice that although the arms start to roll-up a small amount, they mainly continue to grow more straight into the low-density side beyond the vortex layer and greatly disturb the layer structure.Notice also that the width of the high-density arms tends to become smaller over time as they grow in the boundary normal direction.Figure 2(e) shows the growth of the 1D power spectra (k x ) of each U iy modes.We show that in the early phase until t $ 50-60 X i À1 , the m x ¼ 6 mode clearly dominates, but after that the relative power of the other modes relatively gets larger corresponding to the evolution of the highdensity arms, and it becomes more difficult to identify one dominant mode.As will be described in detail in Sec.IV A, this additional evolution of the high-density arms is caused by the Rayleigh-Taylor instability (RTI) driven by the centrifugal force at the rippled surface.
In addition, as the large-scale surface waves develop, smallerscale fluctuations grow, especially near the edges of the waves where the density gradient is enhanced [see, for example, a white arrow in Fig. 2(a)].These fluctuations are seen mainly on the low-density side of the gradient layer and penetrate deep into the lower-density region.
As a result of the fluctuations, the gradient layer, where plasmas between the two sides across the layer are mixed, is substantially broadened near the edges of the RTI arms.As will be described in detail in Sec.IV B, these fluctuations result from the lower-hybrid drift instability (LHDI) driven by the thin density gradient layer formed at the edge of the primary surface waves/vortices.

B. Comparison with the MMS observations
To confirm the realism of the simulation, we compare the simulation results with the MMS observations of the large-scale evolution  In the virtual observations, we assume that the spacecraft crossed the vortex layer in the direction opposite to the propagation direction of the KH waves (i.e., Àx direction) at a fixed time around the transition between the linear and non-linear growth phase of the KHI (t ¼ 60 X i

À1
). Orbit-1 crossed near the magnetosheath side edge of the bulges of the primary KH waves and vortices, while the orbit-2 crossed near the center part of the waves and vortices.For orbit-2 [Figs. 3(e)-3(h)  corresponds to $1/5k KH , which is roughly consistent with the distance between the orbits-1 and 2 in the boundary normal direction (3d i $ 1/5k KH ).We have confirmed that similar consistencies between the simulation and observations are seen in a range t $ 60-65 X i À1 (i.e., near the later linear or early non-linear growth phase of the KHI).All of these consistencies indicate that the KH waves near the later linear or early non-linear growth phase as seen in the present simulation at t ¼ 60 X i À1 likely occur at the duskflank magnetopause during this interval.Note that we see the lowdensity and negative B L values more frequently and for a longer time after the crossing of the green vertical line in Figs.3(i)-3(l), indicating that the spacecraft moved deeper into the magnetosphere.See Blasl et al. (unpublished) for more details of the comparisons between the present simulation and MMS observations.

IV. SECONDARY INSTABILITIES A. Secondary RTI
To investigate the evolution process of the KHI and the subsequent secondary processes in more detail, we performed local 2D and 3D runs featuring one wavelength KH mode with a similar wavelength (15d i ) to the fastest growing mode seen in run-A (listed as runs-B and C in Table I, respectively).In these runs, the in-plane field is initially set to be uniform.Figures 4(a)-4(d) show the time evolution of the ion density contours for the 2D run (run-C).As also seen in Fig. 2, for run-A, the high-density arm penetrates deep into the low-density side, and the arm becomes narrower as the head of the arm propagates in the y direction.The growth rate of the arm, which is roughly reflected by m x ¼ 2 or 3 at least at the beginning of the arm growth as will be shown in Fig. 5 and the next paragraph, is larger than the growth rate of the primary KHI with m x ¼ 1 [compare the slopes of the curves in Fig. 4(e)].After the rapid growth of the arm [see t > 5a À1 in Fig. 4(e)], the amplitude of the parent KH mode (m x ¼ 1) is substantially reduced down to the level comparable to the m x ¼ 2 þ 3 modes (the sum of the amplitudes of m x ¼ 2 and 3 modes).
Figures 5(a)-5(e) feature the evolution of the high-density arm during t ¼ 3.2a À1 to 5.6a À1 , where a À1 ¼ k KH /V 0 is the time unit for the growth phase of the KHI (Nakamura et al., 2013).As mentioned above, the width of the arm becomes smaller as the arm grows in the boundary normal direction.As highlighted by the yellow curves in Figs.5(a) and 5(b), this corresponds to the decrease in the curvature radius of the head of the high-density arm.The ratio of the radius between t ¼ 3.2a À1 and 3.8a À1 is about 1.7 À1 , and this number is consistent with the inverse of the square of the ratio of the growth rate of the m x ¼ 2 þ 3 mode between t ¼ 3.2a À1 to 3.8a À1 [see Fig. 5(f)].Here, the m x ¼ 2 þ 3 mode roughly corresponds to the width of the arm in the x direction (i.e., the width of the arm is roughly 0.3-0.5 times the wavelength of the primary KH mode) during this time interval.Within this curved boundary layer, the current, which is produced by the anti-parallel south-to-north magnetic field across the layer and is mainly carried by ions, flows in the direction along the layer [see arrows in Figs.5(a)-5(e)].Given that the strength of the current near the density surface stays almost the same between t ¼ 3.2a À1 to 3.8a À1 , the relation between the curvature radius and growth rate of the head is consistent with the expected growth rate of the RTI excited by the centrifugal force from the curved flow (Chandrasekhar, 1961;Nakamura and Daughton, 2014); c RT / U r 2 =r À Á 0:5 , where U r is the rotating flow speed and r is the curvature radius of the flow.After the saturation of the m x ¼ 2 þ 3 mode, the mushroom-like structure forms near the head of the arm [see white arrows in Figs.5(e) as well as 4(c)], supporting that the RTI is secondarily induced at the density surface bent by the primary KHI, and forms the high-density arm penetrating into the low-density side.It should be emphasized here that since the current flows in the direction largely tilted from the KHI plane ($ the equatorial plane) for the northward IMF case, this RTI physics would be unique for the southward IMF case in which the current flows nearly along the background shear flow.

B. 3D effects
Figures 6(a)-6(f) show the time evolution of the ion density and the x-component of the electric field in the electron frame x for the local 3D run (run-B) in the x-y plane at z ¼ 0. Similar to the 2D case in the x-y plane (see Fig. 4), the highdensity arm penetrates into the low-density side, and the size of the arm head decreases as the head moves deeper into the low-density side.However, as shown in Fig. 6(g), in the 3D case, the high-density arm grows in the direction somewhat tilted from the x-y plane (vortex plane), forming a large-scale wavy structure of the arm head in the zdirection (i.e., finite k z is produced).This tilt angle corresponds to the direction perpendicular to the magnetic field measured at the arm head, as shown in Fig. 6(h).These features are consistent with the 3D evolution of the RTI in the direction satisfying k Á B $ 0 in which the growth rate of the RTI is maximized (Chandrasekhar, 1961).Note that in Fig. 6(h) we also see a weak enhancement of the wave power at around k x $ 10-20, which reflects the small-scale fluctuations generated by the LHDI as will be explained in the next paragraph.
In the 3D case, it is also notable that small-scale fluctuations of the electric field, which are dominantly seen in the perpendicular components, grow mainly on the low-density side of the density gradient layer, as seen in Figs.6(d 7(c), the most strongly growing mode of the small-scale fluctuations at t ¼ 4a À1 is m x $ 18-23 corresponding to the wavelength $6.5-8 d e in the x direction.Figure 7(d) shows the zoomed-in-view in the x-y plane of E x 0 in the region marked in Fig. 6(d), where the fluctuations occur most strongly.The positive-negative pattern of the fluctuations whose projected wavelength in the x direction is around 7d e is seen as expected from the power spectrum.The tilt angle of the edge layer in the x-y plane is about h edge $ 35 , indicating that the actual wavelength of the dominant mode is around k $ 8.0-9.5 d e , which is in a range of the hybrid-kinetic scale based on the local ion and electron temperatures and the magnetic field strength [k ?ðq i q e Þ 1=2 $ 1:4-1:7].The power spectrum [Fig.7(c)] also shows that the enhanced power of the fluctuations spreads nearly up to the electron kinetic scale (k ?q e $ 1 corresponding to m x $ 50-55).These features (the perpendicular electric field fluctuations consisting of modes from the hybrid-kinetic to electron kinetic scales) are consistent with the ones predicted from the past linear analyses for the LHDI (Daughton, 2003).In the later time [see Fig. 7(b)], the in-plane magnetic field is more strongly compressed near the edge of the vortex/ high-density arm, leading to h perp being larger and the strong power  seen more widely within the range k z /k x < tan(h perp ) max .Notice that as the peaks are scattered, the wavelength of the most strongly growing mode becomes somewhat (1.5-2 times) longer [see the blue curve in Fig. 7(c)].

Physics of Plasmas
To understand the onset condition of these small-scale modes more quantitatively, a linear dispersion relation solver for fully kinetic plasma with a perpendicular drift, which has recently been developed in Umeda and Nakamura (2018), is applied by inputting parameters obtained from the simulation.In this solver, the dispersion relation is derived in the electromagnetic and fully kinetic regime for ions and electrons that drift in the direction perpendicular to the magnetic field, which can include ion and electron cyclotron resonances unlike the conventional approximation by Davidson et al. (1977) [see Umeda and Nakamura (2018) for more details about this linear dispersion relation].This solver requires m i /m e and local values of the ion-to-electron temperature ratio T i /T e , x pe /x ce , the ion and electron drift velocities relative to the thermal speeds V d /V t , and the ratio between the thermal speeds and the speed of light V t /c for ions and electrons as input parameters to obtain the dispersion relation.Figure 7(f) shows the results of the linear analysis in which the input parameters are obtained as local values at the location marked in Fig. 4(b).Note that the parameters are obtained at the location where the small-scale fluctuations are most strongly seen in the 3D run (run-B) but from the corresponding 2D run (run-C) in which the fluctuations are suppressed and thus we can obtain stable background values.The obtained parameters are m i /m e ¼ 100, T i /T e ¼ 3.7, x pe /x ce ¼ 1.8, V di /V ti ¼ 0, V de /V te ¼ 0.16, V ti /c ¼ 0.05, and V te /c ¼ 0.26 (b $ 2.1).The dispersion relation demonstrates that a quasi-perpendicular wave mode with a wave normal angle of / $ 87 relative to the ambient magnetic field has the maximum growth rate c $ 0.06x LHR at a wave number m LH ¼ k ?ðq i q e Þ 1=2 $ 1:6, where x (right) from t ¼ 2 to 8a À1 (the same time interval as Fig. 4) for the local 3D run (run-B) at z ¼ 0. Arrows show the in-plane components of the current density.(g) 3D view of the electron density surface of n e ¼ 0.6n 0 at t ¼ 6a À1 for run-B with the electron density contours in the x-y plane at z ¼ 0, x-z plane at y ¼ 0, and y-z plane at x ¼ 0 and randomly distributed magnetic field lines.(h) 2D power spectra (k x , k z ) of U iy at t ¼ 6a À1 around the center of the boundary (y ¼ 0 6 5d i ) for run-B.The dotted line indicates the angle of the perpendicular direction to the magnetic field in the x-z plane h perp at the head of the high-density arm.x LHR is the lower-hybrid frequency.Considering the tilt angle of the gradient layer in the x-y plane h edge $ 35 as discussed above, the dominant mode seen in the simulation is in excellent agreement with the fastest growing mode obtained from the linear analysis as seen in the red curve in Fig. 7(c).This clearly indicates that the E ?fluctuations seen in the simulation are caused by the perpendicular plasma drifts and the related growth of the LHDI.Notice that the dispersion relation for a mode with / $ 89 in Fig. 7(f) shows a significantly low frequency below the ion cyclotron frequency (<0.1xLHR ), which indicates that some modes could be largely affected by the ion gyro-motion in the present simulation.This could be because when the ion-to-electron mass ratio is low, as in the present simulation (m i /m e ¼ 100), the gyro-radius of a large portion of ions can be comparable to or smaller than the wavelength of the electric field fluctuations (i.e., hybrid-kinetic scale), and hence their motions can largely affect the evolution of the waves in addition to electron motions.For example, in the present case, the ratio between the ion gyro-radius and the wavelength of the fluctuations

Physics of Plasmas
$ 0.7 m LH .This number becomes larger as the mass ratio increases and would greatly exceed one (i.e., the ion kinetic effects would be negligible) in a realistic (m i /m e ¼ 1836) case.Indeed, the linear analyses for larger mass ratio cases with m i /m e ¼ 256 [Fig.7(g)] and m i /m e ¼ 1836 [Fig.7(h)] show that the most unstable modes would be in a range of the pure LHDI without being significantly affected by the ion gyro-motion (i.e., a IC ) 1).Here, the input parameters to the linear dispersion relation solver for these high m i /m e cases are set to be the same as the case for Fig. 7(f) except for higher m i /m e and corresponding larger x pe /x ce and V te /c.An additional simulation for m i /m e ¼ 256, in which the initial settings are the same as run-B except for smaller L z (¼ 2d i ) and 1.6 [¼ (256/100) 0.5 ] times smaller grid-spacing, shows that the dominant mode for m i /m e ¼ 256 is indeed m LH $ 3 [Fig.7 b) indicate the angle of the perpendicular direction to the magnetic field in the x-z plane h perp at the head of the high-density arm, the maximum and mean h perp on the low-density side of the density gradient layer defined as the region with 0.2 < n e /n 0 < 0.4, and the wavelengths for k ?q e ¼ 1 and k ?q i q e ð Þ of the field fluctuations at the trailing edge of the KH waves, shown in Fig. 10 in Blasl et al. (unpublished) ($15:56 UT), are m i /m e ¼ 1836, T i /T e ¼ 11.6, x pe /x ce ¼ 21.4,V di /V ti ¼ 0, V de /V te ¼ 0.12, V ti /c ¼ 0.0012, and V te /c ¼ 0.015 (b $ 2.5). Figure 7(i) shows the results of the linear analysis for this MMS event.Similar to the above simulation cases, a quasi-perpendicular wave mode with / $ 87 relative to the ambient magnetic field has the maximum growth rate c $ 0.2x LHR at a wave number k ?ðq i q e Þ 1=2 $ 3.This wave number corresponds to k $ 56 km.Given the simulation result that the wavelength of the dominant modes tends to become 1.5-2 times longer as the LHDI develops [see blue curve in Fig. 7(c)], this predicted wavelength of the fastest growing mode is in reasonable agreement with the estimation from the observed field fluctuations in this MMS event ($80-100 km) [see Blasl et al. (unpublished) for more details of the MMS observations of the small-scale field fluctuations between the KH waves].Together with the above quantitative consistencies between the simulation and the MMS event on the large-scale evolution of the primary KHI/RTI shown in Sec.III B, these results naturally suggest that the LHDIdriven turbulence at the edges of the KHI/RTI waves/vortices as seen in the present simulation may occur in the Earth's magnetopause during the southward IMF.

C. Turbulent reconnection
The top panels in Fig. 8 show the B z component in the x-y plane at t ¼ 5.4a À1 (early non-linear growth phase of the KHI) in the 3D case (run-B), the corresponding 2D case (run-C), and the 2D case without the initial in-plane magnetic field (run-D).The initial

ARTICLE
scitation.org/journal/phpanti-parallel magnetic field forms the thin primary current layer (corresponding to the density gradient layer) along the edge of the largescale surface wave in all cases.However, the small-scale fluctuations induced by the LHDI disturb and broaden the upper positive-B z side (low-density side) of the layer for runs-B and D. Note that for run-C, the in-plane magnetic field, which is locally compressed and enhanced along the edge of the surface wave [see black curves in Fig. 8(d)], prevents the small-scale fluctuations from growing (i.e., the quasiperpendicular unstable LHDI modes do not exist in the 2D simulation plane), while for run-B, a similar enhancement of the in-plane field is seen but the LHDI can grow obliquely in the direction nearly perpendicular to the magnetic field as shown in Sec.IV C. The middle panels in Fig. 8 show the J z component for the three runs, while the bottom panels show the time evolution of the global reconnection rate computed from the integrated magnetic flux that crosses the (mixing) surfaces of the layer on the low-density and positive-B z (top) and the high-density and negative-B z (bottom) sides (Daughton et al., 2014).
For runs-B and C, the compressed in-plane field produces strong outof-plane currents (J z ) mainly in the low-density side of the primary gradient layer [Figs. 8(b) and 8(e)].For run-B, the magnetic flux rapidly flows into the low-density side of the primary layer [Fig.8(c)].
These results indicate that magnetic reconnection occurs within the low-density side of the layer where the strong J z (i.e., the large magnetic shear) is produced by the compressed in-plane field as mentioned above.Note that since the surfaces of the z-component of the vector potential do not exactly correspond to the in-plane field lines in 3D, the J z variations for run-B [Fig.8 .This is because the unstable mode of the tearing instability satisfying k Á B ¼ 0 cannot be found in the 2D simulation planes (i.e., there is no anti-parallel, in-plane component of the magnetic field).
It is notable that the rate of the flux inflow on the high-density side for run-B [see the red curve in Fig. 8(c)] is about ten times smaller than that on the low-density side (R $ 0.5).This indicates that in this 3D case, reconnection occurs mostly on the low-density side and hardly occurs between the two sides across the primary layer (i.e., between the northward and southward magnetic field lines across the magnetopause).This could be because the small-scale fluctuations broaden the low-density side of the layer and prevent the inflowing flux on the low-density side from reaching the high-density side across the layer.This result is not seen in recent 3D fully kinetic simulations in the southward IMF case (Nakamura et al., 2020a), in which the initial density jump is set to be weak (n 10 /n 20 ¼ 3.0), no significant evolution of the LHDI fluctuations is seen, and reconnection strongly develops across the layer.Note that a weak enhancement of the flux inflow rate on the high-density side is seen not only for run-B but also for run-D in which reconnection cannot occur in the 2D simulation plane [see the red curves in Figs.8(c) and 8(i)].These weak enhancements on the high-density side result from the mixing caused by the LHDI fluctuations.

D. Plasma transport and mixing
Figure 9(a) shows the time evolution of the thickness in the boundary normal (y) direction of the region where the profile in the x-direction contains the density variation larger than 0.3(n 10 À n 20 ), for runs B-D.This thickness roughly corresponds to the penetration distance of the high-density arm into the low-density side.For all runs, we see a similar continuous penetration of the arm, corresponding to a similar evolution of the large-scale KHI/RTI (see also top panels in Fig. 8), at a speed 0.3-0.4V0 until t $ 6a À1 .Note that for runs B-D, the penetration is suppressed at t $ 6a À1 by the effects of the simulation boundaries in the y-direction, while for run-A, in which the system size is set to be larger than that for runs B-D, the penetration continues even after t $ 6a À1 without being affected by the boundaries [see the black curve in Fig. 9(a)].
Figure 9(b) shows the time evolution of the averaged thickness in the y direction of the mixing region defined as the region with jF e j < 0.99 for runs B-D.We clearly see that for runs-B and D, in which the LHDI turbulence occurs, the mixing region expands much more efficiently than run-C.This indicates that the LHDI turbulence dominantly contributes to the mixing between the high-and low-density sides in the vortex layer, while the RTI and the resulting penetration of the high-density arm, which contribute to the convective transport of high-density plasmas into the low-density side as shown in Fig. 9(a), do not significantly contribute to the plasma mixing between the two sides.Similar evolutions of the high-density arm for all runs shown in Fig. 9(a) indicate that the local mixing by the LHDI turbulence does not significantly disturb the large-scale evolution of the KHI and RTI-i.e., the primary KHI and the subsequent secondary RTI are the dominant processes to control the large-scale evolution of the shear layer.As described in Sec.IV A, during this large-scale evolution of the layer, the high-density arm becomes narrower as the arm develops in the boundary normal direction.This suggests that when the spacecraft crosses the vortex layer, the observed dominant mode would be different at different locations (relative to the high-density arm) and/or at different times (depending on the growth phases of the KHI and RTI).Figures 10(a) and 10(b) show k x spectra for run-B of the total pressure for each y (corresponding to the spectra from the virtual observations in which the virtual probe crosses the layer in the x direction at each y) at t ¼ 3a À1 and t ¼ 6a À1 (before and after the onset of the RTI).At t ¼ 3a À1 , a dominant mode of the KHI (m x ¼ 1) is clearly seen near y $ 0, while at t ¼ 6a À1 , smaller-scale modes become visibly stronger especially at the low-density side (y > 0) and it is no longer easy to identify one dominant mode.Notice that this flattening of the spectra corresponds to the formation of intermittent and irregular variation patterns of the surface waves as shown in Figs.3(a)-3(h) for run-A.

Physics of Plasmas
Figures 10(c) and 10(d) show the k x spectra integrated over y of the total pressure and the y component of the ion bulk velocity.At t ¼ 3a À1 , the amplitude of the primary KH (m x ¼ 1) mode is more than one order of magnitude larger than that of the smaller-scale modes, while at t ¼ 6a À1 and later times, the difference in the amplitude between the KHI and smaller-scale modes becomes much smaller.Figure 10(f) shows the time evolution of the ratio of the amplitude of the primary KH (m x ¼ 1) mode and the smaller-scale FIG. 10. [(a) and (b)] 1D power spectra (k x ) for run-B of the total pressure P t for each y averaged in z at t ¼ 3a À1 (linear phase) and t ¼ 6a À1 (early non-linear phase).[(c) and (d)] The same spectra as (a) and (b) but integrated in y of (c) P t and (d) U iy at t ¼ 3, 6, and 9a À1 .[(e) and (f)] Time evolution for run-B of (e) the thickness of the region with the large-scale density fluctuations [the same as the red curve in Fig. 9(a)], and (f) the ratio of the P t power between m x ¼ 1 and m x ¼ 2-8 within the region having the moderate density fluctuations.The black dashed curves in (e) and (f) show the results calculated from the 3D run of the September 8, 2015 event employed in Nakamura (2021), in which the KH waves were observed during northward IMF.(m x ¼ 2-8) modes.The ratio becomes below 10 (i.e., the amplitudes of the primary KH and smaller scale modes become within an order of magnitude) at t $ 3.5-4a À1 , indicating that it quickly becomes difficult to identify the primary KH mode just after the onset of the secondary RTI (t $ 3a À1 ).As shown in Fig. 10(e) [the same as the red curve in Fig. 9(a)], the decrease in the ratio roughly corresponds to the evolution of the high-density arm.These results may explain the low observational probability of the KH waves during southward IMF periods (e.g., Kavosi and Raeder, 2015)-i.e., during southward IMF, even when the KHI is unstable at the magnetopause, the evolution of the secondary RTI and the resulting penetration of the high-density arm could make it difficult to identify the periodic KH waves.

V. SUMMARY AND DISCUSSION
A. Summary Based on the 2D and 3D fully kinetic simulations of the MMS event on September 23, 2017, in which signatures of the KH waves were observed during a southward IMF interval as shown in our companion paper Blasl et al. (unpublished), we investigated the evolution process of the KH waves and related secondary waves/instabilities during the southward IMF.The main results of this paper are summarized as follows: 1.The KHI is indeed unstable and forms surface waves at the magnetopause boundary layer (velocity shear layer) under the parameters obtained from the MMS observations during southward IMF on September 23, 2017.2. As the primary KHI grows, the ion current, which is produced by the anti-parallel (south-north) magnetic field across the surface and flows in the direction nearly parallel to the background shear flow (i.e., in the equatorial plane), is bent in the boundary normal direction and drives the secondary RTI. 3. The growth of the RTI results in the penetration of high-density arms into the low-density (magnetospheric) side, which disturbs the vortex motion and significantly reduces the observational probability of the primary KH waves.4. At the edge of the KH waves and the high-density arms, where the thin density gradient layer forms, the LHDI is induced and produces the small-scale (hybrid kinetic-scale) fluctuations in the perpendicular component of the electric field.The LHDI grows mainly in the low-density side of the gradient layer, leading to the diffusion of the layer.5.The secondary RTI and the resulting penetration of the highdensity arms cause plasma transport into the low-density side, while the LHDI contributes to the plasma mixing along the edge of the primary KH waves and the high-density arms.6.Although signatures of magnetic reconnection (inflowing magnetic flux) are observed within the LHDI turbulence, reconnection does not significantly grow across the primary shear layer.7. The large-scale variations of the surface waves by the KHI and RTI and the small-scale fluctuations by the LHDI are reasonably consistent with those seen in the MMS observations [see Blasl et al. (unpublished) for details of the MMS observations].

B. IMF dependence
A key difference in the profile across the low-latitude magnetopause between the northward and southward IMF cases is the direction (and the amplitude) of the current produced by the magnetic shear between the magnetosheath and magnetospheric sides compared to the background shear flow (the equatorial plane).During the northward IMF, the magnetopause current is weaker and flows in the direction largely tilted from the shear flow, while during southward IMF, stronger current flows in the direction closer to the shear flow.As shown in Sec.IV A, during southward IMF, this strong in-plane current causes the secondary RTI at the boundary layer bent by the KH waves, leading to the deep penetration of the high-density arm into the magnetospheric side.Here, the ratio of the growth rates of the primary KHI and the secondary RTI can roughly be estimated in the incompressible MHD regime (Chandrasekhar, 1961) as where k KH , k RT , k KH , and k RT are the wavenumber and wavelength of the KHI and RTI, respectively, r is the curvature radius of the surface wave, and V MP is the parallel (to the shear flow) velocity component of ions that partly carry the magnetopause current.Assuming that the current is carried only by ions with typical parameters near the flank magnetopause as the layer thickness L MP $ 10 3 km, the magnetic field varying from À20 to 20 nT across the layer (i.e., DB z $ 40 nT) and the density near the center of the layer n MP $ 1 cm À3 , V MP is roughly estimated as V MP $ (DB z /L MP )/(l 0 en MP ) $ 200 km/s, which is comparable to the amplitude of the shear flow near the flank magnetopause (i.e., V MP $ V 0 ).Since k RT and r, both of which depend on k KH , are only a few times smaller than k KH and comparable to k KH , respectively, as shown in Sec.IV A, the growth rate of the secondary RTI would commonly be comparable to that of the primary KHI at the flank magnetopause.Thus, the strong growth of the RTI and the resulting deep penetration of the high-density arm as seen in the present simulations could actually occur at the flank magnetopause during southward IMF.The black dashed curves in Figs.10(e) and 10(f) show the results from the 3D run of an MMS observation event of KH waves during northward IMF on September 8, 2015 (Nakamura et al., 2017a;2017b;Nakamura, 2021).In this northward IMF case, the vortex-induced reconnection quickly destroys the structure of the primary KH vortex (Nakamura et al., 2017b), leading to the quick decrease in the relative amplitude of the primary KH mode as seen in Fig. 10(f).Interestingly, the decrease rate in the present simulation for southward IMF is somewhat larger than the simulation for northward IMF [compare red solid and black dashed curves in Fig. 10(f)].This could be because of the quicker evolution of the high-density arm due to the RTI for the southward IMF case [compare red solid and black dashed curves in Fig. 10(e)].This result indicates that it would be rather difficult to identify the KH waves from the periodicity of the surface waves during southward IMF even when the unstable condition for the KHI is satisfied.In addition, notice that the evolution of the high-density arm for the southward IMF leads to intermittent and irregular variation patterns of the plasma and field parameters as shown in Figs.3(a)-3(h).Similar variation patterns were observed in the MMS event treated in this paper [see Figs.3(i)-3(l)] as well as the past observation of the KH waves during southward IMF by Cluster (Hwang et al., 2011).Combined with similar estimated growth rates of the RTI and KHI, the intermittent and irregular variation patterns of the surface waves would be a common feature of KH waves observed during southward IMF. C. Some remarks on future work

Physics of Plasmas
Recent 3D kinetic simulations of the magnetopause KHI showed that reconnection induced by the magnetic shear across the magnetopause quickly decays the non-linear vortex structure for both northward (Nakamura et al., 2017a;2017b) and southward (Nakamura et al., 2020a) IMF cases.On the other hand, in the present simulation which has a density jump across the magnetopause that is more than two times higher than the jumps employed in these recent simulations, the secondary RTI and LHDI driven by the large density gradient are the main processes that control the non-linear vortex structure, and reconnection occurs only locally within the LHDI turbulence.These results suggest that while reconnection plays a key role in controlling the non-linear KH vortex structure for both northward and southward IMFs if the density jump across the magnetopause is sufficiently small, as the density jump becomes larger, the secondary RTI and LHDI become more active and can play a more significant role in controlling the vortex structure.We expect the density jump to be generally larger for southward IMF conditions when the LLBL is thinner or less prominent than for northward IMF (Mitchell et al., 1987).An additional survey on the density ratio across the magnetopause would lead to a more systematic understanding of the secondary reconnection, RTI and LHDI processes and their dependence on the density ratio across the magnetopause.
In the present 3D simulation, as shown in Sec.IV A, the secondary RTI dominates plasma transport into the low-density side and strongly disturbs the large-scale vortex structure, while as shown in Sec.IV B, the secondary LHDI broadens the low-density side of the density gradient layer where plasmas are mixed across the layer.Assuming that the primary KHI, whose wavelength and phase speed are k KH $ 10 4-5 km and V ph $ V 0 $ 200-300 km/s, respectively, starts growing near the dayside-to-flank magnetopause, since the secondary processes are driven in the non-linear growth phase of the primary KHI, which corresponds to Dx $ V 0 Á 4a À1 $ 4k KH $ a few to 50 Earth radii from the onset location of the primary KHI, the vortexinduced diffusive plasma transport and mixing would be active at the flank-to-tail magnetopause during southward IMF.
However, to more quantitatively understand these diffusive processes, further simulations under more realistic conditions would be required.For example, as shown in Fig. 7, since the growth rate of the LHDI and the related ion-kinetic effects depend on the ion-to-electron mass ratio, to understand the actual roles of the LHDI turbulence (as well as reconnection within the turbulence), larger-scale simulations with higher mass ratios would need to be performed.In addition, although the evolution of the secondary RTI would not be significantly affected by the mass ratio, as shown in Fig. 10(a), the depth that the high-density arm can penetrate into the low-density side depends on the simulation system size.Thus, to estimate the actual transport rate and understand the role of the RTI, larger-scale simulations would also need to be performed.Furthermore, since the pure (i.e., non-KHI-related) reconnection process and related phenomena, such as flux transfer events frequently occur at the low-latitude magnetopause during southward IMF (e.g., Fuselier, 2021 and references therein), to more comprehensively understand the local physics of the KHI, it would also be required to consider the global coupling with the other processes.Although it is difficult to include these global physics in the fully kinetic regime, global hybrid (fluid-electron and kinetic-ions) simulations, which recently became applicable to the whole Earth's magnetosphere (e.g., Palmroth et al., 2018;Guo et al., 2021), may provide important support for treating such global inter-process couplings.
As discussed above, the estimated ion velocity component produced by the magnetopause current V MP ($200 km/s) is close to the background velocity shear across the magnetopause (V MP $ V 0 ).Although V MP is basically dominant within the current layer and the shearing flow component is less effective within the current layer, the shearing flow can modify the total flow profile to some degree.At the dusk-side magnetopause as in the present MMS event, the direction of the shearing flow is close to that of the ion current, and therefore the shearing flow can enhance the total flow and strengthen the growth of the RTI, while on the dawn-side, the direction of the shearing flow is mostly opposite from the ion current and therefore the RTI can be weakened.This effect would increase with increasing down tail distance, as the velocity shear becomes stronger.Thus, there may be a local-time difference in the plasma transport rate by the KHI/RTI during southward IMF.Systematically surveying the local-time dependence would also be important to more practically understand the effects of the magnetopause KHI during southward IMF.
Finally, it should be noted that the present simulations started from a Harris-type current layer in which the total pressure within the layer is sustained by the plasma pressure and the magnetic field strength within the layer is weaker than that in the background regions.Although no significant decrease in the magnetic field strength near the current layer crossings of this MMS event [for example, see Fig. 3 in Blasl et al. (unpublished)] may indicate that the observed KH waves were driven at a non-Harris-type layer, it is difficult to know the exact initial equilibrium conditions from the observations.Past observational and analytical studies of the Earth's magnetopause crossings showed that in addition to the Harris-type layer, the force-free configuration of the current layer, in which the magnetic field strength is nearly constant across the layer, also exists at the magnetopause (Panov et al., 2011).Investigating these different types of the initial current layers would be required to more systematically understand the evolution process of the magnetopause KH waves during southward IMF.

VI. CONCLUSIONS
We have performed a series of 2D and 3D fully kinetic simulations of an MMS observation event on September 23, 2017 in which the KH waves were observed at the dusk-flank magnetopause during southward IMF.This is the first numerical challenge to investigate the magnetopause KHI under realistic parameters for the southward IMF.The simulations demonstrate (i) that the KHI is unstable under the parameters of this MMS event; (ii) that the growth of the secondary RTI and the subsequent penetration of high-density arms into the low-density side disturb the primary vortex structure and produce intermittent and irregular variations of the surface waves, leading to a low observational probability of the primary KH waves; and (iii) that the secondary LHDI induced near the vortex edge causes the efficient plasma mixing across the magnetopause.The large-scale variations of the surface waves and the small-scale fluctuations from the LHDI are reasonably consistent with the MMS observations.These results indicate that the multi-scale secondary processes and the resulting plasma transport and mixing as shown in the present simulations may actually occur at the flank-to-tail magnetopause during southward IMF.

FIG. 1 .
FIG. 1.Time evolution of the profiles in the y-direction at x ¼ L x /2 for the largescale 2D run (run-A) of (a) the ion density n i , (b) the out-of-plane component of the magnetic field B z , and (c) the x-component of the ion bulk velocity U ix from t ¼ 0 to 40 X i À1 .
FIG. 3. (a)-(h) Virtual observations of the KH waves for run-A, made along the orbits-1 and 2 shown in Figs.2(b) and (i)-(l) the in situ observations by the MMS1 spacecraft for 1 h from 15:30 UT of the ion density, L (¼ z) component of the magnetic field, N (¼ Ày) component of the ion bulk velocity and the total pressure.The time interval in (i)-(l) is the same as the one in Fig. 3 in Blasl et al. (unpublished).

FIG. 4 .
FIG. 4. (a)-(d) Time evolution of color contours in the x-y plane of n i from t ¼ 2 to 8a À1 for the local 2D run (run-C), where a À1 ¼ k KH /V 0 is the time unit for the growth phase of the KHI (Nakamura et al., 2013).(e) Time evolution of the 1D power spectra (k x ) of U iy modes (m x ¼ 1 and m x ¼ 2 þ 3) obtained by calculating the power of each mode at each y and then picking up the maximum value for each mode.Vertical dotted lines in (e) indicate the times shown in (a)-(d).
)-6(f).These small-scale electric field fluctuations produce corresponding density fluctuations and broaden the low-density side of the layer [compare Figs.6(a)-6(c) and 6(d)-6(f)].The amplitude of the fluctuations is locally enhanced in the thin gradient layer compressed by the non-linearly developed primary surface wave and vortex [see the right-side of the edge layer of the primary wave/arm in Figs.6(d)-6(f)].The 3D view of the density surface in Fig. 6(g) shows that the wavy patterns of the small-scale density fluctuations are nearly aligned with the local magnetic field lines, indicating that the wavevector of the fluctuations develops in the direction nearly perpendicular to the local magnetic field.This point is clearly seen in the power spectrum at t ¼ 4a À1 [Fig.7(a)].The power of the electric field fluctuations is enhanced dominantly in the direction nearly perpendicular to the mean magnetic field near the low-density side of the density gradient layer [k z /k x $ tan(h perp )], which roughly corresponds to the direction perpendicular to the local magnetic field in the region where the fluctuations are enhanced.As seen in the red curve in Fig.

FIG. 5 .
FIG. 5. (a)-(e) Time evolution of n i contours for run-C from t ¼ 3.2 to 5.6a À1 .Arrows show the in-plane components of the current density.The white and red dashed curves in (a) and (b) show the density surface with (n 10 þ n 20 )/2 and a portion of the fitted circles near the top of the primary wave structures, respectively.r in (a) and (b) indicate the curvature radius for the yellow curves.(f) Time evolution of the 1D power spectra (k x ) of U iy modes m x ¼ 1 and m x ¼ 2 þ 3. Vertical lines indicate the times shown in (a)-(e).

FIG. 6 .
FIG. 6. (a)-(f) Time evolution of color contours of n i (left) and the x-component of the electric field in the electron frame E x 0¼ (E þ U e Â B)x (right) from t ¼ 2 to 8a À1 (the same time interval as Fig.4) for the local 3D run (run-B) at z ¼ 0. Arrows show the in-plane components of the current density.(g) 3D view of the electron density surface of n e ¼ 0.6n 0 at t ¼ 6a À1 for run-B with the electron density contours in the x-y plane at z ¼ 0, x-z plane at y ¼ 0, and y-z plane at x ¼ 0 and randomly distributed magnetic field lines.(h) 2D power spectra (k x , k z ) of U iy at t ¼ 6a À1 around the center of the boundary (y ¼ 0 6 5d i ) for run-B.The dotted line indicates the angle of the perpendicular direction to the magnetic field in the x-z plane h perp at the head of the high-density arm.
(e)] as predicted from the linear analysis [Fig.7(g)], supporting the predictions from the linear analyses.Based on the above consistencies between the simulation and the linear analysis, we further applied the linear dispersion relation solver to real parameters obtained from the MMS observation event treated in this paper andBlasl et al. (unpublished).The input parameters, based on the local values near the time interval with large amplitudes

FIG. 7 .
FIG. 7. (a)-(c) 2D power spectra (k x , k z ) of E x 0 and sum of the spectra in k z at t ¼ 4 and 6a À1 around the center of the boundary (y ¼ 0 6 5d i ) for run-B.The dotted lines in (a) and (b) indicate the angle of the perpendicular direction to the magnetic field in the x-z plane h perp at the head of the high-density arm, the maximum and mean h perp on the low-density side of the density gradient layer defined as the region with 0.2 < n e /n 0 < 0.4, and the wavelengths for k ?q e ¼ 1 and k ?q i q e ð Þ 1=2 ¼ 1.(d) The zoomed-in view of the E x 0 contour in the region marked in Fig. 6(d).(e) The same as (d), but for an additional run with m i /m e ¼ 256.(f)-(i) Linear dispersion relations at the edge of the KH waves derived by Umeda and Nakamura (2018).Panels show wave modes with wave normal angles quasi-perpendicular to the ambient magnetic field (/ ¼ 70 À 89 ) for parameters, obtained from the simulation near the compressed density gradient layer at t ¼ 4a À1 marked in Fig. 4(b) (f), similar to (f) but for m i /m e ¼ 256 (g) and m i /m e ¼ 1836 (h), and obtained from the MMS observations on September 23, 2017 Blasl et al. (unpublished).
FIG. 7. (a)-(c) 2D power spectra (k x , k z ) of E x 0 and sum of the spectra in k z at t ¼ 4 and 6a À1 around the center of the boundary (y ¼ 0 6 5d i ) for run-B.The dotted lines in (a) and (b) indicate the angle of the perpendicular direction to the magnetic field in the x-z plane h perp at the head of the high-density arm, the maximum and mean h perp on the low-density side of the density gradient layer defined as the region with 0.2 < n e /n 0 < 0.4, and the wavelengths for k ?q e ¼ 1 and k ?q i q e ð Þ 1=2 ¼ 1.(d) The zoomed-in view of the E x 0 contour in the region marked in Fig. 6(d).(e) The same as (d), but for an additional run with m i /m e ¼ 256.(f)-(i) Linear dispersion relations at the edge of the KH waves derived by Umeda and Nakamura (2018).Panels show wave modes with wave normal angles quasi-perpendicular to the ambient magnetic field (/ ¼ 70 À 89 ) for parameters, obtained from the simulation near the compressed density gradient layer at t ¼ 4a À1 marked in Fig. 4(b) (f), similar to (f) but for m i /m e ¼ 256 (g) and m i /m e ¼ 1836 (h), and obtained from the MMS observations on September 23, 2017 Blasl et al. (unpublished).

FIG. 8 .
FIG. 8. (Top and middle) Color contours of B z and J z at t ¼ 5.4a À1 in the x-y plane (z ¼ 0 for the 3D run), and (bottom) time evolutions of the normalized, global reconnection rate R (Daughton et al., 2014) computed from the integrated magnetic flux that crosses the mixing surfaces, defined as jF e j ¼ 0.99, on the low-density and positive-B z (top) and high-density and negative-B z (bottom) sides for runs-B-D.Here F e ¼ (n e1 À n e2 )/(n e1 þ n e2 ) is the mixing measure.The black curves for runs-B and C show the surfaces of the z-component of the vector potential corresponding to the in-plane magnetic field lines in 2D.The vertical lines in the bottom plots indicate the time shown in the top plots.
(b)] are not exactly aligned with the surfaces of the potential [black curves in Fig. 8(b)], especially near the head of the high-density arm where the magnetic field lines are three-dimensionally twisted [Fig.6(g)].Note also that we see no significant enhancement of the inflowing flux in the 2D cases [Figs.8(f) and 8(i)]

FIG. 9 .
FIG. 9. Time evolutions of (a) the length in the boundary normal (y) direction of the region with the large-scale density fluctuations (i.e., the high-density arm), defined as the region where the profile in the x-direction contains the density variation larger than 0.3(n 10 À n 20 ), and (b) the averaged thickness in y of the mixing region with jF e j < 0.99 for runs B-D.The black curve in (a) shows the result from run-A in which the normalized units k KH and a À1 ¼ k KH /V 0 are employed for the fastest growing mode m x ¼ 6.The vertical lines in (a) and (b) indicate the time shown in Fig. 8.
ARTICLEscitation.org/journal/php E. Detection probability of the primary KH mode

TABLE I .
Key differences among simulation runs employed in this paper.