In a cold plasma, basic equations of an ideal MHD plasma are ρ ∂ ∂ t + V ⋅ ∇ V − ν V = 1 μ 0 ∇ × B × B , (3) ∂ B ∂ t = ∇ × V × B , (4) where ν is a collisional frequency that is introduced to damp waves propagating outside the region of interest, which effectively imposes outgoing boundary conditions. It should be noted that collisional effects play no role in the stability of the primary or secondary KH instabilities that we analyze in the rest of this article.

We assume that a field variable consists of background equilibrium (0) and small perturbation (1) components (B = B ₀ + B ₁, ρ = ρ ₀ + ρ ₁, and V = V ₀ + V ₁), and a shear flow ( V 0 ( x ) = V 0 ( x ) y ̂ ) and a uniform background magnetic field ( B 0 = B 0 z ̂ ) lie in the y- and z-directions, respectively. Then, the perturbed quantities can be Fourier analyzed in the y- and z-directions (∂/∂y → ik _y and ∂/∂z → ik _∥, where k _y and k _∥ are wavenumbers in the y and field-aligned (z) directions). Thus, the linearized MHD wave equations are ρ 0 ∂ ∂ t + i k y V 0 V 1 x + ν V 1 x = B 0 μ 0 i k ∥ B 1 x − ∂ B 1 z ∂ x , (5) ρ 0 ∂ ∂ t + i k y V 0 V 1 y + ν V 1 y = B 0 μ 0 i k ∥ B 1 y − i k y B 1 z − ρ 0 V 1 x ∂ V 0 ∂ x , (6) ∂ ∂ t + i k y V 0 B 1 x = i k ∥ B 0 V 1 x , (7) ∂ ∂ t + i k y V 0 B 1 y = i k ∥ B 0 V 1 y + B 1 x ∂ V 0 ∂ x , (8) ∂ ∂ t + i k y V 0 B 1 z = − i k y B 0 V 1 y − B 0 ∂ V 1 x ∂ x . (9)

In Section 3, we solve the spectrum of eigenmodes of these equations in slab geometry, while in Section 4, we solve these equations using a finite-difference time-domain method.

To proceed with the spectral analysis, we define an auxiliary set of variables including the fluid displacement (ξ) V 1 ≡ ∂ ∂ t + V 0 ⋅ ∇ ξ , (10) the total pressure perturbation (p), and compressibility (ψ), p ≡ B 0 B 1 z , (11) ψ ≡ ∇ ⋅ V 1 . (12)

Taking the Fourier transform in time ( ∂ ∂ t → − i ω ) and ignoring the collision term (ν → 0), Equations 5–9 become μ 0 ρ 0 ω ̃ 2 ξ x = − i B 0 k ∥ B 1 x + ∂ p ∂ x , (13) μ 0 ρ 0 ω ̃ 2 ξ y = − i ω ̃ μ 0 ρ 0 ξ x ∂ V 0 ∂ x − i B 0 k ∥ B 1 y + i k y p , (14) ω ̃ B 1 x = i ω ̃ B 0 k ∥ ξ x , (15) ω ̃ B 1 y = i ω ̃ B 0 k ∥ ξ y + i B 1 x ∂ V 0 ∂ x , (16) ω ̃ B 1 z = − i B 0 ψ , (17) where ω ̃ = ω − k y V 0 .

Then, Equations 11–17 can be reduced to two coupled first-order differential equations, d p d x = μ 0 ρ 0 ω ̃ 2 − k ‖ 2 V A 2 ξ x , (18) μ 0 ρ 0 d ξ x d x = − ω ̃ 2 − k y 2 V A 2 + k ‖ 2 V A 2 ω ̃ 2 − k ‖ 2 V A 2 p B 0 2 . (19)

We solve Equations 18 and 19 to analyze the eigenmode frequency in Section 3.

