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

The Kelvin–Helmholtz instability (KHI) becomes unstable when the plasma shear flow is super-Alfvénic 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 non-linear 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 vortex-induced 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énic 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 small-scale 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 VIR-related 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.

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-to-magnetosphere), 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₁₀, n₂₀), 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₀) function (Nakamura and Daughton, 2014) as $n1y=n1021−tanhzL0$ and $n2y=n2021+tanhzL0$⁠, where L₀ = 2.5d_i is the initial half thickness of the shear layer and d_i = c/ω_pi is the ion inertial length based on n₀ = n₁₀. 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) = $Bx,z1021−tanhyL0+Bx,z2021+tanhyL0$⁠. Additional electron 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₁₀/n₂₀ = 8.0, B_z10 = −B₀, B_z20 = B₀, U_x10 = V₀/2, and U_x20 = −V₀/2, where n₁₀ = 8 cm⁻³, B₀ = 12 nT, |V₀| = 290 km/s = 3.0V_A (V_A: Alfvén speed based on n₀ and B₀). 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.2B₀. 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₀. In this paper, we show representative simulation runs with (B_x10, B_x20) = (0.2B₀, 0) and (0.17B₀, 0.17B₀) (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 β₁ = 1.5, the ratio between the electron plasma frequency and the gyrofrequency is ω_pe/Ω_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¹⁰ (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¹¹ particles for 3D, and L_x × L_y = 15d_i × 30d_i = 864 × 1728 cells with a total of 6.0 × 10⁹ particles for 2D. In these local runs, to initiate the one wavelength KH mode, we added an initial weak flow perturbation as $δUiy=δUey=0.02V0exp−0.5(y/L0)2sin(2πx/Lx)$⁠. 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⁴ cores of MareNostrum at Barcelona Supercomputing Center (BSC) for more than 10² 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 (λ_KH ∼ 10^4–5 km ∼10^2–3 d_i) of the KH waves in the MMS observations [Blasl et al. (unpublished)].

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 high-density 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.

Figures 2(a)–2(d) show the time evolution of the ion density for run-A. At t = 50 Ω_i⁻¹, 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 high-density 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 Ω_i⁻¹, 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 high-density 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, smaller-scale 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.

To confirm the realism of the simulation, we compare the simulation results with the MMS observations of the large-scale evolution of the KH waves. Figures 3(a)–3(h) show the virtual observations along the orbits-1 and 2 marked in Fig. 2(b), while Figs. 3(i)–3(l) show the MMS observations during 15:30–16:30 UT on September 23, 2017. 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 Ω_i⁻¹). 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)], we clearly see regular patterns of the KH waves in which the density enhancement aligns with the positive to negative B_L (= B_z) oscillation, and the minimum U_N and maximum P_t are seen near the start point of these density and B_L variations [see the variations near the vertical line in Figs. 3(e)–3(h)]. For the orbit-1, these variations of the KH waves cannot clearly be seen, especially for the pressure [Fig. 3(d)], and instead a clear density drop aligns with positive B_L and U_N peaks [see the variations near the vertical line in Figs. 3(a)–3(d)].

