The Kelvin-Helmholtz Instability From the Perspective of Hybrid Simulations

The flow shear-driven Kelvin-Helmholtz (KH) instability is ubiquitous in planetary magnetospheres. At Earth these surface waves are important along the dawn and dusk flanks of the magnetopause boundary while at Jupiter and Saturn the entire dayside magnetopause boundary can exhibit KH activity due to corotational flows in the magnetosphere. Kelvin-Helmholtz waves can be a major ingredient in the so-called viscous-like interaction with the solar wind. In this paper, we review the KH instability from the perspective of hybrid (kinetic ions, fluid electrons) simulations. Many of the simulations are based on parameters typically found at Saturn’s magnetopause boundary, but the results can be generally applied to any KH-unstable situation. The focus of the discussion is on the ion kinetic scale and implications for mass, momentum, and energy transport at the magnetopause boundary.


INTRODUCTION
From the perspective of hydrodynamics, any flow shear with a deformed interface develops pressure gradients following from the Bernoulli principle. A positive feedback yields further deformation leading to the Kelvin-Helmholtz (KH) instability where the boundary interface can be rolled up into a vortex structure. The planetary magnetopause boundary is a prime example of a potentially KHunstable interface, provided that there is sufficient kinetic energy in the flow to overcome magnetic tension forces. Thus, the KH instability has generated considerable attention as a mechanism for facilitating a viscous-like interaction of the solar wind with planetary magnetospheres (review by Johnson et al. (2014) and references therein).
While magnetopause boundary processes garner much of the attention, it should be noted that a discontinuous flow shear with sufficient kinetic energy density to overcome magnetic tension will be KH unstable (Chandrasekhar, 1961). For the simple case of an incompressible plasma with a discontinuous flow shear, the condition for KH instability is where v denotes velocity, ρ is the mass density, B is the magnetic field, k is the wave vector, and the indices refer to the two regions across the shear flow boundary. The corresponding linear growth rate of the wave is given by rates will be a function of the magnetosonic Mach number with most unstable modes occurring when the wavelength is comparable to the width of the shear layer, or 2kL o ∼ 1.
The flow shear on the dawn and dusk flanks of the terrestrial magnetopause boundary are inherently KH unstable, leading to global-scale KH vortices advecting down the magnetotail (Figure 1). In the broader context of planetary magnetospheres, KH waves are ubiquitous and are thought to play an even more important role at the giant planet magnetospheres where the viscous-like processes are thought to play a critical role in the solar wind interaction (Delamere and Bagenal, 2010). Figure 2, from a multifluid Lyon-Fedder-Mobarry simulation of Jupiter's magnetosphere, shows KH waves on the dayside magnetpause boundary due to internal corotational flows (Zhang et al., 2018). It should be noted that KH waves are not found exclusively at the magnetopause FIGURE 1 | GAMERA simulation of Earth's magnetosphere (Sorathia et al., 2020). The color depicts the residual magnetic field (with the dipole subtracted). boundary. For example, secondary KH instabilities could occur at the boundaries of injection channels associated with, e.g., tail reconnection events.
The KH instability is rife with kinetic-scale physics. In fact, KH vortices could be considered as a microcosm for kinetic-scale plasma physics. This papers explores the wealth of KH-related plasma processes with specific application to planetary magnetospheres. While the goal is to explore implications for the terrestrial magnetosphere, numerous examples from the giant planet magnetospheres will be discussed to further understand the kinetic aspects of the KH instability. Specifically, we discuss KH growth characteristics, magnetic field topology, heavy ion effects, momentum transport, particle transport, energy tranport and heating, and electron energization due to parallel electric fields associated with driven KH reconnection.