Using Equations 31 and 32, we calculate the eigenfrequency (ω), growth rate (γ), and amplitude ratio between magnetic compressional component (A _p ) for various widths of the shear transition layer, k _y d = 0, 0.25, 0.75, and 2.0, as shown in Figure 2. For these plots, the plasma densities in region I and region III are assumed to be N _0I = 5 × 10⁶/m³ and N _0III = 5 × 10⁵/m³, and the background magnetic field strength is B ₀ = 25nT. We also specify an angle of propagation (ϕ) with respect to the ambient magnetic field, ϕ = tan⁻¹(k _y /k _‖) = 80° and k y 2 + k ‖ 2 = π / ( 2 R E ) . For complete stability analysis, this angle would be varied to determine the maximum growth rate for a given Mach number. The upper panels of Figure 2 are the calculated real (black and red) and imaginary (blue, growth rate γ) frequencies as functions of the Alfvén Mach number (M _A ≡ V _0I /V _AI ), and the lower panels plot the amplitude ratio A _p of unstable wave modes.

In the absence of the shear transition layer (k _y d = 0) as shown in Figure 2A, forward and backward propagating fast waves, which have positive and negative frequencies at M _A = 0, occur when M _A is small (Taroyan and Erdélyi, 2002). These waves are stable until M _A reaches the threshold of the KH instability, M A s = tan − 1 ( ϕ ) 2 ( 1 + V A I / V A I I I ) . For M _A > M _As marked as a gray-shaded region in Figure 2A, the waves develop a complex frequency and become unstable. For the given range of M _A , ω and γ increase linearly with M _A . Because the characteristics of this wave mode are the same as the typical KH waves (e.g., Johnson et al., 2014), this wave corresponds to primary KH waves (hereafter PKHW). In this figure, we also found a shear Alfvén wave mode at ω = ω _AI = k _‖ V _AI . The fast and shear Alfvén waves cross each other near M _A ∼ 0.45, but the coupling of the two wave modes does not occur.

Introducing a finite width of the shear transition layer significantly changes the wave dispersion relations. In Figure 2, the PKHWs also occur for the cases of k _y d ≠ 0. The M _A threshold decreases from 0.83 for k _y d = 0 to 0.35 for k _y d = 2.0. Overall, the wave frequency ω decreases, while the growth rate γ increases as k _y d increases. For example, for M _A = 0.85, ω/ω _AI = (4.35, 3.68, 2.95, 2.48) and γ/ω _AI = (0.35, 1.5, 1.85, 1.89) when k _y d = (0.0, 0.25, 0.75, 2.0). Thus, when a shear transition layer is included, lower frequency PKHWs are excited with a stronger growth rate and lower M _A threshold.

For k _y d = 0.25 in Figure 2B, coupling between the backward propagating fast and shear Alfvén waves occurs near ω/ω _AI ∼ 1, and unstable waves also appear for 0.355 ≤ M _A ≤ 0.47 (shaded yellow in Figure 2B). These waves correspond to the secondary KH waves (hereafter SKHW) (Turkakin et al., 2013). In this case, the SKHWs are clearly separated from the PKHWs and have lower ω, lower γ, and lower M _A threshold than the PKHWs. The compressional amplitude ratio (A _p ) in the lower panel shows significant differences between PKHWs and SKHWs; A _p ≪ 1 for the PKHWs and A _p ∼ 1 for the SKHWs. Therefore, for SKHWs, the amplitude of the instability is similar at both the V ₀ and V _A interfaces, indicating a spreading of wave power over a more extended region, while the PKHWs are localized about the V ₀ interface. When the Mach number is low, it is expected that only the SKHWs would be excited.

When the V _A interface is further away from the V ₀ interface (k _y d = 0.75), as shown in Figure 2C, the PKHW and SKHW modes merge near M _A ∼ 0.47. Although ω and γ monotonically increase as a function of M _A , the KH waves have similar behavior to the SKHW (A _p ∼ 1 and ω ∼ ω _AI ) at smaller M _A and the PKHWs (A _p < 1 and ω ≫ ω _AI ) at larger M _A . Thus, the waves may still be divided into the semi-SKHW marked as a light yellow-shaded region and PKHW marked as a gray-shaded region in Figure 2C.

For k _y d = 2, as shown in Figure 2D, only a single unstable wave mode corresponding to the PKHWs occurs localized at the V ₀ interface. The V _A profile can be treated as a constant at the V ₀ interface and the M _A threshold becomes M _As ∼ 2 tan⁻¹(ϕ) = 0.359. The threshold occurs near ω/ω _AI ∼ 1; thus, the wave frequencies are always higher than ω _AI .