Interestingly, variations similar to the virtual observations of both orbits-1 and 2 are seen in the MMS observations shown in Figs. 3(i)–3(l). During this one-hour interval, MMS repeatedly encountered positive–negative variations of B_L accompanied by density variations, indicating multiple encounters of the magnetopause surface waves or the oscillations of the magnetopause in the boundary normal directions. The observed variation patterns of the boundary normal vectors at these density changes strongly suggest that these variations were most likely caused by the surface waves [see Fig. 3 in Blasl et al. (unpublished)]. During the early part of this interval, MMS observed the high-density plasmas for a longer time, while at a later time, MMS observed the low-density plasmas for a longer time and more frequently encountered wave-like repeated variations, indicating that the spacecraft moved from the magnetosheath side toward the center of the layer where the surface waves were active. The variation patterns of the density, B_N, V_N, and P_t during the early and later intervals are quantitatively consistent with the virtual observations for the orbits-1 and 2 in the normalized units of the simulation (B₀ = 12 nT, n₀ = 8 cm³), respectively [compare the variations near the blue vertical line in Figs. 3(a)–3(d) and the green vertical line in Figs. 3(e)–3(h) and those near the blue and green vertical lines in Figs. 3(i)–3(l), respectively]. Furthermore, given the averaged V_N during 15:40–16:00 UT (the transition between the early and later intervals) ⟨V_N⟩ ∼ +9.5 km/s, and assuming that the magnetopause moved toward the magnetosheath side in the boundary normal direction at this speed, it is inferred that the spacecraft relatively moved about 10⁴ km toward the magnetospheric side during 15:40–16:00 UT. Since the estimated wavelength of the observed surface waves is about 5 × 10⁴ km, this distance corresponds to ∼1/5λ_KH, which is roughly consistent with the distance between the orbits-1 and 2 in the boundary normal direction (3d_i ∼ 1/5λ_KH). We have confirmed that similar consistencies between the simulation and observations are seen in a range t ∼ 60–65 Ω_i⁻¹ (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 Ω_i⁻¹ likely occur at the dusk-flank magnetopause during this interval. Note that we see the low-density 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.

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 > 5α⁻¹ 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.2α⁻¹ to 5.6α⁻¹, where α⁻¹ = λ_KH/V₀ 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.2α⁻¹ and 3.8α⁻¹ is about 1.7⁻¹, 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.2α⁻¹ to 3.8α⁻¹ [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.2α⁻¹ to 3.8α⁻¹, 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); $γRT∝Ur2/r0.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.

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 E_x′ = (⁠$E+Ue×B$⁠)_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 high-density 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 z-direction (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)–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 = 4α⁻¹ [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(θ_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. 7(c), the most strongly growing mode of the small-scale fluctuations at t = 4α⁻¹ 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′ 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 θ_edge ∼ 35$°$⁠, indicating that the actual wavelength of the dominant mode is around λ ∼ 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⊥(ρiρ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⊥ρ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 θ_perp being larger and the strong power seen more widely within the range k_z/k_x < tan(θ_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)].

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, ω_pe/ω_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, ω_pe/ω_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 (β ∼ 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 γ ∼ 0.06ω_LHR at a wave number $mLH=k⊥(ρiρe)1/2 ∼ 1.6$⁠, where ω_LHR is the lower-hybrid frequency. Considering the tilt angle of the gradient layer in the x–y plane θ_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.1ω_LHR), 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 $αIC=ρi/[2π(ρiρe)1/2/mLH]=(mLH/2π)· (mi/me)1/4·(Ti/Te)1/4$ ∼ 0.7 $mLH$⁠. 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., $α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 ω_pe/ω_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 $mLH$ ∼ 3 [Fig. 7(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 and Blasl et al. (unpublished). The input parameters, based on the local values near the time interval with large amplitudes 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, ω_pe/ω_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 (β ∼ 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 γ ∼ 0.2ω_LHR at a wave number $k⊥(ρiρe)1/2 ∼ 3$⁠. This wave number corresponds to λ ∼ 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 LHDI-driven 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.

The top panels in Fig. 8 show the B_z component in the x–y plane at t = 5.4α⁻¹ (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 anti-parallel magnetic field forms the thin primary current layer (corresponding to the density gradient layer) along the edge of the large-scale 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 quasi-perpendicular 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 out-of-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(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)]. 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₁₀/n₂₀ = 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.

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₁₀ − n₂₀), 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.4V₀ until t ∼ 6α⁻¹. Note that for runs B–D, the penetration is suppressed at t ∼ 6α⁻¹ 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 ∼ 6α⁻¹ 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 |F_e| < 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 = 3α⁻¹ and t = 6α⁻¹ (before and after the onset of the RTI). At t = 3α⁻¹, a dominant mode of the KHI (m_x = 1) is clearly seen near y ∼ 0, while at t = 6α⁻¹, 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.

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 = 3α⁻¹, 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 = 6α⁻¹ 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 (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–4α⁻¹, indicating that it quickly becomes difficult to identify the primary KH mode just after the onset of the secondary RTI (t ∼ 3α⁻¹). 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.

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 high-density 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].

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, λ_KH, and λ_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³ km, the magnetic field varying from −20 to 20 nT across the layer (i.e., ΔB_z ∼ 40 nT) and the density near the center of the layer n_MP ∼ 1 cm⁻³, V_MP is roughly estimated as V_MP ∼ (ΔB_z/L_MP)/(μ₀en_MP) ∼ 200 km/s, which is comparable to the amplitude of the shear flow near the flank magnetopause (i.e., V_MP ∼ V₀). Since λ_RT and r, both of which depend on λ_KH, are only a few times smaller than λ_KH and comparable to λ_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.

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 λ_KH ∼ 10^4–5km and V_ph ∼ V₀ ∼ 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 Δx ∼ V₀$·$ 4α⁻¹ ∼ 4λ_KH ∼ a few to 50 Earth radii from the onset location of the primary KHI, the vortex-induced 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₀). 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.

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.

This work was supported by the Austrian Research Fund (FWF) (No. P32175-N27). For the simulations employed in this paper, we acknowledge PRACE for awarding us access to MareNostrum at Barcelona Supercomputing Center (BSC), Spain. A part of the simulation data were analyzed with resources at the Space Research Institute of Austrian Academy of Sciences. The observation data employed in this paper were obtained from the MMS spacecraft and are publically available via NASA resources and the Science Data Center at CU/LASP (https://lasp.colorado.edu/mms/sdc/public/). The work by H.H. was supported by JSPS Grant-in-aid for Scientific Research KAKENHI JP21K03504. J.E.S. was supported by the Royal Society University Research Fellowship (No. URF\R1\201286). The work by T.U. was supported by JSPS Grant-in-aid for Scientific Research KAKENHI JP19H01868. Y.-H.L. and S.A.P. are grateful for support from grants NSF-DoE 1902867 and NASA MMS 80NSSC18K0289.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The data that support the findings of this study are openly available in Zenodo at https://doi.org/10.5281/zenodo.5725490 [Nakamura (2021)].