HYBRID SIMULATIONS
The basis of our discussion will be hybrid (kinetic ions, fluid electrons) simulations of the KH instability Delamere et al., 2018;Delamere et al., 2021). Hybrid simulations push ions according to the Lorentz force equation using the Boris (1970) method. The electric and magnetic fields are ordered on a Yee grid to facilitate the correct field topology for curl operations (Yee, 1966) and to guarantee a divergence-free magnetic field. The magnetic field is updated with Faraday's law using a secondorder, predictor-corrector method. The electric field is determined from the electron momentum equation (massless electrons) that contains an ion-electron collision term and an electron pressure term. The electron fluid is assumed to be isothermal in time and space for numerical simplicity. Ampere's law is used to calculate the electron bulk flow used in the electric field calculation.
With the addition of the electron pressure term, the hybrid simulation is capable of modeling kinetic Alfvén waves (KAW). The dispersion relation for the KAW is where ρ i is the ion Larmor radius, v A is the Alfvén speed. At ion inertial scales, KAW modes capture Alfvén and ion-ion hybrid resonances of the fast magnetosonic/whistler branch. In addition, mode conversion of the compressional fast mode waves to KAW will be captured in hybrid simulations (Lin et al., 2010;Lin et al., 2012). To elucidate the ion kinetic effects, the simulations have been typically conducted for plasma β ∼ 1 where both ion inertial and gyroradius effects are important. By expressing the plasma β as the ratio of thermal pressure to magnetic field pressure, it can be shown that β/2 ρ 2 i /λ 2 i where ρ i is the thermal ion gyroradius and λ i is the ion inertial length. For electron kinetics, we limit the studies to test particle simulations of electron motion along the magnetic field subject to parallel electric fields.
All simulations are based on typical conditions for Saturn's magnetopause boundary. The magnetospheric (corotational) flows are in the +x direction, the +z is normal to the magnetopause boundary from magnetosheath to magnetosphere, and +y is the direction of the magnetic field. Protons are used by default and admixtures of heavy ions (where relevant) were added to the magnetospheric side of the boundary. The typical density is 0.4 cm −3 with an ion inertial length of 360 km. The magnetic field was of the order of 5 nT with a small (e.g., 2%) in-plane component included in the x direction. The velocity shear was of the form where v A is the Alfvén speed, and L o ∼ λ i . The grid resolution was ∼ 0.5λ i . Boundary conditions perpendicular to the flow shear boundary were open while the remaining boundary conditions were periodic. We note that periodic boundary conditions are a limitation for understanding global scale evolution of the instability where field lines can be tied to the ionospheric boundary. Nevertheless, these local simulations should be considered valid during the initial non-linear evolution of the KH instability. Typical spatial domains spanned a few planetary radii.

KH GROWTH CHARACTERISTICS
The fastest growing wavelength of KH modes is typically at the ion kinetic scale. Boundaries in space plasma can be comparable to the ion inertial length (e.g., magnetopause boundary); therefore, the wavelength of fastest growing wavelength will be of the order of 4πλ i . Compared to magnetospheric scales, this is surprisingly small, i.e., a small fraction of the planet's radius. The slower growing, larger-scale modes can subsume the small-scale modes leading to an apparent inverse cascade. Eventually the KH waves can reach global scale with average KH wavelengths ∼ 5 R E on the flanks of Earth's magnetopause boundary (Otto and Fairfield, 2000).
In this paper we consider very simple configurations that promote an inverse cascade. But it is important to note that, e.g., density inhomogeneity can lead to secondary Rayleigh-Taylor instability that can disrupt vortex formation (Matsumoto and Hoshino, 2004). Faganello et al. (2008) showed that in the case where the magnetic field is strictly perpendicular to the flow, density inhomogeneity can lead to both inverse and direct cascades dominating the dynamics, though when an in-plane magnetic field component is added, the secondary instabilities can be stabilized (Cowee et al., 2010;Tenerani et al., 2010). For magnetospheric systems, the high latitude region is typically KH stable and the inverse cascade stops with the vortex scale is comparable to the height of the unstable region (Takagi et al., 2006).
One of the challenges for global-scale simulations is resolving the thin magnetopause boundary. With insufficient resolution and/or excessive numerical diffusion, the fastest growing mode could approach global-scale and the growth rates could be slower than the solar wind advection time scale, prohibiting the onset of KH waves. By using higher order numerics, the Lyon-Fedder-Mobarry model (and later incarnation as the Grid Agnostic MHD for Extended Research Applications or GAMERA simulation (Zhang et al., 2019)) has generated KH modes in both terrestrial (Merkin et al., 2013)  simulations (Zhang et al., 2018). Figure 1 shows a terrestrial GAMERA simulation with the KH waves reaching global-scale on the dawn and dusk flanks (Sorathia et al., 2020). Despite the success of global-scale simulations to model the KH instability, a three-dimensional, meso-scale simulation is required to capture the kinetic-scale processes to understand the interaction between the spectrum of modes that can be present. Typically, a mesoscale (few planetary radii) simulation will capture the inverse cascade and saturate at the m 1 mode that is fixed by the spatial domain. Figure 3 illustrates the inverse cascade from a twodimensional hybrid simulation . The cascade from an m 20 (short wavelength) mode is apparent, with final state of the system in the m 1 mode (longest wavelength supported by the simulation domain). An unresolved issue is the extent to which the slower-growing, smaller-m modes grow independently and subsume the larger m modes, or whether there is a true inverse cascade in which information is transferred from one scale to another. In section 8 we note that in the presence of heavy ions, the protons can separately exhibit a forward cascade while the heavy ions participate in an inverse cascade. This is a topic for more detailed investigation in future studies.