It is also useful to examine how ω and A _p depend on M _A and k _y d. Figure 3A,B shows contour plots of ω normalized to 1) ω _KH0 and 2) ω _AI , respectively. In this figure, two wave modes are clearly organized by ranges of M _A ; the PKHW for M _A > 0.47 and the SKHW for 0.355 ≤ M _A ≤ 0.47. Red and magenta lines in Figure 3A represent the M _A threshold for the PKHW and SKHW, respectively. The M _A threshold of the PKHWs decreases and the upper M _A limit of the SKHWs increases as k _y d increases. The thresholds merge near k _y d ∼ 0.534 and M _A ∼ 0.47. Thus, for M _A > 0.47, a single wave mode appears (see Figure 2C); however, wave characteristics at lower and higher M _A are significantly different.

The PKHWs show that all parameters (ω/ω _KH0 , ω/ω _AI , and A _p ) have a strong dependence on k _y d, and they decrease as k _y d increases. For most M _A , ω/ω _KH0 ∼ 1 and 1 < ω/ω _KH0 < 2. Because both ω and ω _KH0 increase proportionally to M _A , ω/ω _KH0 has less dependence on M _A . However, because ω _AI does not depend on k _y and ω increases as M _A increases, ω/ω _AI depends on both M _A and k _y d. For the given conditions, ω/ω _AI is maximized when k _y d is small and M _A is large. Figure 3C shows A _p < 1 for M _A ≥ 0.47, except k _y d → 0. Thus, it shows that the PKHWs are almost always dominant at the V ₀ interface. For k _y d → 0, a strong amplitude of the pressure term occurs at the secondary interface. However, this increase in A _p is not an indicator of a separate instability, but rather it simply indicates that the decay of the wave power from the V ₀ interface to the V _A interface reduced as the shear layer vanishes.

On the other hand, the eigenmode frequency of the SKHWs is comparable to ω _KH0 and ω _AI (0.9 ≤ ω/ω _KH0(AI) ≤ 1.2) in the entire range of k _y d and M _A because this wave mode appears due to the coupling between shear Alfvén mode and the fast compressional waves (thus, ω _KH0(AI) ∼ ω _AI ). For the entire range of k _y d, A _p is always close to or even higher than 1. These results suggest that the KH instability occurs at both the V ₀ and V _A interfaces with almost the same amplitude even though the interfaces are well separated.

The eigenmode calculations can be summarized as follows: the PKHWs are localized at the V ₀ interface having a higher frequency than ω _A in the magnetosheath for faster shear flow velocity, while the SKHWs can be detected at both the V ₀ and V _A interfaces with similar wave frequency to ω _AI in the magnetosheath for slower shear flow velocity.

We first examine the KH waves in a plasma where V _A does not vary in space. In this simulation, a hyperbolic tangent V ₀ profile along with constant V _A was adopted in the wave code: V 0 x = V 0 I 2 1 − tanh x a , (33) where V _0I is the flow velocity in region I, and this profile characterizes the V ₀ discontinuity in a scale length a, as shown in Figure 1B. One of the primary differences between the background profile used in the time-dependent analysis, and the previously discussed eigenmode analysis (a → 0) is the fact that the discontinuous profile has been smoothed.

We assume that the length of the simulation box is L _x ∼ 45/k _y . Since the KH surface wave is expected to not fully decay by the time it reaches the edge of the simulation domain in the x-direction, we add an absorption layer near the boundary (∼30/k _y ) in the simulation box to prevent reflection. An initial perturbation is launched as a compressional component of V _1x at the source location (k _y x _source ∼ − 7.5) in region I (i.e., magnetosheath). This source is assumed to have a narrow spatial width (k _y δ _source = 0.093) and to include broadband frequencies, V 1 x ( x , t ) = exp − 1.5 t 2 t K H 2 exp − ( x − x s o u r c e ) 2 δ s o u r c e 2 , where k _KH = 2π/ω _KH0 . The simulation is run from t = 0 to t = 5.6t _KH , and all components of B ₁ and V ₁ at each time step are stored during the simulation run time. The background densities in the magnetosheath (region I) and the magnetosphere (region III), the background magnetic field strength, k _y , and k _z are the same as in the eigenmode analysis of Section 3.

