In this study, we present the first observations from the Magnetospheric Multiscale (MMS) mission of the Kelvin–Helmholtz instability (KHI) at the dusk-flank magnetopause during southward interplanetary magnetic field conditions on September 23, 2017. The instability criterion for the KHI was fulfilled for the plasma parameters observed throughout the event. An analysis of the boundary normal vectors based on the application of the timing method onto the magnetic field and the electron density data and the minimum directional derivative method onto the magnetic field data shows signatures of surface waves in the plane parallel to the velocity shear. A comparison to 2D fully kinetic simulations demonstrates reasonable consistencies with the formation of surface waves generated by the KH instability, as well as the structures of rolled-up KH waves. The observations further indicated low density faster than sheath plasma as an indicator of rolled-up vortices, which is also consistent with the simulations. All of these results show that the observed waves and vortices are most likely generated by the KH instability. High-time resolution MMS measurements further demonstrate kinetic-scale electric field fluctuations on the low-density side of the edges of surface waves. Detailed comparisons with the simulations suggest that the observed fluctuations are generated by the lower-hybrid drift instability excited by the density gradient at the edges of these surface waves. These secondary effects can lead to a flattening of the edge layers, indicating the connection between kinetic and larger scales within the KH waves and vortices.
The Kelvin–Helmholtz instability (KHI) is a velocity shear driven instability and is considered one of the main candidates for energy and momentum transfer across Earth's low-latitude magnetopause. It has been frequently observed during northward interplanetary magnetic field (IMF) conditions in the flank regions of Earth's magnetopause.1–5 The KHI forms sinusoidal surface waves that travel anti-sunward along the magnetopause. These waves can develop into large-scale rolled-up vortices and efficiently mix plasmas from the two sides of the boundary.6–11 KH vortices are therefore, besides magnetic reconnection, a mechanism for the transport of plasma into the magnetosphere.12 The onset and evolution of this instability is governed by a criterion formulated in ideal MHD for incompressible plasmas and an infinitely thin boundary layer13 subject to a super-Alfvénic shear flow based on the magnetic field component parallel to the wave vector of the KHI given by
where k is the KH wave vector and ρ1, ρ2, , and are the ion densities, ion bulk velocities, and magnetic fields on the two sides of the boundary (1: magnetosheath and 2: magnetosphere), respectively.
Since the magnetic field component in the direction perpendicular to the wave/vortex plane ( north–south direction of Earth's magnetic field) is perpendicular to an idealized wave vector k along the magnetosheath flow direction, the instability criterion will not be affected by a change in the sign of the IMF . However, the KHI has been much less frequently observed during southward IMF conditions.4 The reason why the observational probability of KH waves and vortices is much lower during southward IMF is still under debate. Based on 3D fully kinetic simulations, it was recently suggested that during southward IMF, the KH vortex-induced reconnection (VIR)8,10 quickly destroys the wave and vortex structures and reduces the observational probability of clear signatures of KH waves.14 However, there have been no observational studies to confirm such theoretical predictions. Hwang et al.15 provided the first observations of non-linear KH vortices during southward IMF from Cluster data, which showed the irregular and temporally intermittent behavior of KH waves during southward IMF. Only two studies showed from in situ data that the KH instability can indeed grow during southward IMF conditions at Earth's magnetopause.15,16 These studies were performed using Cluster and THEMIS data.
Considering that recent Magnetospheric Multiscale (MMS) observations and supporting kinetic simulations during northward IMF5,10,17,18 showed the importance of small-scale physics in MHD-scale KH waves and vortices, in this study we present the first high-resolution MMS observations of the KHI during southward IMF. The structures of these KH waves and vortices were found on September 23, 2017, between 15:30 and 16:30 UT.
In our companion paper by Nakamura et al.,19 we report the results of 2D and 3D fully kinetic particle-in-cell (PIC) simulations, performed for the plasma parameters observed during this MMS event. The simulations successfully confirmed the onset and evolution of KH waves and vortices under the observed conditions. They further demonstrated the excitement of secondary instabilities such as the lower-hybrid drift instability (LHDI) and the Rayleigh–Taylor instability (RTI) near the edges of KH structures, where high density gradients can form. While the LHDI can lead to a diffusion of the boundary layer and an enhancement of plasma mixing across the boundary, the RTI can significantly deform the KH structures. Consequently, they can both lead to non-ideal signatures in in situ data, and can help explain the lower observational probability of KHIs during southward IMF. Based on the high time resolution of the MMS mission, this paper shows the simultaneous observations of large-scale KH structures and kinetic-scale fluctuations that can be interpreted as being produced by the LHDI, enabling a multi-scale study of the KHI during southward IMF.
This paper is organized as follows: Sec. II introduces the data sets employed in this work. Section III is split into an overview of the plasma parameters in the solar wind that influenced our KH event and the identification of KH waves in the MMS event. Based on a comparison of the MMS data to the simulation results, Sec. IV shows a more detailed analyses of non-linearly rolled-up vortex structures. Section V then analyzes the observed secondary instability (LHDI). Finally, Sec. VI summarizes the main results of this paper and Sec. VII discusses some additional remarks and an outlook.
II. DATA SETS
In this work, we use data from the plasma and field instruments of the Magnetospheric Multiscale (MMS) mission.20 Magnetic field data are provided from the Fluxgate Magnetometers (FGMs)21 in fast survey and burst mode, at 16 and 128 Hz, respectively. Furthermore, the magnetic field data are used from the Fluxgate–Searchcoil Merged (FSM)22,23 data set, providing a resolution of 8192 Hz. It is a combination of the FGM and search-coil magnetometer (SCM)24 data. The electric field data are obtained from Electric Field Double Probes (EDPs),25,26 in fast and burst mode, at 32 and 8192 samples/s, respectively.
For the particle moments, we employed data from Fast Plasma Investigation (FPI) instruments,27 which calculate the moments of the measured electron (ion) distribution functions in burst mode at a cadence of 30 (150 ms). In the fast mode, the data are sampled every 4.5 s. FPI data can also be generated at a cadence of 7.5 ms (Ref. 28) and were employed in the detailed analysis of plasma waves in Sec. V. If not stated otherwise, only data from the MMS1 spacecraft are shown in the plots.
To analyze the solar wind conditions that influenced our MMS event, data from the Time History of Events and Macroscale Interactions during Substorms (THEMIS) mission,29 as well as data from the OMNI data set,30 are used in this work. The OMNI data set is comprised from observations of several solar wind observing spacecraft, which are time-shifted to the bow shock nose.
The High Resolution OMNI (HRO) data are employed with a cadence of 1 vector/min. For the THEMIS spacecraft, we use data from the Fluxgate Magnetometers (FGMs)31 for the magnetic field at a spin-resolution (FGS) of 1/3 Hz. The ion and electron particle data and moments are obtained from an electrostatic analyzer (ESA)32 at the same resolution.
We transformed the data from the geocentric solar ecliptic (GSE) coordinate system to aberrated GSE (AGSE) coordinates to account for the aberration of solar wind in the rest frame of the Earth. Then, the data were further transformed to a boundary normal (LMN) coordinate system introduced by Russell and Elphic.33 It is defined such that it has its N axis perpendicular to an unperturbed magnetopause boundary, determined from the Shue et al. magnetopause model34 and is directed outward from the Earth. The M axis can be obtained from the cross-product between N and the Geocentric Solar Magnetospheric (GSM) z axis and approximately corresponds to the direction opposing the magnetosheath flow along the dusk flank. The L-component completes the orthogonal right-hand coordinate system, approximately corresponding to the north–south () direction. The new directions have the GSE coordinates of , and .
Finally, the observational data were compared to results from multiple runs of 2D and 3D fully kinetic PIC simulations, performed for the plasma parameters observed during the KH event and reported in our companion paper.19
A. Overview of event
On September 23, 2017, the four spacecraft from NASA's Magnetospheric Multiscale mission20 observed multiple magnetopause crossings between 12:00 and 20:00 UT. During this interval, they were located at the dusk flank of the magnetopause at around [x, y, z] [−6, 10, 21] Earth radii (RE) in GSE coordinates with an inter-spacecraft separation of approximately 20 km. Concurrently, two spacecraft from NASA's THEMIS mission (THEMIS-A and THEMIS-E) were located slightly dawnward of the dayside magnetopause at around [12, −6, −2] in GSE coordinates. The locations of the MMS and THEMIS spacecraft are marked in Fig. 1, where the two thick black lines mark a standard magnetopause and bow shock, estimated from the Roelof and Sibeck35 model and the Fairfield36 model, respectively. Due to their positions near the subsolar point behind the nose of the modeled bow shock, the THEMIS-A/E spacecraft observations are used as a reference for the magnetosheath conditions that influenced the magnetopause dynamics at the time of detection by the MMS spacecraft. Together with the OMNI data set, the evolution of the solar wind and its influence on the magnetopause can be estimated. Note that although THEMIS was located in a different sector (dawnside), the component in the THEMIS data would be useful to check the north–south orientation of the magnetosheath magnetic field in this MMS event in addition to the solar wind OMNI data.
Figure 2 shows observations from the MMS1 and THEMIS-A/E spacecraft, together with magnetic field data from the OMNI data set between 12:00 and 20:00 UT. In panel (a), the energy spectrum of ions, observed by MMS, is shown. Together with panels (b)–(d), it shows multiple transitions between the magnetosphere and the magnetosheath. They are characterized by high-energy, low-density, relatively stagnant plasma with northward (+) magnetic field orientation and lower-energy, high-density, and fast anti-sunward flowing plasma with northward- or southward-oriented magnetic field orientation, respectively. The interval focused on in this work (15:30 to 16:30 UT) is marked by vertical dashed lines and shows dominant southward magnetic field in the magnetosheath intervals. Panels (e)–(g) show the magnetic fields observed by the spacecraft from the OMNI data set in the solar wind and by THEMIS-A/E in the magnetosheath near the subsolar point, respectively. The data are time-shifted to the MMS event by 10 min, assuming an average flow speed along the magnetopause of around 315 km/s and not taking into account the curvature of the magnetopause. This allows for a rough direct comparison of OMNI, THEMIS, and MMS data. Main tendencies, such as an increase in at around 17:30 UT, can be seen at the same time in panels (d)–(g), which supports the reliability of this estimation. When comparing the components of panels (d)–(g), these observations confirm the negative (southward) orientation of the magnetic field of the solar wind and the magnetosheath plasma during the selected interval (15:30 to 16:30 UT).
Although there are several intervals with a southward-oriented magnetic field in panels (e)–(g) before and after the selected interval, they mostly show dominant or components [cf. panels (f) and (g)] or no clear KH signatures, which is why they were not considered in the present study. Furthermore, during most of these intervals, the magnetosheath magnetic field observed by MMS was directed northward.
B. Identification of KH waves
Figure 3 shows the in situ data observed by the MMS1 spacecraft between 15:30 and 16:30 UT. The first two panels show the energy spectra of ions and electrons, respectively. Together with the ion and electron densities in panel (c) and the ion and electron temperatures in panels (g) and (h), they indicate multiple crossings of the magnetopause between the cold and dense magnetosheath and the hot and tenuous magnetospheric plasma. Panel (c) indicates a higher density encounter around 15:57 UT, where the MMS1 spacecraft observed almost pure magnetosheath plasma, as can be seen in panel (a). This encounter of a density of around 12 cm−3 can be estimated as the background density of the magnetosheath during the time of the event. In all the other encounters of the magnetosheath, especially between 16:00 and 16:20 UT, the ion and electron densities are somewhat lowered, indicating a mixing of plasmas across the boundary layer or a crossing line of the spacecraft that was located somewhat closer to the center of the magnetopause layer. Also, since the energy spectrum shows a more frequent encounter of magnetospheric populations at a later time between 15:30 and 16:30 UT, this further indicates that the magnetopause moved outward during the interval. This can be estimated to be a distance of about 104 km between 15:40 and 16:00 UT, given the average observed velocity perpendicular to the magnetopause during this time of km/s [panel (e)].
Panel (d) shows the three components of the magnetic field in LMN coordinates. The L-component varies between positive and negative values in the magnetospheric and magnetosheath intervals, consistent with the southward magnetosheath magnetic field conditions. The M-component of the magnetic field, which can suppress the growth of KH waves, was relatively weak compared to the dominant in both the magnetosheath and magnetospheric intervals, respectively (cf. 15:58 and 16:02 UT). The component of the magnetic field is also significantly lower than the dominant , but it shows positive-to-negative variations before the crossings of the trailing edges (marked by vertical black dashed lines), which could be interpreted as rippled magnetic field lines following the surface wave structure close to the edges on the low-density side of the boundary layer.
In panel (e), clear fast-flowing magnetosheath plasma can be observed in the (−M) direction, while the magnetospheric plasma is relatively stagnant. The velocity component in the L-direction is relatively weak and was therefore neglected in further analysis. However, the velocity component perpendicular to the magnetopause (+N) shows clear positive-to-negative variations before the trailing edge crossings (3–7) and only weakly before crossings (9) and (11) marked above Fig. 3. These variations might originate in the relatively stagnant magnetospheric plasma following the wave structure at the leading and trailing edges, respectively. This feature is only observed in a limited set of crossings, probably due to a missing of the center of the layer as mentioned before. A similar behavior was seen in the N-component of the magnetic field, also showing less clear variation for waves (1) and (2), during which the spacecraft were most likely located on the magnetosheath side of the layer.
To evaluate the instability criterion of the KHI, in Table I the plasma parameters for the ion density, ion bulk velocity, and magnetic field are collected for the two sides of the boundary layer between 15:30 and 16:30 UT. Note that since there is a larger uncertainty for and this in-plane magnetic field component has a great impact on the evolution of the KHI, the error ranges of these components are relatively large. They were obtained from a comparison to simulation results, which revealed an enhancement of the background near the edges of the waves and vortices by a factor of 2 due to the wave/vortex motion. The background BM in our MMS event can then be estimated by dividing the observed enhancements of near the edges between 15:30 and 16:30 UT by a factor of 2 and using the results as an error range.
|Parameter .||Magnetosheath (1) .||Magnetosphere (2) .|
|n (cm−3)||12 ± 1||0.6 ± 0.2|
|v (km/s)||(0, 0, 0)|
|B (nT)||(7 ± 2, 0 ± 4, 0)|
|Parameter .||Magnetosheath (1) .||Magnetosphere (2) .|
|n (cm−3)||12 ± 1||0.6 ± 0.2|
|v (km/s)||(0, 0, 0)|
|B (nT)||(7 ± 2, 0 ± 4, 0)|
The magnetospheric velocity and all vector components in the N-direction were set to zero due to their small values and their weak influence on the result of the instability criterion, which is why only a 2D configuration is considered here. For the instability criterion, the Alfvén Mach number was formulated as follows:
which is obtained by dividing the velocity-dependent term of Eq. (1) by the term dependent on the magnetic field. These two terms represent the velocity shear and magnetic tension present in the system. As shown in previous studies,37,38 this relation can be used as a reference for the onset probability of the KHI, where corresponds to a super-Alfvénic shear flow. To estimate , the angle of the wave vector relative to the M-axis has to be determined. This angle was obtained from a maximum growth rate analysis as described in the following.
The instability criterion [cf. Eq. (1)] can be probed for the angle that maximizes the growth rate for certain input parameters by using the wave vector in LMN coordinates and varying the angle θ. By inputting all combinations of the background plasma parameters of Table I into Eq. (1), we can estimate the angle range that maximizes the growth rate of the KHI during our event for all possible combinations. When including angles of the wave vector from to measured from the positive M axis, i.e., only negative M-components (magnetosheath flow direction), we find an angle range of maximum positive growth rate in the range . When inputting the plasma parameters from Table I together with the estimated angle range that yields the maximum growth rate, the results for the Alfvén Mach number easily exceed the discussed relation . This suggests that the onset condition of the KHI under the observed conditions was easily fulfilled during this southward IMF event observed by MMS.
Together with the variations of , the total (magnetic and plasma) pressure, shown in panel (f) of Fig. 3, can also be used as a parameter to identify KH waves and vortices.39,40 As KH waves grow, the plasma flows are converged toward the trailing edge, increasing the total pressure at this location, while near the center of the waves, the pressure becomes lower. Interestingly, the total pressure shown in panel (f) not only confirms these basic signatures, but also the differences between the discussed intervals with different spacecraft locations relative to the center of the structures. At waves (5) and (7), the very clear KH signatures of a maximum total pressure at the trailing edge and the minimum total pressure near the leading edge or center are confirmed. Between crossings (1) and (2), the variations are significantly smaller, which will be studied in more detail in Sec. IV. This can be related to the offset of the crossing line relative to the center, which will also be discussed from simulation results in Sec. IV A.
The boundary normal vectors, characterizing the local orientation of the magnetopause boundary layer at each crossing, were determined from three different methods to increase their reliability. We applied the timing method41,42 on the electron density data and the timing and Minimum Directional Derivative (MDD)43 methods on the magnetic field data, all in burst mode, respectively.
The multi-spacecraft timing method assumes a one-dimensional discontinuity moving at constant velocity in space with fixed spacecraft positions. Then, by measuring the time differences between the encounters of this discontinuity at each spacecraft and by considering their respective positions, it is possible to determine the orientation (and thus the boundary normal direction) and the velocity of the discontinuity by cross correlation.
The MDD method is based on multi-spacecraft magnetic field measurements to determine the dimensionality and three characteristic directions of a structure in space. In the first step, the magnetic field gradient tensor is determined at every time instance of a defined time interval. Then, the eigenvalues and eigenvectors of the symmetrical matrix are calculated and the maximum, intermediate, and minimum eigenvalues are hereafter known as λmax, λmid, and λmin, respectively. From a comparison of these eigenvalues, the dimensionality of the structure can be found, where corresponds to a one-dimensional structure, to a two-dimensional structure, and to a three-dimensional structure. If the structure can be identified as one-dimensional, the eigenvector corresponding to λmax is defined as the boundary normal direction. Therefore, the MDD method can be used to support the results from the timing method for the boundary normal direction, but not for the estimation of the boundary normal velocity.
The results from the discussed methods are collected in Table S-I given in the supplementary material. The time intervals used for the timing method are also given in Table S-I and are selected individually for each crossing. By varying all time intervals by a few seconds at each crossing, a convergent result for the boundary normal directions can be obtained to reduce errors. The MDD method was applied by utilizing the implemented analysis tool in the GUI of SPEDAS.44,45 For the intervals of each boundary crossing where a one-dimensional structure was identified, we calculated the mean value of the corresponding eigenvector for this interval as implemented in the software.
Panel (m) of Fig. 3 depicts the results for the boundary normal vectors in the M-N plane. Each vector obtained from the timing method (and thus carrying information about the local velocity of the boundary) was normalized individually (cf. Table S-I) to be comparable to the normalized results from the MDD method. Since the lengths of the vectors in the two-dimensional plane might differ from unity, we added a corresponding reference in the panel. At each crossing, the local orientation of the magnetopause boundary layer in the M-N plane can be studied with the help of these boundary normal vectors at the leading and trailing edges, respectively.
Comparing the theoretical prediction of to the boundary normal vectors in the L-M plane shown in Fig. 4, we find that the general trend of the boundary normal vectors qualitatively coincides with the theoretical prediction (marked by red vectors). Note that these comparisons are done in LMN coordinates, which contain an uncertainty, since they are computed on a global scale for a magnetopause model. Nevertheless, the consistencies between the observation data in LMN coordinates and the local simulations of this event shown in Sec. IV strongly suggest the adequacy of this LMN system.
The boundary normal vectors determined from different methods and shown in panel (m) of Fig. 3 are in qualitative agreement with each other at all crossings, except for crossings (3, 8, and 10), where they strongly deviate from each other. However, these crossings were affected by a strongly diffused boundary layer crossing together with large variations (3) and considerably short magnetospheric intervals before the trailing edge crossings (8, 10), respectively, which can lower the reliability of the determined boundary normal vectors. Overall, the boundary normal vectors indicate an alternating orientation of the boundary in the N-direction at the leading and trailing edges (light and dark gray dashed vertical lines), corresponding to sinusoidal waves being crossed through their center by the spacecraft. As an example, waves (7–9) show clear wave-like patterns in the boundary normal vectors and can therefore be regarded as rippled surface waves. In combination with the positive-to-negative variations of and , these signatures indicate the detection of magnetopause surface waves. In contrast to this, the boundary normal vectors of waves (6 and 10) point into the same N-direction at the leading and trailing edges, respectively, which may indicate waves that transition into their non-linear stage, steepening the leading edge.46 However, the interpretation of boundary normal vectors can be misleading, when the line of crossing is not taken into account and can indicate wrong stages of development. This will be discussed in more detail in Sec. IV and Fig. 8.
C. Properties of the observed KH waves
Hwang et al.15 have shown that the Kelvin–Helmholtz instability shows an irregular and temporally intermittent behavior during southward IMF conditions. To test this feature, panel (l) of Fig. 3 shows a wavelet analysis of the ion temperature , which reveals a change in the wave period over time. The main wave periods lie in the range of 250–400 s. A more detailed estimation was done in Fig. 5, where the power spectra of the total pressure and the ion temperature are shown in panels (a) and (b), respectively. The peaks of the power spectra, observed by MMS1–4, indeed lie in the same range with values of 238 and 400 s, respectively. However, multiple peaks with comparable power can be observed, strongly indicating an irregular behavior. The first peak of around 1 mHz ( min) reflects the interval near crossings (1) and (2) in Fig. 3, during which the spacecraft were most likely located far from the wave/vortex center on the magnetosheath side at the time of crossing and only crossed the KH structure partially. Therefore, it is excluded from the analysis. Also, when computing the time period between the leading and trailing edges of each wave/vortex, we find strong variations in the wave period.
Note that the superposition of several wave frequencies observed in Fig. 5 might also originate from the varying crossing line of the spacecraft, which is supported by the magnetopause moving outward during our event. However, previous studies of KH waves during northward IMF5 made by the MMS spacecraft showed a more steady behavior of the KHI during northward IMF. Therefore, it is strongly suggested that the KH instability develops an irregular structure during southward IMF, as also suggested by Hwang et al.,15 contributing to the lower observational probability during these conditions. Since there were no significant variations in the density, dynamic pressure, flow speed, and magnetic field intensity of the solar wind before and during the event, the variations observed in Fig. 5 were most likely not induced by solar wind driven external variations of the magnetopause.
To estimate the phase speed of the observed KH waves, the boundary normal vectors determined from the timing method on the magnetic field and containing the boundary normal velocity, shown in panel (m) of Fig. 3 and collected in Table S-I, are deprojected onto the L-M plane, since the wave vector of the KHI is assumed to lie in this plane. The concept behind deprojection47 can be seen in Fig. 6, which shows that the local speed determined from the timing method is directed normal to the local boundary layer. Assuming that the rippled surfaces propagate in the (−M) direction, as the surface waves propagate along the boundary, we take the deprojected M-component of the timing velocity as the wave propagation speed of the surface waves. Deprojecting all boundary normal vectors determined from Timing (B) onto the L-M plane and taking the mean value of the M-components, the phase speed can be estimated to km/s. Together with the estimated wave period of s (cf. Fig. 5), the wavelength of the KH waves can approximately be estimated to km, which is comparable to the values seen in previous studies for northward IMF cases.1
Note that the range of observed wave periods in the present study is approximately twice as large as in a previous study on the KHI during southward IMF by Hwang et al.15 In addition, the magnetosheath flow velocity v0 reported in the study by Hwang et al. exceeds the one observed in the present study (cf. Fig. 3) by a factor of 2. The sampling positions are given as and in GSE coordinates. By estimating the growth time of the KHI as with the distance between the onset and observation of the KHI and further normalizing by the growth time of the KHI used in previous comparison studies10,48 , we find . Since is comparable in both studies (assuming a similar onset location of the KHI in both cases), we estimate and find comparable wavelengths, suggesting a similar growth phase of the KHI in both studies (late linear to early non-linear stage).
To check if the assumption of neglecting the L-components is well justified, we plotted the boundary normal vectors in the L-M plane in Fig. 4 together with the theoretical prediction of . Overall, the vectors seem to be approximately within this range.
IV. COMPARISON TO SIMULATION
We compare the observed MMS event of the KHI with data from a 2D fully kinetic PIC simulation performed by Nakamura et al.19 In the simulation setup, the (X, Y, Z) coordinates approximately correspond to the magnetosheath flow (−M), anti-boundary normal (−N), and the south-to-north direction (L) of the LMN coordinate system. Denoting the magnetosheath and the magnetosphere as 1 and 2, respectively, the parameters used to set up the simulation were with cm−3, nT and nT, and , with , with the local Alfvén velocity. Note that the two half-spaces are set to the frame with half the velocity of the velocity shear. The parameters were obtained from the interval between 15:33:00 and 15:34:15 UT, corresponding to the KH structure (1) in Fig. 3.
Although some different sets of the initial in-plane field ( and ) are tested in our companion paper,19 in this paper we especially focus on a large-scale 2D run in which nT and . The system size of this run is cells with a total of simulated particles. More details on the simulation setup can be found in our companion paper.19
The simulation demonstrates a strong evolution of the KHI, which supports that the magnetopause would indeed be unstable to the formation of KH waves under the conditions of the present MMS event. Figure 7 shows an overview of the results from this simulation in the early non-linear growth stage of the KHI (). We see that clear surface waves generated by the KHI have already developed at this time and some of these waves begin to form rolled-up vortex structures. Three virtual spacecraft orbits are depicted in the plot by the horizontal dotted lines. The direction of the spacecraft crossings is set from right-to-left, since the KH structures travel into the (−M) direction in the realistic setup.
Additionally, a zoom-in view of the area marked by the blue box is shown in panel (b). Two vortex structures and the virtual spacecraft crossing orbit-2 are shown together with schematic boundary normal vectors of the first vortex. This vortex structure will be studied in more detail in the following.
A. KH vortex detection
Figure 8 shows MMS observations between 15:33:00 and 15:34:15 UT in burst resolution and a comparison to virtual observations (orbit-2) from the simulation. In this plot, the normalizations for the simulation data are cm−3 for the ion density, nT for the magnetic field, km/s based on and for the ion velocity, and nPa for the total pressure, corresponding to . The time steps are given as multiples of the ion gyro-period () and the (X, Y, Z) coordinates in multiples of the ion inertial length ().
Panels (a) and (b) show the ion energy spectrum and density, with the crossing of the leading edge of the vortex at around 15:33:28 UT and the trailing edge at around 15:34:00 UT, marked by vertical dashed lines. Interestingly, at around 15:33:45 to 15:33:50 UT, a short encounter of magnetosheath-like higher-density and colder plasma was observed on the low-density side of the vortex. In panel (c), the L-component of the magnetic field, indicating the north–south component, also decreases during this short encounter of the magnetosheath-like plasma. All of these signatures most likely indicate the crossing of a vortex arm, where the high-density plasma is dragged into the center of the vortex.
In panel (a), a high-energy, low-density plasma population was observed between approximately 15:33:28 and 15:33:38, which shows the same ion density as magnetospheric plasma and a strong northward magnetic field orientation. This comes together with an ion bulk velocity exceeding the magnetosheath velocity of km/s to peak values of around −400 km/s. As shown in previous studies,1,50 this Low-Density Faster Than Sheath (LDFTS) plasma can be an indicator for the observation of a vortex structure (see Sec. IV B for more details on the LDFTS signatures).
Comparing these signatures to the virtual crossing of a vortex from the simulation on orbit-2 of Fig. 7, marked by vertical dashed lines, we find many consistencies. In Fig. 8, panels (f) and (g) show a simultaneous increase in the ion density and decrease in () at the crossing of a vortex arm, encountered at around in Fig. 7(b). Furthermore, panel (h) confirms the detection of LDFTS plasma, however located closer to the vortex arm and trailing edge of the structure than in the MMS data.
It is notable that the total pressure observed in MMS data [Fig. 8(e)] is relatively constant. This is in contrast to an expected maximum of the total pressure at the trailing edge, as seen at crossings (5) or (7) in Fig. 3(f). Interestingly, this almost constant total pressure can also be observed in the virtual crossing of the simulated vortex [Fig. 8(i)]. When compared to Fig. 7(b), the marked vortex is crossed on its magnetosheath side. The center of the structure is somewhat missed and results in a loss of clear signatures.
Note that the variation in ion density in the simulation [Fig. 8(f)] is more gradual in the leading edge than in the trailing edge in contrast to the MMS observations [Fig. 8(b)]. A thinner trailing edge can more easily form because the vortex motion tends to compress the trailing edge more strongly due to converging plasma flows as shown in previous studies46 and observed in the present simulation. Therefore, the observed thicker trailing edge seen in panel (b) is a unique feature of this observation interval. Our companion paper suggested that the lower hybrid drift instability (LHDI), generated by the density gradient across the edge layers, can more or less diffuse the layers.19 In Sec. V, we show the LHDI activity and the relation to a diffuse edge layer by plasma mixing across the boundary. However, the smaller ion-to-electron mass ratio used in the simulation and the large difference between the simulation box size and the real size of observed KH structures leads to differences in the relative scales of the primary (KHI) and secondary (LHDI) effects between the simulation and the observations. Due to this larger difference in scales between the LHDI and the primary KHI in observations, in this observation interval, the LHDI can have more time to evolve during the growth of the KHI and therefore diffuse the boundary layer more strongly. It may, however, not reach the same stage in the simulation to strongly diffuse the trailing edge layer, which might explain the difference between the simulation and the observations.
B. LDFTS plasma
A detailed look at the detection of LDFTS plasma shown in Fig. 8 is provided in Fig. 9. Following the ideas of Takagi et al.,50 we plotted the ion velocity in the M-direction vs the ion density, for both the observational and simulation data. The observational data for this study are taken from the whole interval between 15:30 and 16:30 UT, shown in Fig. 3. The simulation data are taken from the whole simulation domain of Fig. 7(a), including the regions above , in which the three marked virtual crossings (orbits 1–3) are marked as green, red, and blue dots in Fig. 9(b), respectively.
To perform this LDFTS analysis, we set the boundary for low-density plasma to the value of cm−3 ( of the maximum ion density during the interval from 15:30 to 16:30 UT). The magnetosheath velocity is defined as the mean velocity of all data points of , which are observed during intervals with . Then, the standard deviation σ of the mean value is added to and the larger absolute value ( km/s is taken to be the cutoff for the magnetosheath velocity. In summary, the data points exceeding the estimated magnetosheath velocity and lying below the threshold of are defined to be of LDFTS origin.
Figure 9(a) shows the results for the MMS observations between 15:30 and 16:30 UT, with cm−3 and km/s. The red squares indicate plasma of low density exceeding the magnetosheath velocity and considered to be of LDFTS origin. When checking the times of detection of these data points, they coincide with the data points observed at the leading edge of the vortex structure, discussed in Fig. 8(d). Therefore, the low-density fast-flowing plasma detected during the event shown in Fig. 8 was most likely LDFTS plasma originating from vortex motion. The fact that we only observed a small fraction of LDFTS plasma in our event [cf. Fig. 9(a)] is most likely due to the much lower threshold for low-density plasmas compared to previous studies3 and the possible strong deformation of vortices due to the RTI as suggested by our companion paper.19 Since the RTI can lead to a strong penetration of vortex arms perpendicular to the magnetopause boundary, it can disturb the rolling up of vortices and therefore reduce the LDFTS plasmas.
Panel (b) of Fig. 9 shows the same analysis for the simulation data, with the boundary values of MMS data normalized to the simulation setup. For this simulation setup, the density ratio between the low- and high-density sides was set to be . Therefore, normalizing the low-density boundary cm−3 by , we obtain for the analysis of simulation data. The magnetosheath velocity in the simulation was set to be .
When comparing the data points obtained from different crossings, it becomes clear that crossings on the magnetosheath side of the vortex structures are much more likely to show LDFTS signatures [cf. orbits 1 and 2, green and red dots in subplot (b)]. A crossing through the center of the vortices or on the magnetospheric side [cf. orbit-3, blue dots in subplot (b)] does not produce LDFTS signatures, considering our boundary values.
The preferred observation of LDFTS plasma for crossings on the magnetosheath side of vortices further strengthens our arguments of a magnetosheath side crossing line in the observed vortex structure in MMS data in Fig. 8.
V. LOWER HYBRID WAVES
Recent studies have shown the interplay between the LHDI and the KHI from 2D fully kinetic simulations using parameters typical for the magnetopause of Mercury.51 Depending on the density gradient and velocity shear across the boundary, the LHDI and the KHI can compete for their respective growth. If the density gradient is strong enough, the LHDI disturbs the structure of the shear layer and can even suppress the onset and evolution of the KHI. For this situation to happen, the boundary layer width should not exceed a few ion inertial lengths, which is usually not fulfilled at Earth's magnetopause. However, when large-scale KH vortices form at Earth's magnetopause, they locally compress the boundary layer at the trailing edges due to vortex motion and can thus, in combination with a large enough density gradient, lead to a strong growth of local lower hybrid waves.
Utilizing the unique high-resolution capability of the MMS mission, we were able to confirm enhanced wave activity on the low-density sides of the trailing edges of the KH waves and vortices shown in Fig. 3. A first indication of plasma waves is given in panels (i)–(k), where strong fluctuations in the electric field near the trailing edges (indicated by darker vertical dashed lines) can be seen together with enhanced wave power in the electric field in panel (j). The white lines in panels (j) and (k) indicate the local lower hybrid frequency estimated from the relation
where , and correspond to the ion cyclotron, electron cyclotron, and ion plasma frequencies, respectively. The maximum wave power can be observed around this frequency. Enhanced wave activity can also be observed near the leading edges of some waves. However, due to a weaker density gradient and less compressed boundary layer at the leading edges, the lower hybrid wave activity is weaker in this region. Interestingly, panel (k) shows almost no wave activity in the magnetic field in this frequency range, indicating that the observed waves are almost electrostatic.
Since plasma waves are often characterized by their relative orientation to the magnetic field, we transformed the electric field data to a field-aligned coordinate (FAC) system. It is defined by the unit vectors , with being the ion bulk velocity and B the magnetic field. The transformed electric field data correspond to the perpendicular and parallel components . The ion velocity was chosen for this transformation to minimize the noise present in the electron velocity data. To keep the resolution of the burst electric field data (8192 samples/s), the burst velocity data, which have a lower resolution, are linearly interpolated to the electric field data resolution. Each electric field component, namely, , and Epara, was transformed individually. For analyzing the wave activity, we selected the KH structure (4) of Fig. 3 due to its large density gradient and clear transition between the two sides of the boundary, which is shown in Fig. 10.
A. MMS event analysis
Figure 10 shows a typical crossing of the trailing edge of a KH wave during 15:55:40 to 15:56:25 UT. In panels (a) and (b), which show the energy spectrum of ions and electrons, respectively, we see a transition from the magnetosphere-like to the magnetosheath-like conditions at the crossing of the boundary layer on the low-density side ( 15:56:08 to 15:56:10 UT). This corresponds to the enhancement of the ion and electron densities in panel (c) that exceed the background magnetospheric density. In panel (d), a southward turning of the magnetic field can be observed across the boundary layer in conjunction with a negative-to-positive variation of , which was identified as rippled magnetic field lines following the wave structure close to the boundary layer in Fig. 3 of this paper. The ion and electron velocities in panels (e) and (f) show a basically quiet behavior with relatively stronger fluctuations in the electron velocity components, especially on the low-density side.
Figures 10(g) and 10(h) show the perpendicular and parallel components of the electric field, respectively. We see increased activity in the electric field components perpendicular to the magnetic field during the interval marked by two thick vertical dashed lines. The amplitudes of these fluctuations exceed the fluctuations in the parallel electric field by more than two orders. This is in good agreement with a characteristic feature of the LHDI.52 Interestingly, the electric field fluctuations are not only seen in the region with large gradients in the ion density and the magnetic field, but extend far into the low-density side of the boundary. These features are in reasonable agreement with results from our 3D simulation shown in Fig. 6 of our companion paper19 and introduce turbulence far into the low-density side of the boundary layer.
To clearly mark the transition from the low- to the high-density side across the current sheet, we plotted the current density components (GSE) calculated using the curlometer technique53,54 in panel (l), where the peaks of the individual components at around 15:56:11 mark the crossing of the current layer of the magnetopause. The interval marked by the light yellow box near the boundary layer will be used for further analysis of the waves closest to the boundary layer in Fig. 11.
Panel (m) shows the wave power spectrum of the electric field, indicating increased wave power near the lower hybrid frequency (marked by the black line). In panel (n), the power spectrum of the magnetic field shows almost no wave activity near the lower hybrid frequency, characterizing the observed plasma waves as electrostatic, which is in good agreement with previous studies of this characteristic property of the LHDI.55 This has important consequences on our previous analysis of the boundary normal vectors based on the timing and MDD methods, both utilizing the magnetic field data—i.e., the electrostatic character of the small-scale waves would not significantly disturb the magnetic field, enabling the timing and MDD methods to identify the larger-scale primary KH waves. Although we observe an enhancement in the wave power of the magnetic field in panel (n) around 15:56:10 UT, the amplitudes of fluctuations in the electric field relative to the background electric field strength still exceed the relative amplitudes of the magnetic field fluctuations by more than two orders of magnitude. Furthermore, no clear signatures for magnetic reconnection were found during this interval.
Panel (i) shows a gradient in the electron pressure near the boundary layer, coinciding with the onset of wave activity as well as the density gradient. As seen in panels (e) and (f), since the current (dominated in the x-components) is carried mainly by electrons near the gradient layer, this electron pressure gradient could serve as an energy source for the observed lower hybrid waves.56
Lower hybrid waves are often associated with particle heating.55 Panels (j) and (k) show the parallel and perpendicular components of ion and electron temperatures, respectively. Although a slight enhancement of the parallel electron temperature can be seen during the interval of increased wave activity and an increased perpendicular component of the ion temperature far from the boundary can be observed, no clear evidence for strong particle heating can be found.
Finally, note the local peak of the ion temperature in the yellow shaded region together with a plateau of the ion and electron temperatures between 15:56:08 and 15:56:10 UT. These signatures are seen in most trailing edge crossings that feature lower hybrid wave activity. Although lower hybrid waves can lead to a heating of the plasma, these signatures are most likely associated with mixing of plasmas due to the LHDI. As can be seen in panels (a)–(c) of Fig. 10, mixed plasma populations are observed in the energy spectra of ions and electrons together with slightly increased densities in the respective interval. At the time of local peak in the ion temperature, there is a weak negative enhancement of the densities in panel (c). A comparison to the 3D simulation results showed a similar local peak of the ion and/or electron temperatures together with a slightly decreased density due to the mixing of plasmas. This is caused by different mixing fractions of the hot and tenuous magnetospheric and the cold and dense magnetosheath plasma, which can locally decrease (increase) and thus lead to the increased (decreased) local temperatures.
Although we only show one crossing in detail, panels (j) and (k) of Fig. 3 clearly show that similar lower hybrid wave signatures can be found at almost all other crossings, showing the ubiquity of the LHDI at the edges of the KHI. We analyzed many of the crossings in the same detail as we did in Fig. 10 and found similar signatures throughout the KH event, however mostly concentrated near the low-density side of the trailing edges. This comes as no surprise, since the trailing edges are mostly much more compressed by KH waves and we find stronger density gradients, than at the leading edges of KH structures.
B. Linear dispersion comparison
In our companion paper, a linear dispersion relation solver in the fully kinetic regime (see Umeda and Nakamura57 for more details about this linear solver) was applied to electric field fluctuations within KH waves.19 It utilized both simulation and observation parameters from the present studies to test the LHDI criterion and determine the onset conditions of these small-scale fluctuations. The results showed that the wave normal angles of the fastest growing modes in both simulation and observation feature a strong perpendicular orientation with respect to the local magnetic field. In one run, the input parameters were set to best fit most of the LHDI events in our KH interval. The results indicated a wavelength of the lower hybrid waves in the present MMS observations of km at onset and an increase by a factor 1.5–2 as the LHDI develops. Overall, the estimated wavelength lies in the range of 50–100 km, depending on the state of evolution of the lower hybrid waves. The estimated phase velocity and frequency of the observed waves were found to be km/s and Hz.
When comparing results from the linear dispersion analysis to our MMS event, we need to take into account a possible Doppler shift and write , where the observed frequency depends on the real frequency of the waves (in the rest frame) and the background flow velocity of electrons . If the flow velocity term largely exceeds the rest frame term, then we can apply the Taylor hypothesis58 and neglect the rest frame term.
To this end, we transformed the electron velocity data near the boundary layers of several crossings that feature clear LHDI signatures and consistent boundary normal vectors from all three methods, to the direction perpendicular to the local boundary normal vector n and the magnetic field B, i.e., . This most likely corresponds to the direction of the wave vector of the lower hybrid waves, since they will move along the boundary layer nearly perpendicular to B. The flow velocity along the boundary layer perpendicular to B was found to always lie in the range of 200–300 km/s [cf. Fig. 10(f)], clearly exceeding the estimated phase velocity from the linear dispersion analysis. Utilizing the Taylor hypothesis, we can write and estimate the wavelength of the observed lower hybrid waves in MMS data by inputting the observed frequency.
To obtain a reliable frequency of the observed waves near the boundary layer, we plotted the electric field fluctuations observed by MMS1–MMS4 during the interval marked by the light-yellow box of Fig. 10 in panels (a)–(d) of Fig. 11. We can see that clear fluctuations are visible during the interval marked by vertical dashed lines in Fig. 11. For this interval, we performed a power spectrum analysis on the data of all four MMS spacecraft and plotted the results in Fig. 12. The main peaks of all four spacecraft are located within a range of Hz. Plugging these numbers into our dispersion relation, we obtain a wavelength of the observed lower hybrid waves in the range of km. This is in excellent agreement with the wavelength estimated from the linear dispersion analysis.
Note that the estimated wavelength of the observed lower hybrid waves ( 60–100 km) is a few times larger than the inter-spacecraft separation ( 20 km). From this rather small difference, we expect to see different phases (waveforms) of lower hybrid waves among the four MMS spacecraft. This is confirmed to some degree in Fig. 11, where the four spacecraft indeed observed different waveforms. These different waveforms make it difficult to apply the timing method to the electric field data to determine the wave vectors of the observed lower hybrid waves directly.
Finally, panel (e) of Fig. 11 shows a detailed plot of the electric field fluctuations in the x-direction (GSE) observed by MMS1, together with the ion and electron convection terms and . The data for this analysis were obtained from the new highest-resolution FPI data28 with a cadence of 7.5 ms. The electric field data were downsampled to the same cadence as the electron moments by linear interpolation and the convection terms were scaled to the peak of the electric field during this interval. A clear trend of the electron convection term following the electric field fluctuations can be seen and thus , which is in good agreement with previous studies of magnetopause lower hybrid waves by Graham et al.55
VI. SUMMARY AND DISCUSSION
In this work, we showed the first MMS observations of the Kelvin–Helmholtz instability during southward IMF. The event was observed on September 23, 2017, at the duskside magnetopause and solar wind data from the OMNI data set and the THEMIS-A/E spacecraft confirmed the southward orientation and the dominant component in the solar wind and the magnetosheath during the event. The instability criterion for the KHI was fulfilled throughout the event for the local plasma parameters collected in Table I.
A boundary normal vector analysis showed a wavy structure of the magnetopause throughout the event and indicated, together with the detection of a KH vortex, that the stage of development of the KHI was most likely in the late linear to the early non-linear stage. Note that the KH waves most likely had their onset in a more upstream region and the local plasma parameters observed by MMS would be more or less globally changed at the onset location.2,59 However, the agreement with the boundary normal vector analysis and with the simulation results discussed in Sec. IV suggests that the estimation of the onset conditions and the related maximum growth rate is reasonable given the local parameters observed by MMS. The orientation of the boundary normal vectors in the wave vector plane was also in qualitative agreement with results from a maximum growth rate analysis.
Clear variations of indicative parameters for the KHI, e.g., the total pressure or the normal component of the ion velocity, , were only observed in a limited time range. The reason was found to be the line of crossing of the spacecraft relative to the center of the KH structures and was confirmed from a comparison to 2D simulation results. Furthermore, a power spectrum analysis confirmed the irregular behavior of the KHI during southward IMF, leading to a lower observational rate during these conditions, as shown in a previous study by Hwang et al.15
Comparisons with a 2D fully kinetic PIC simulation, performed for the plasma parameters observed during our MMS event, confirmed the detection of a Kelvin–Helmholtz vortex in MMS data between 15:33:00 and 15:34:15 UT together with the crossing of a vortex arm. Signatures consisting of the relatively constant total pressure, the detection of LDFTS plasma, and matching boundary normal vectors between MMS data and simulation all strongly indicated a crossing of the spacecraft on the magnetosheath side of the vortex, where LDFTS plasma is preferentially observed, as discussed in Sec. IV.
As shown in Sec. V, lower hybrid waves were found to develop at the trailing edges of the KH waves and vortices due to a strong density gradient at the boundary layer, which can extend deep into the low-density side of the structures. They are characterized by electric field fluctuations perpendicular to the magnetic field, strong gradients in the electron density and the electron pressure and . We found that the wave power of the perpendicular electric field fluctuations exceeded that of the parallel electric field by an order of magnitude. Further, the waves were found to be mostly electrostatic, making our boundary normal vector estimations (using the magnetic field) more reliable. The frequency of the lower hybrid waves was found to lie in the range of Hz and their wavelength in the range of km, in good agreement with the estimations from a linear dispersion analysis performed in our companion paper.19 Signatures for particle mixing across the boundary were found together with enhanced lower hybrid wave activity.
The secondary instabilities observed in our event (LHDI) and our companion paper (LHDI, RTI)19 were found to significantly disturb the evolution of KH structures from the early to later non-linear stage. They can lead to a disturbed evolution of the KH structures, to diffuse boundary layers and to non-ideal signatures in in situ data. Their early onset and strong disturbance may help explain the low observational rate of KHIs during southward IMF.
The presented data on the first MMS event of the KHI during southward IMF enables the advancement of research on this phenomenon focusing on both large-scale and high-resolution small-scale kinetic effects. Therefore, it constitutes an important extension of the existing in situ data set of the KHI. The detection of lower hybrid waves at the edges of KH structures serves as a first step for a more extensive research of this event, which will be introduced in Sec. VII.
VII. ADDITIONAL REMARKS
Considering past kinetic studies of the KHI during northward IMF and the subsequent small-scale physics such as turbulence,17,60,61 reconnection,10,62 and kinetic waves,18,63,64 this MMS event can serve as an ideal case study for such high-resolution research studies of the KHI during southward IMF.
The 3D simulation in our companion paper showed signatures of magnetic reconnection within the LHDI turbulence excited at the edge layers of primary KH waves and vortices.19 Although the slight enhancement of the parallel electron temperature seen within the lower hybrid wave interval in the present MMS observations (Fig. 10) may result from reconnection, we have not found any clear signatures of reconnection during the interval. The occurrence of reconnection and related or not-related particle heating within the KH-induced lower hybrid waves would be an important future research target.
The simulations also showed that the intermittent and irregular structures of the KH layer and the resulting lower observational probability of the periodic primary KH waves result from the Rayleigh–Taylor instability (RTI). It deforms the KH structures and leads to a deeper penetration of high-density plasma into the magnetosphere. This way, the vortices cannot roll up and the vortex arms grow in the direction perpendicular to the magnetopause. Especially during southward IMF, the onset of the RTI can be strongly enhanced, since the additional current flow in the x-direction (due to the strong magnetic field gradient in the z-direction) can act as an additional centrifugal force driving the evolution of the RTI arm heads. This deformation can significantly lower the observational rate of the KHI during southward IMF. Although the observed irregular and intermittent behavior of surface waves indicates the evolution of RTI arms during this MMS event, to directly identify the large-scale RTI effects in our data, an analysis method that can account for the small spacecraft separation of around 20 km during our event has to be established.
As a final remark, during the same magnetopause crossing event on September 23, 2017, we found multiple encounters of magnetosheath-like and magnetosphere-like plasmas even during northward IMF intervals before and after the southward IMF interval studied in this paper. They showed less variation in the time period between the crossings and less wave activity in the electric field near the boundary layer, maybe indicating further differences between northward and southward IMF conditions not only concerning the primary KH waves but also the secondary instabilities. Although we have not identified clear KH waves during these intervals, comparing the different IMF conditions during the same series of magnetopause crossings would also be an important future research target of this MMS event.
See supplementary material for detailed results of the methods determining the boundary normal vectors in Table S-I.
This work was supported by the Austrian Science Fund (FWF), Grant No. P32175-N27. For the simulations employed in this paper, we acknowledge PRACE for giving us access to MareNostrum at the Barcelona Supercomputing Center (BSC), Spain. A part of the simulation data was analyzed with resources at the Space Research Institute of the Austrian Academy of Sciences. The work by H.H. was supported by JSPS Grant-in-aid for Scientific Research KAKENHI 21K03504. J.E.S. was supported by the Royal Society University Research Fellowship URF\R1\201286. The observational data employed in this paper were obtained from the MMS spacecraft and are publicly available via NASA resources and the Science Data Center at CU/LASP (https://lasp.colorado.edu/mms/sdc/public/). They were analyzed using the SPEDAS44 software package for IDL (http://spedas.org/blog/). The wave-power spectra in Fig. 10 were obtained using the wavpol.pro routine of SPEDAS.
We thank the entire MMS team for their work on the mission. We are also thankful to the Science Data Center for providing us with the highest-resolution FPI data for this KH event. Further, we thank Quanqi Shi and his Ph.D. student Shichen Bai for their help on the usage of the MDD method, Stefan Eriksson for his help on the power spectrum analysis of the KHI, and Paul Edwin Curtis and Kim Hyuntae for valuable input on the data analysis and interpretation of the results.
Conflict of Interest
The authors have no conflicts to disclose.
The data that support the findings of this study are available from the corresponding author upon reasonable request.