MAGNETIC FIELD TOPOLOGY
Magnetic reconnection plays a critical role in the Kelvin-Helmholtz instability and requires a three dimensional treatment. We start with a discussion of strong guide field reconnection. The standard two-dimensional model (e.g., Sweet and Parker and/or Petschek) of a simple current sheet formed by strictly anti-parallel magnetic fields is misleading because these models neglect the more likely existence of a guide magnetic field in the ignorable, out-of-plane direction ( Figure 4). Magnetic flux conservation in the two-dimensional configuration requires a constant out-of-plane electric field. But note that in the presence of an out-of-plane guide field, the electric field has a parallel component to the magnetic field in a small region of space near the X-line. Schindler et al. (1988) showed that magnetic reconnection occurs if and only if the integral of the parallel electric field is nonzero on a measurable set of field lines in the diffusion region. The parallel potential also gives the rate of reconnected flux . Auroral acceleration in planetary magnetospheres occurs close to the planet on (highly) dipolar magnetic field lines. Extreme guide field reconnection can operate in this case and Seyler (1990) FIGURE 3 | Growth of KH, showing inverse cascade .
FIGURE 4 | Schematic of extreme guide field reconnection and its relation to discrete auroral arc structure (Chaston et al., 2015). Reconnection involves the horizontal magnetic field components.
Frontiers in Astronomy and Space Sciences | www.frontiersin.org December 2021 | Volume 8 | Article 801824 and Otto and Birk (1993) have suggested that thin, discrete aurora arcs are governed by magnetic reconnection with small magnetic shear. Chaston et al. (2015) recently presented the first experimental evidence of the role of extreme guide field reconnection in the auroral acceleration region. Figure 4 summarizes the time sequence of the auroral reconnection process where red and magenta field lines reconnect, giving rise to islands of vorticity in thin auroral arcs (Chaston et al., 2015). Similar reconnection processes are expected to occur in KH vortices. The magnetic field structure associated with KH vortices leads to intense currents due to the twisting of the magnetic field. The net result is driven reconnection which can lead to magnetic field line threading of the magnetopause boundary. In twodimensions, the vortex structure can be defined by magnetic barriers. Figure 5 shows the magnitude of the in-plane component 1 of the magnetic field in a 2-D hybrid simulation. As the surface wave develops, magnetic flux accumulates along the spine and vortex boundary. In the limiting case, where v sh is the flow shear. These magnetic barriers are consistent with the flow deflection that is required for a surface wave. The ongoing vortical motion can bring oppositely directed magnetic fields in close proximity, resulting in reconnection of the in-plane magnetic field components (Nykyri and Otto, 2001).
While a simple 2-D picture of KH illustrates general KH properties, many interesting aspects are revealed in three dimensions. When the antiparallel components of the magnetic field in the presence of flow shear are considered in a 2-D plane, KH and reconnection are mutually exclusive. If the flow shear is sub-Alfvénic, spontaneous reconnection can operate (and KH is stabilized), while if the flow is super-Alfvénic, KH is destabilized and spontaneous reconnection is suppressed as information transport along the magnetic field is truncated (La Belle-Hamer et al., 1995;Faganello et al., 2010). We note that during the KH growth, driven reconnection can occur as current layers are pinched during vortex rollup with a rate that is independent of prescribed mechanism (Knoll and Chacón, 2002;Nakamura and Fujimoto, 2006). In three dimensions, on the other hand, KH and reconnection can interact. Ma et al. (2014a) showed that for an initial KH driver, reconnection can be triggered along the spine region as the current layer thins and intensifies. On the other hand, Ma et al. (2014b) showed that reconnection flows for a preexisting current layer at the magnetopause boundary can trigger the KH instability. Ultimately, the twisting of the magnetic field can lead to strong guide field reconnection in pairs of reconnection sites (Otto, 2007;Faganello et al., 2012Faganello et al., , 2014. Ma et al. (2017) showed that a KH-active magnetopause boundary, the reconnection sites form in the midlatitudes and asynchronous reconnection can lead to open flux generation. The double reconnection process can also effectively transport plasma across the magnetopause boundary.
During the non-linear rollup of KH vortices, multiple current sheets (and associated magnetic barriers and/or filaments) are folded together. Therefore, a transect of the vortex region yields a highly variable magnetic field with many bipolar signatures, typically in component pairs. Delamere et al. (2013) used these bipolar signatures to identify potential KH active regions near Saturn's magnetopause boundary, though any vortical flow associated with, e.g., radial transport can produce similar bipolar signatures (Stauffer et al., 2019).
In three-dimensional hybrid simulations, the complexity of the magnetic field topology is further enhanced during the inverse cascade. A key aspect of hybrid simulations is that the finite particle statistics provide the "stochastic" seed for surface wave formation. The phase of the surface wave at a given point along the magnetic field is therefore random. Provided that the simulation domain along the unperturbed magnetic field is sufficiently long (e.g., perpendicular to parallel spatial scale proportional to the Alfvén Mach number of the flow shear), then the phase of the surface waves can be a function of distance along the magnetic field. In other words, the system supports k ∼ k ⊥ where k ⊥ is the KH wavenumber of the fastest growing mode. As a result, twisted nodes form between out-of-phase surface waves, promoting strong guide field reconnection. At each stage of the inverse cascade, twisted node formation promotes reconnection. As a result, multiple reconnection sites can be identified on an instantaneous field line trace, adding further complexity to the magnetic field topology. Figure 6 shows field line tracing from a 3-D hybrid simulation. The image slices and field line color show the magnitude of the parallel electric field, E ′ E /(v A B 0 ), or the expected strong guide field reconnection rate for the case where the electron pressure term was neglected (i.e., T e 0). In a high Alfvén Mach number flow, the in-plane magnetic field perturbations can approach B 0 , so the normalization uses B 0 as an upper limit (lower limit for the reconnection rate). The sample magnetic field lines originate from a localized region on the magnetosheath side of the boundary. The field lines trace to a region along the spine in the second slice (y ∼ 100 λ i ) in a manner similar to that described by Otto and Birk (1993) in a current striation model for understanding formation of long thin auroral structures. Note the strong parallel electric field in this region (E ′ ∼ 0.2). Subsequent field lines continue to be altered by parallel electric fields until the mapping is spread over a significant portion of the KH vortex. In addition, the field lines thread the magnetopause boundary (defined by 50% mixing); therefore, momentum transport (i.e., Maxwell stress) across the magnetopause boundary occurs Delamere et al., 2018). Future simulations should compare the mapping with field-line tying boundary conditions (i.e., to mimic the ionospheric boundaries) following Ma et al. (2017).