Figure 4 shows the time evolution of the magnetic compressional component (B _1z ) in the x-direction for M _A = 1 and k _y a = 0.025. Two vertical lines represent the source location (k _y a = − 7.5) and the V ₀ interface (k _y x = 0), and thick dashed lines represent Alfvén speed (V _A ). Since the initial wave packet includes broadband frequencies, the wave packet disperses in time and space. Leftward propagating waves reach a strong collisional layer near the boundary (k _y x < − 10.5) and are totally absorbed. Rightward propagating waves reach the V ₀ interface at k _y x = 0 around t/t _KH = 0.7, and they partially reflect from the interface due to a steepened density gradient. The rest of the waves penetrates the V ₀ interface and reach the collisional layer (k _y x > 3.0). Once the magnetic field and velocities are perturbed near the interface, an unstable wave mode begins to grow at around t/t _KH ∼ 1.2. Unlike the initial perturbation, these waves decay in the x-direction rather than propagate. The wave amplitude in Figure 4 saturates at ± 100.

We focus on the surface waves at k _y x = 0 and determine the growth rate, wave frequency, and polarization. Time histories of B _1z and B _1y at x = 0 in Figure 5A rapidly grow in time; thus, the sinusoidal wave form is not clearly seen. However, the wave growth term can be removed from the time histories using the magnetic (U _B ), kinetic energy (U _V ), or total energy (U = U _B + U _V ). We plot U _tot (t) = ∑ _x U(x, t) in the simulation box in Figure 6A. Early in the simulation period (t/t _KH < 1.2) U _tot is quasi-stable; however, once an unstable waves generated, it increases linearly. The wave magnetic field with a constant growth rate γ can be written as B 1 x , t ∼ b 1 x , t exp γ x , t t , (34) and because of the magnetic energy U _B ∝ |B|², the wave growth rate γ in each grid point can be estimated from γ x , t ∼ ∂ ∂ t ln U B x , t . (35)

We also confirmed that γ calculated using either the magnetic (U _B ) or kinetic energies (U _V ) are identical; thus, γ ( x , t ) = ∂ ∂ t ln ( U ( x , t ) ) = ∂ ∂ t ln ( U B ( x , t ) ) = ∂ ∂ t ln ( U V ( x , t ) ) . Furthermore, once the initial wave vanishes near the boundary, only the localized surface waves (such as KH waves) remain in the simulation domain; thus, γ also can be calculated using U _tot : γ t = ∂ ∂ t ln U t o t t . (36)

When a boundary has a finite thickness, the normalized growth rate (Γ ≡ 2aγ/V _0I ) becomes a function of normalized boundary width (2k _y a) (Miura and Pritchett, 1982). To illustrate, the time evolution of Γ(t) is plotted in Figure 6B for k _y d = 0.025 and M _A = 1, and it converges to ∼ 0.0083 . Therefore, for these parameters, the normalized growth rate can be estimated as Γ = 0.0083.

Once the growth rate is determined, the wave components (and polarization) can be obtained from b 1 t ∼ B 1 t / exp γ t − t 0 t > t 0 , (37) where t ₀ is the time at which the wave growth begins. Figures 5C,D show that b _1z and b _1y have clear sinusoidal structures with a single frequency. The wave spectra of b _1z and b _1y in Figures 5E,F confirm that the single peak corresponds the KH wave frequency, ω K H 0 = 1 2 k y V 0 . In this manner, we can determine both the real and imaginary components of the frequency, which can be compared with the eigenmode analysis.

For code validation, we also compared the simulation results with prior analytical results in Miura and Pritchett (1982). Figure 7 shows the growth rate, Γ, as a function of (a) 2k _y a for M _A = 1 and (b) as a function of M _A for 2k _y a = 1. In this figure, the prior analytic results (gray lines) and our simulations (red stars) show excellent agreement with each other. For M _A = 1 in Figure 7A, wave growth only occurs for a limited value of the normalized boundary width 0 < 2k _y a < 1.8 and maximizes near 2k _y a = 0.8. For 2k _y a = 1, the maximum Γ occurs for M _A → 0 and has a value of 0.144, as predicted from Miura and Pritchett (1982). The growth rate decreases as M _A increases and no KH wave arises for M _A > 1.6. Therefore, the new MHD wave code successfully demonstrates KH waves and benchmarking comparisons of the simulations with previous analytical results validate the code accuracy.