MOMENTUM TRANSPORT
Momentum transport at the magnetopause boundary occurs as a result of the KH instability. The classic Chapman-Ferraro magnetopause boundary is the idealized case where the normal component of the magnetic field is zero and F T ⃡ · da 0, where T ⃡ is the Maxwell stress tensor with units of momentum flux. Due to the magnetic field line threading of the magnetopause boundary, momentum can be transported via Maxwell stresses (BB/μo). Reynolds stresses (ρuu) can also be generated as a result of the flow shear and the vortical motion.
Assuming an initial pressure balanced configuration, an estimate of the momentum transfer rates will be made by considering the steady flux conservative form of the momentum equation, i.e., This equation will be integrated over any arbitrary volume and rewritten as a surface integral using the divergence theorem (i.e., ∇ · T ⃡ dV T ⃡ · da). If +t is the tangential sheath flow (tailward) direction and +n is the direction normal to the magnetopause boundary, then the sheath flow is modified (i.e., reduced) due to Maxwell shear stresses at the magnetopause boundary when BB ≈ − B t B ntn ≠ 0. The Reynolds stress are (ρu t u n )tn and account for plasma source terms due to the direct transport of plasma across the magnetopause boundary. Figure 7 shows an example of the (normalized) Maxwell and Reynolds stresses at the magnetopause boundary in a 3-D hybrid simulation where, e.g., Txz is x momentum transported in the z direction. Using the formalism of Miura (1984) and we can evaluate the "anomalous viscosity", v ano ρ −1 0 〈B t B n /μ o − ρv t v n 〉 (d〈v t 〉/dn) −1 , to formally quantify KH contribution to the socalled "viscous-like" interaction. Delamere et al. (2018) estimated the flow reduction in Saturn's dawnside magnetosheath due to Maxwell (T M ) and Reynolds stresses (T R ) from 3-D hybrid simulations. Figure 7 shows the midplane of the simulation (i.e., magnetopause boundary) for both T M and T R where the color scale indicates relative momentum fluxes. The upper limit for combined stress is roughly 4 × 10 -13 N/m 2 and this yields a flow reduction that is consistent with observations . Again, it should be emphasized that significant magnetic field line threading of the magnetopause boundary occurs and in general the Maxwell stresses (T M xz ) were found to be larger than the Reynolds stresses.

PARTICLE TRANSPORT
There are two mechanisms that facilitate mass transport. In the magnetohydrodynamic limit, magnetic reconnection is required, involving the flow-aligned components in both 2-D and 3-D cases (Fadanelli et al., 2018;Sisti et al., 2019). Hybrid simulations, on the other hand, show both magnetic reconnection operating as well as diffusive transport due to large gyroradius effects and waveparticle interactions (Johnson and Cheng, 1997;Chen, 1999;Chaston et al., 2009). In a 2-D comparison of MHD + test particles vs. hybrid simulations with periodic boundary conditions, Ma et al. (2019) showed that both simulations produced similar transport rates, but that particle mixing by finite gyroradius is the dominant process in the hybrid simulations. The derived diffusion coefficients are typically ∼ 10 9 to 1010 m 2 /s (Cowee et al., 2009(Cowee et al., , 2010Delamere et al., 2011), andCowee et al. (2009) indicated that the time dependent diffusion coefficients can describe a non-classical "superdiffusion" process.
Hybrid simulations allow for a simple calculation of the diffusion coefficient. Initially, the particles are labeled as either magnetospheric or magnetosheath particles. During the evolution of the KH instability, the ions are mixed on the grid with fully mixed grid cells containing 50% of both populations, or a mixing of 0.5. The width of the mixed cell region is L mix (t) where values above, e.g., 0.1 are considered mixed. The diffusion coefficient is D (t) dL mix (t) 2 / dt. Figure 8 compares mixing/diffusion in 2-D and 3-D hybrid simulations. In both cases the mixing is very strong in the KH active region, but the 3-D simulations tend to produce a narrower diffused region (boundary layer) due to the magnetic tension generated by out-of-phase surface waves during the inverse cascade. Presumably when the KH modes reach global scale and periodic boundary conditions are no longer valid, the system will evolve with coherent vortical structure as seen in the global simulations.
Boundary layer formation can occur rapidly when growth rates exceed the solar wind/magnetospheric advection time. For an initial boundary width of the order of the ion inertial length, the growth rates can exceed c 0.01 s −1 , or a growth time scale of minutes . Eventually the boundary layer will thin such that the onset of KH waves is episodic. Burkholder et al. (2019) found an occurrence frequency ∼ 18% of active KH waves in the 10-12 LT sector of Saturn's magnetopause boundary where KH activity is expected to be a maximum due to the strong flow shears. During this episodic activity, particle transport should be observed, particularly in the superthermal ion population. Indeed, Sergis et al. (2013) reported on water group (i.e., originating from Enceladus in Saturn's inner magnetosphere) "islands" in the Saturn's magnetosheath for particles > few keV. The observed pitch angle distribution of the W + ions was typically peaked away from 90°with a substantial field-aligned component (i.e., streaming distribution). While these energetic ions have large ion gyroradii and can likely escape across the magnetopause boundary, the transport could be modulated by KH waves. Similar distributions of iogenic ions have been observed by the Juno spacecraft at Jupiter (Mauk et al., 2019). At Earth, energetic oxygen ions tend to be entrained along the dayside magnetopause boundary during stable magnetospheric conditions and are less likely to be scattered compared to protons (Mauk et al., 2019). While much work remains regarding the dependence on mass/charge on energetic particle escape at Earth, Jupiter, and Saturn, we suggest here that KH waves could modulate the escape process, accounting for the episodic water group island observed in Saturn's dayside sheath. As a preliminary test of KH-related transport of energetic ions, we conducted a proton-only, 3-D hybrid simulation that contained a small population (e.g., 1% of total mass density) of energetic protons (4v th ) on the magnetospheric side of the flow shear boundary. During the growth of the KH waves, these particles were found to be transported across the boundary (i.e., to regions dominated by magnetosheath plasma). The velocity distribution was often peaked well away from 90°pitch angle. Figure 9 shows an example of such a streaming distribution with a peak near 130°, similar to the distributions reported by Sergis et al. (2013). We note that streaming distributions have also been measured by the magnetospheric multiscale mission at Earth. Vernisse et al. (2016) suggested that the streaming particles were generated by reconnection induced by KH instabilities. The Themis and Cluster missions have also observed streaming ion distributions far away from the vortex location (Bavassano Cattaneo et al., 2010;Faganello et al., 2014). Figure 10 shows a sample field line trace taken through the KH active region. The top to bottom is the ion pitch angle distribution of the superthermal magnetospheric ions, the ion mixing, the normalized x component of flow, the normalized perpendicular magnetic field perturbations, the parallel electric field, and the integrated parallel electric field (field-aligned potential). In this example, the magnetic field samples a wide range of mixing values, plasma flows. We note that the pitch angle distribution is often not peaked near 90°as would be expected for an unperturbed Maxwellian distribution. The modification of the distribution function as a function of energy is a topic for future study; however, initial studies suggest that transport can be facilitated and/or modulated by KH activity.