In this section, the simulation results include the shear transition layer between the magnetosphere and the magnetosheath, as shown in Figure 1B. In contrast to the results of Section 3, we consider a finite width of the boundary layer. Similar to the V ₀ profile in Equation 33, V _A is assumed to have a hyperbolic tangent profile: V A x = V A I + V A I + V A I I I 2 tanh x − d a . (38)

The two interfaces are separated with width d, although each interface has its own width, a. From the eigenfrequency calculations in Figure 2, we showed that the inclusion of a shear transition layer effectively generates the SKHWs when the shear flow velocity is slow; thus, we ran the simulations for k _y d = 0.25 and 0.75 and M _A < 0.85 to compare with the eigenmode calculation.

We used the time histories of b ₁, which does not include the exponential growth, in order to analyze the real frequency and relative strength of the field components. Figure 8 presents wave spectra of perturbed magnetic field and the Poynting flux for k _y d = 0.25. For M _A = 0.45 in Figure 8A, only the waves at ω/ω _AI ∼ 1.2 have strong amplitude. This frequency is close to the eigenmode frequency of ω/ω _AI = 1.17 in Figure 2. The estimated growth rate near the V ₀ and V _A interfaces are identical with γ/ω _AI = 0.128. This growth rate is also in good agreement with the analytical results of γ/ω _AI = 0.142 in Section 3.

In order to examine the detailed wave properties, we plot spatial structures of the fluctuating magnetic field (b _1x , b _1y , and b _1z ) at ω/ω _AI = 1.2 in the middle row of Figure 8A. In this case, the compressional components (b _1x and b _1z ) maximize at the V _A and V ₀ interfaces and decay in the x-direction away from the interfaces. The b _1z and b _1x amplitudes at the two interfaces are comparable; thus, b _1z (x = d)/b _1z (x = 0) = 1.025. This ratio is almost identical to the amplitude ratio of the pressure A _p = 1.06 from Figure 2B.

On the other hand, the transverse component b _1y is enhanced at three different locations near V 0 k y x = − 0.022 and 0.033) and V _A interfaces (k _y = 0.22), where the wave frequency matches the Alfvén resonance condition (ω _AR ): ω = ω A R ± ≡ k y V 0 ± k ‖ V A .

Due to the finite width of the V ₀ interface near x = 0, ω A R − can be positive at the V ₀ interface; thus, two separate regions of enhanced wave power can occur corresponding to Doppler-shifted resonance with both Alfvén resonances. In this case, b _1y is significantly stronger than b _1z or b _1x , and b _1y /max(b _1z ) ∼ 28 at the V _A interface. Furthermore, strong field-aligned Poynting flux occurs at the interfaces, as shown in the lower panels of Figure 8A. The Poynting flux parallel (S _‖) and perpendicular (S _⊥) to B ₀ show that the wave energy predominantly flows along the magnetic field line at both interfaces. Since we launch compressional waves (with V _1x ) in the magnetosheath, and the growing KH waves are compressional waves, the amplitude enhancement of the magnetic transverse component and intense field-aligned Poynting flux at the interfaces are clear evidences of the mode conversion from the surface KH waves to the surface Alfvén waves.