ENERGY TRANSPORT AND HEATING
KH waves radiate compressional modes into the magnetospheric cavity. Pu and Kivelson (1983) and Delamere et al. (2011) showed that roughly 1% of the incident solar wind power can be radiated into the magnetosphere. As a result, magnetospheric resonant cavity modes can be excited and/or plasma heating can occur. The empirical estimate for the plasma heating rate at Saturn is 30-400 GW (Bagenal and Delamere, 2011), which compares favorably with Poynting flux in generated in hybrid simulations multiplied by the expected KH-active area of the dayside magnetopause boundary . The internal energy transport via the magnetosonic fast mode can be comparable to or even shorter than the solar wind advection time past the magnetosphere (e.g., hours for the giant magnetospheres); therefore, prompt KH modulation of magnetospheric processes should be expected. Though, the details of energy flow partitioning in the magnetospheric system remains unresolved.
Within the KH vortex, on the other hand, significant plasma heating can be expected based on models for the non-linear interaction between counter-propagating Alfvén waves (Burkholder et al., 2020a). Delamere et al. (2021) examined the heating in 3-D KH simulations and found that indeed the ion heating is consistent with the turbulent heating mechanism (Saur, 2004). Assuming a turbulent spectrum of magnetic field fluctuations (δB ⊥ ), the heating rate density for the kinetic Alfvén wave (KAW) is where ρ is the mass density and ρ i is the thermal ion gyroradius. Figure 11 shows the heating rate density for a 3-D hybrid simulation, showing the localization of regions with significant heating. The contours show the magnitude of the in-plane magnetic field. Future studies should investigate, in detail, the interaction of ions with kinetic Alfvén waves to fully understand the heating process.