For the higher M _A case in Figure 8B, the amplitude is maximized at ω = 2.57ω _AI , which is similar to the analytical value of ω = 2.6ω _AI in Section 3. The b _1z component maximizes near x = 0 and the secondary peak near k _y x = 0.25 becomes weaker. The b _1z amplitude ratio between the two interfaces is 0.64, which is in good agreement with the analytical value of A _p = 0.4. The b _1y component shows strong amplitude near k _y x = 0 and k _y x = 0.26. In this case, because the eigenmode frequency is higher than ω A R − , b _1y enhanced only at ω = ω A R + . The amplitude ratio between b _1y and b _1z is much lower than that for the case with M _A = 0.45 in Figure 8A, having b _1y /max(b _1z ) ∼ 6 at x ∼ d. Strong field-aligned flux S _‖ appears at the V ₀ interface in the bottom panels, but S _⊥ becomes stronger than the case of M _A = 0.45. The analytic eigenmode calculations predict that the PKHWs occur under (k _y d, M _A ) = (0.25, 0.75). The simulation results show that the mode conversion from the PKHWs to the surface Alfvén wave still occurs at each interface, but this process is less effective than that from the SKHWs.

For k _y d = 0.75, the waves have a strong amplitude peak near ω/ω _AI = 1.4 and 2.03 for M _A = 0.45 and 0.6, respectively, and the spatial structures of these waves are presented in Figure 9. In this case, the PKHWs and SKHWs are not separated anymore (See Figure 2) and we define the KH waves in the lower M _A as semi-SKHWs in Section 3. For M _A = 0.45 in Figure 9A, b _1z maximizes near x = 0 and a weak secondary peak appears near x = d. On the other hand, three amplitude peaks near k _y x = − 0.025 4, 0.018, and 0.754 appear in b _1y . The power ratio |b _1y /max(b _1z )| at the V ₀ interface is reduced to |b _1y /max(b _1z )| ∼ 7 from 23 for (k _y d, M _A ) = (0.25, 0.45) in Figure 8A. The enhancement of S _‖ is also seen at both interfaces and relatively strong S _⊥ also appears. Near x = 0 at the V ₀ interface, |S _⊥/S _∥| _x=0 is about 0.39, which is almost twice as large as |S _⊥/S _∥| _x=0 = 0.195 for (k _y d, M _A ) = (0.25, 0.45). Therefore, even though the compressional wave behavior of the semi-SKHWs is similar to the SKHWs, the mode conversion from semi-SKHWs becomes much weaker than that from the SKHWs.

For M _A = 0.6 in Figure 9B, b _1z decays along the x-direction from the V ₀ interface and no amplitude bump occurs at the V _A interface. The b _1y component shows a discontinuity at the V _A interface following the compressional Alfvén wave dispersion relation. Therefore, S _⊥ becomes comparable to S _‖ for k _y x > 0.75. In this case, the mode conversion at the V _A interface does not occur, but energy still flows along the magnetic field line at the V ₀ interface.

We also analyzed the cases for M _A for k _y d = 0, 0.25, and 0.75 at various M _A . Figure 10 shows the extracted eigenmode frequency and growth rate from the simulations. In this figure, the red and blue circled lines represent simulations, and the gray lines are taken from the eigenfrequency analysis from Figure 2 in Section 3. Although the boundary thicknesses used in Section 2 (slab) and Section 3 (width a) are different because the inhomogeneity scale length for the numerical simulation is much shorter than the wavelength, the eigenmode analytical and simulation results in ω and γ show excellent agreement.

We also calculate the amplitude ratios between the compressional component at the two interfaces (A _p ) and between transverse (b _1y ) and compressional ( b 1 x 2 + b 1 z 2 ) components at each interface. Due to the finite thickness of the boundary, two wave amplitude peaks can occur within the V ₀ interface as shown in Figures 8, 9, so we average the amplitude near x = 0, if ω = ω A R ± = k y V 0 ± k z V A . The simulated A _p and eigenfrequency calculations show good agreement with each other in Figure 11A,B. In particular, A _p of the SKHWs in Figure 11A are almost identical to the analytic calculations. Thus, these results confirm that the SKHWs occur with nearly the same amplitude at both interfaces, while the PKHWs only happen at the V ₀ interface. The amplitude ratio between the transverse and compressional components in Figure 11C,D suggests that the transverse magnetic component of the SKHWs is dominant. In other words, the mode conversion to the shear Alfvén wave from the SKHWs effectively occurs at the interfaces. The PKHWs in Figure 11C and semi-SKHWs and PKHWs in Figure 11D show that the transverse mode amplitudes are comparable to the compressional mode amplitude; therefore, weaker or no mode conversion occurs under given conditions at the V _A interface.