HEAVY ION EFFECTS
Magnetospheric plasmas are rarely homogeneous and often contain multiple ion species. At Earth, ionospheric outflow can introduce heavy ions into the system (i.e., O + ) (Welling et al., 2015). Bouhram et al. (2005) showed that heavy ions are the dominant contributor to the mass density 30% of the time on the dusk flank and 5% of the time on the dawn flank. In the giant planet magnetospheres, heavy ions are sourced from the moons and play a significant role in magnetospheric dynamics. Magnetospheric heavy ions can also modify the growth rates of the KH instability and at kinetic-scale can significantly modify the structure of the KH modes.  showed, using 2-D hybrid simulations, that: 1) growth rates tend to increase when protons (magnetosheath) are separated from the heavies (magnetosphere); 2) growth rates tend to decrease with equal mixtures of protons and heavies; and 3) decrease when only heavy ions are present on both sides of the boundary. In the latter case, it was suggested that inhibited growth was due to suppression of wavelengths less than the ion inertial length. The dawn/dusk asymmetry of heavy ions at Earth could have implications for KH-related plasma transport.
In the limiting case of magnetospheric heavy ions separated from the magnetosheath protons, the initial onset of the KH modes is distinctly different. Figure 12 shows the results from a 2-D hybrid simulation during the onset of the KH instability for parameters based on Saturn's magnetopause boundary environment (i.e., heavy ions 17 AMU for water group ions). The top region is the magnetosheath containing protons and the bottom region is the magnetosphere containing heavy ions. From left to right and top to bottom the figure shows 1) the magnetic field pressure, 2) plasma pressure, 3) total ion mass density, 4) proton density, 5) heavy ion density, and 6) normalized local entropy. Perhaps the most striking result is the absence of heavy ions in the magnetosheath. The onset of an m ∼ 6-7 mode is clearly present, but the heavy ions do not participate in the fluid-scale wave motion because the wavelengths of the m > 6 modes are less than the heavy ion inertial length (e.g., cannot generate the compressional magnetosonic fast mode). Thus, the heavies cannot participate in the wave motion. On the other hand, the protons are affected by the surface wave and move in the vortical flows. The appearance of protons on the magnetospheric side of the boundary is not balanced by heavy ions in the sheath, leaving voids in plasma pressure. The local entropy (p/ρ c ) increases because of the apparent dominance of the inverse relation with plasma mass density (ρ). In terms of long-term wave evolution, the heavy ions exhibited the expected inverse cascade while the protons showed a forward cascade, saturating at m 5 ( Figure 13). These simulations are highly idealized but illustrate the potential importance of heavy ions in modifying the KH structure at kinetic scale. Future studies should examine the transfer entropy from conditional mutual information theory to further understand these forward and inverse cascades.
In a 2-D study of density asymmetry with heavy ions, Meier (2012) found that mixing was reduced by increasing the mass density across the magnetopause boundary. The heavy magnetospheric ion (16 amu) fraction was fixed at 25% (with a mixture of protons) while magnetosheath to magnetosphere number density ratio was varied between 2 and 8 (the corresponding mass density ratio varied between 0.5 and 2). Pressure balance was achieved by modifying the magnetospheric proton temperature accordingly. In addition to reduced mixing (presumably due to large gyroradius effects in the non-linear rollup phase), the vortices developed on the low density (magnetospheric) side of the boundary. Thus, the density asymmetry may expedite the formation of mixed boundary layers in the giant magnetosphere as these asymmetric density conditions are often found at Jupiter and Saturn.
Heavy ions can affect mass transport at the magnetopause due to wave-particle interactions. While we have not explicitly investigated these effects in KH hybrid simulations, we note that transport coefficients can be derived following standard methods from the FIGURE 10 | Magnetic field line trace through KH active region in hybrid simulation, showing various properties. From top to bottom is "superthermal" ions, plasma mixing, plasma flow, in-plane magnetic field normalized to the background magnetic field, parallel electric field, and field-aligned potential.
Frontiers in Astronomy and Space Sciences | www.frontiersin.org December 2021 | Volume 8 | Article 801824 drift-kinetic or the gyrokinetic equation (Johnson and Cheng, 1997;Chen, 1999;Chaston et al., 2009). In particular, some of the dominant mechanisms include contributions from Landau damping, magnetic field gradient drift, and transit time damping associated with KAW-related magnetic bottles. The presence of heavy ions affects the diffusion coefficients in two ways. First, increased ion mass leads to larger parallel electric fields. This effect is necessary because the ion polarization drift is larger for heavy ions, leading to increased charge separation and thus the parallel electric field required to maintain quasineutrality. Similarly, the ion density fluctuations associated with the polarization drift enhances the magnetic field compressions (particularly in a high β plasma), enhancing the transit time damping. Second, the heavy ions reduce the Alfvén speed which increases the fraction of resonant particles in the wave-particle interaction. Future studies should compare hybrid simulation results with theoretical expectations for the effects on heavy ions on transport.
FIGURE 11 | Turbulent heating rate density (Delamere et al., 2021). Contours indicate the magnitude of the in-plane magnetic field. While electron kinetic effects cannot be addressed directly in hybrid simulations, it is possible to study the effects of parallel electric fields using electron test particle simulations. The parallel electric fields have been evaluated in both Rayleigh-Taylor (Stauffer et al., 2019) and KH simulations and are found to be independent of the specified resistivity. Estimates for the magnitude of the parallel electric fields (Section 4) based on local conditions in the simulations show that the reconnection rates approach or exceed 0.1, consistent with expectation for driven fast reconnection.
Using an ensemble of field line traces (i.e., similar to Figure 6), we conducted electron test particle simulations for parallel motion, i.e., m dv dt qE We initialized a 20 eV Maxwellian distribution on the field line ensemble. Figure 14 shows the phase space (v vs. y) at a snapshot in time for a sample field line, showing field-aligned streaming (acceleration) as well as trapped populations. Assuming that the evolution of the non-linear KH vortex modifies the parallel electric field configuration in time, we conducted an ensemble average to determine the likely steady state electron velocity distribution. For Saturn-like conditions, the final distribution was consistent with a 120 eV Maxwellian. Curiously, this is consistent with the thermal electron temperature in Saturn's outer magnetosphere (Neupane, 2021), but further study would be required to understand how internal radial transport could affect the thermal electron temperature. Figure 15 shows the initial distribution (red) with a 20 eV Maxwellian fit and the final distribution (blue) with a 120 eV Maxwellian fit. A snapshot of any given field line would likely reveal bidirectional electron beams. Such beams are commonly observed in KH active regions of the terrestrial magnetopause boundary (Burkholder et al., 2020b). Nykyri et al. (2021) suggested that, among other mechanisms, magnetic reconnection could be a candidate for explaining the > 100 eV counter-streaming electrons observed by the Magnetosphere Multiscale mission. An example topic for further study is the role of oblique whistler waves in electron energization. Finally, the addition of heavy ions would enhance the parallel electric fields as noted in Section 8, leading to more effective electron energization.
FIGURE 13 | KH wave mode cascade with heavy ions for the simulation shown in Figure 12. The heavy ions exhibit an inverse cascade to m 1 while the protons show a forward cascade to m ∼ 5.
FIGURE 14 | Electron test particle simulations using the sample magnetic field line trace. The distribution in Figure 15 is an ensemble average of many sample field lines taken from the KH active region.
Frontiers in Astronomy and Space Sciences | www.frontiersin.org December 2021 | Volume 8 | Article 801824 Considerable progress has been made with regard to ion kineticscale processes associated with the Kelvin-Helmholtz instability.
Here we summarize the key advances in understanding and on directions for future studies.
• The growth of KH modes in our simulations can be described as a inverse cascade. For typical magnetopause boundary thickness (few ion inertial lengths), the fastest growing modes are relatively small scale (less than planetary radii) and, under certain conditions, could evolve to global scales (few planetary radii) as the structures advect down the tail flanks. • The inverse cascade results in smaller scale structures being subsumed by larger scale structure, resulting in complex magnetic field topologies with a patchy network of strong guide field reconnection sites. • Maxwell stresses associated with magnetic field line threading of the magnetopause boundary lead to substantial momentum transport and contribute to viscous-like tangential drag at the boundary. • Plasma transport is facilitated by the patchy network of reconnection sites and diffusive transport is possible due, in part, to finite gyroradius effects. Energetic particles can be transported across the boundary with a strong field-aligned component to the velocity distribution. • Turbulent magnetic field spectra can lead to substantial ion heating rates due to the non-linear interaction between counter propagating Alfvén waves. Future work should focus on the heating mechanism due to the interaction between particles and the kinetic Alfvén waves (Johnson and Cheng, 2001).
• Heavy ions can alter the growth rates, transport rates, and modify ion kinetic-scale structure due to competing ion scales during KH wave evolution. In the limiting case, heavy ions contained entirely on, e.g., the magnetospheric side of the boundary will be unaffected by KH wavelengths less than the heavy ion inertial length. • Bidirectional electron beams can be expected from the parallel electric fields associated with the strong guide field reconnection.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
PD contributed to all aspects of this article. Author NB provided the field line tracing analysis. XM contributed to the general discussion of the interaction between reconnection of KH modes in three dimensions. JJ contributed to the analysis of heavy ion effects and turbulent heating.