High-resolution Simulations of the Inner Heliosphere in Search of the Kelvin–Helmholtz Waves

The Kelvin–Helmholtz instability (KHI) can be generated at velocity shears in plasmas. While shears are abundant in the solar wind, whether they generate KHI in situ remains an open question, because of the lack of models that can simultaneously resolve the global structure of the expanding solar wind and the local structure of much smaller-scale velocity shears. In this paper, we use the Grid Agnostic MHD for Extended Research Applications model whose high resolving power, in combination with a highly refined spatial grid, allowed us to extend the simulation from global scales roughly into the first decade of the inertial range (∼1.5 × 105 km, which we refer to as mesoscale). We employ this computational capability to extract from the simulation the local properties of radial and azimuthal solar wind velocity shears and investigate their KH stability using a linear dispersion relation, which includes both the finite width of the shear and plasma compressibility. We find that radial shears, which dominate the global structure of the inner heliosphere, are stabilized by compressibility. However, depending on the local Alfvén speed, sound speed, shear thickness, and the strength of the stabilizing azimuthal magnetic field, the azimuthal shears generated inside stream interaction regions could be KH-unstable. While our highly resolved simulation allowed us to analyze the local properties of the velocity shears, its resolution was still insufficient to confirm the instability. We argue that even higher resolution simulations are required to reproduce in situ generation of KHI at velocity shears in the solar wind.


Introduction
The inner heliosphere is a prominently turbulent plasma medium. It is generally believed that turbulence in this region of space plays a critical role in heating and accelerating the solar wind (Tu & Marsch 1995;Usmanov et al. 2011). Observations near Earth show that the solar wind turbulence spectrum has a wide range, extending over eight orders of magnitude (e.g., Kiyani et al. 2015). Turbulent fluctuations cascade from global heliospheric scales on the order of 1 astronomical unit (au) ≈ 150 × 10 6 km and larger down to smaller scales and begin dissipating at ion kinetic scales ( ( )  100 km). Over the years, different types of numerical simulations have been developed and have provided an essential tool for understanding the physics of solar wind plasma turbulence in the inner heliosphere. A large number of magnetohydrodynamic (MHD) simulations have been performed to investigate solar wind turbulence using simulation boxes with periodic boundary conditions (e.g., Matthaeus et al. 1991;Roberts et al. 1991;Shaikh & Zank 2010;Meyrand et al. 2016;Andrés et al. 2017). For example, Yang et al. (2019) studied the structure of waves and fluctuations in the solar wind using compressible MHD turbulence simulations in a threedimensional (3D) Cartesian box with periodic boundaries. MHD simulations using such idealized setups (i.e., periodic boundaries and Cartesian geometry) are capable of covering a large portion of the inertial range, and they can thus properly describe the energy spectrum and turbulence cascade under these simplified conditions. However, a significant complication in the study of solar wind turbulence arises from solar wind expansion, which is essential to turbulence evolution. To address this problem, an approach to modeling the spherical expansion was developed by Grappin et al. (1993) and Grappin & Velli (1996). These so-called expanding box model (EBM) simulations track a box co-moving with the radial mean flow. This approach has allowed studies of the nonlinear evolution of waves, turbulence, switchbacks, and stream structures in a more realistic representation of the inner heliospheric conditions (Grappin et al. 1993;Grappin & Velli 1996;Tenerani & Velli 2017;Shi et al. 2020).
The wealth of dedicated MHD simulations has significantly enhanced our knowledge about the fundamental physics of solar wind turbulence during the past decades. However, because of the large size of the inner heliosphere and the enormous range of scales over which turbulent fluctuations develop, all existing approaches inevitably make simplifying assumptions and by design are not intended to cover the entire inner heliosphere domain or to include realistic data-driven boundary conditions near the Sun. An independent track of solar wind model development has been ongoing in parallel, which has focused specifically on this last aspect, i.e., covering the global inner heliospheric domain and using data-driven boundary conditions (e.g., Riley et al. 2001;Merkin et al. 2016b;Wang et al. 2020). These models have largely been applied to studies of large-scale solar wind structures such as corotating interaction regions (CIRs) and coronal mass ejections (CMEs), and address both fundamental physical questions and space weather applications (e.g., Shiota et al. 2014;Feng et al. 2015;Wu et al. 2018;Pomoell & Poedts 2018;Werner et al. 2019;Iwai et al. 2021).
One can think of these two major types of solar wind models (i.e., dedicated turbulence models and global inner heliosphere models) as operating in different subranges of the solar wind turbulence spectrum. The global models remain largely in the injection range, and the smallest scales they resolve (in the radial direction) are on the order of the solar radius (7 × 10 5 km) or a relatively large fraction thereof, i.e., they hardly enter the inertial range (see below). At the same time, the dedicated models are designed to study the turbulent cascade and cover the inertial range, but they do not extend into the injection range where the realism of the global structure and data-driven boundary conditions become important. Given the growth of available computing power and the availability of MHD simulations with high resolving power, can we begin to bridge the gap between these two approaches? This is the primary question we address in this paper. More specifically, we consider whether MHD instabilities can be generated in situ in the solar wind given data-driven coronal boundary conditions, realistic global structures such as stream interactions, and a sufficiently high-resolution model, extending the simulation further into the inertial range than ever attempted before.
Velocity shears abound in the solar wind (e.g., Borovsky & Denton 2010). These velocity shears can lead to the generation of the Kelvin-Helmholtz (KH) instability (Coleman 1968;Helmholtz 1868;Kelvin 1871;Chandrasekhar 1961). Many observations show that KH waves frequently appear in various space plasma domains, including where the solar wind interacts with the Earthʼs magnetosphere (e.g., Fairfield et al. 2000;Hasegawa et al. 2004;Eriksson et al. 2016;Kieokaew & Foullon 2019), at magnetopauses of other planets like Mercury and Saturn (Slavin et al. 2008;Masters et al. 2010), on the flanks of CMEs (e.g., Foullon et al. 2011;Möstl et al. 2013;Johnson et al. 2014), and where magnetic flux interacts with solar prominences (e.g., Berger et al. 2010Berger et al. , 2017Ryutova et al. 2010;Hillier & Polito 2018). Recently, Solar Orbiter observed KH waves in the inner heliosphere at 0.69 au for the first time (Kieokaew et al. 2021). The KH instability was observed in a shear layer in the slow solar wind close to the heliospheric current sheet. In this paper, we focus on studying the stability of KH waves at velocity shears in the inner heliosphere.
There has been a wide range of analytical studies and numerical simulations investigating properties of the KH instability. Early investigations concentrated on the stability of linear KH waves for a tangential discontinuity with zero thickness in an incompressible flow (Chandrasekhar 1961). The effect of compressibility in the presence of magnetic fields was introduced in the case of a tangential discontinuity by Sen (1963Sen ( , 1964 and Fejer (1964). They showed that depending on the velocity shear speed relative to the magnetosonic speed, compressibility can have a stabilizing effect on KH waves. Compressibility introduces an upper limit for the range of velocity shears that can generate KH waves. Another limitation on the flow instability has been introduced by the magnetic field component parallel to the flow direction (Lau & Liu 1980;Miura & Pritchett 1982). The velocity shear should be large enough to overcome the stabilizing effect of the magnetic field line tension. Lerche (1966) emphasized the importance of considering the finite thickness of shear layers. Their paper showed that the MHD analysis of a zero-thickness tangential discontinuity leads to an inconsistency whereby the MHD approach fails for short wavelengths, thus indicating that the finite thickness of the velocity shear layer should not be neglected. Finally, Miura & Pritchett (1982) considered both the effects of a finite thickness of a shear layer and plasma compressibility on linear stability of KH waves.
In the present paper, we use a global 3D MHD model of the inner heliosphere, the Grid Agnostic MHD for Extended Research Applications (GAMERA; Zhang et al. 2019), to simulate the inner heliosphere starting from realistic datadriven coronal boundary conditions. The high resolving power of the model, in combination with a highly refined grid, allows us to extend our simulation from the global scale roughly into the first decade of the inertial range of solar wind turbulence (down to 1.5 × 10 5 km). Structures of this size are well in the range of mesoscales as defined by Viall et al. (2021). This capability allows us to start addressing the question of KH stability of velocity shears developed in the vicinity of stream interactions in the simulation using the theoretical framework of Miura & Pritchett (1982) as a guide. In Section 2, we describe the GAMERA model and the simulation setup used to model the inner heliosphere. Section 3 presents the results of our simulation, including analyzing the stability of the KH waves in the inner heliosphere using the Miura & Pritchett (1982) compressible linear theory. We compute the growth rate and the e-folding distance of the KH instability in our simulation and discuss the probability of generating KH waves in the solar wind. In Section 4, we summarize and conclude the paper.

Model
GAMERA is a reinvention of the Lyon-Fedder-Mobarry (LFM) MHD code (Lyon et al. 2004), which is well-known in the space physics community for its rather unique numerical algorithm featuring high-order reconstruction in arbitrary nonorthogonal curvilinear coordinates and a constrained transport scheme fulfilling the ∇ · B = 0 condition to machine precision. The LFM code was initially developed for global simulations of the terrestrial magnetosphere and has been used for decades for magnetospheric research (e.g., Merkin et al. 2019, and references therein). More recently, LFM was applied to global simulations of the inner heliosphere (Merkin et al. 2011;Pahud et al. 2012;Merkin et al. 2016aMerkin et al. , 2016b. GAMERA was written from scratch with the vision to retain the LFM MHD numerics, with some improvements, while adapting the code to modern and future supercomputer architectures. A detailed description of the underlying numerical methods was recently given by Zhang et al. (2019). A magnetospheric adaptation of GAMERA has recently been used to carry out simulations of ballooning-interchange instability in the Earth's magnetosphere, which required an unprecedentedly high resolution enabled by our numerical algorithms (Sorathia et al. 2020). In this paper we leverage the same MHD numerics to perform solar wind simulations at a similarly unprecedented resolution, using GAMERA adapted to the plasma and magnetic field conditions of the inner heliosphere.
Our simulation setup is similar to other models of the inner heliosphere, and it has been described elsewhere (Merkin et al. 2011(Merkin et al. , 2016b. As a summary, the inner boundary conditions in the MHD simulation, placed here at 21.5 R s (i.e., ∼0.1 au), are taken from the solution of the Wang−Sheeley−Arge (WSA) model of the corona as follows. The coronal magnetic field is initially obtained via extrapolation from photospheric magnetograph measurements by combining a potential field model and a heliospheric current sheet model. The Air Force Data Assimilate Photospheric Flux Transport (ADAPT; Arge et al. 2010) model is used to generate maps of the global photospheric magnetic field distribution. The distribution of the solar wind radial speed at 21.5 R s in the WSA model, which coincides with GAMERA's inner boundary, is then derived from a semi-empirical formula relating information about the magnetic field topology, such as the expansion of the magnetic flux tubes and proximity to coronal hole boundaries, to the speed of the solar wind outflow from the corona. The details of this procedure can be found in the work by Arge & Pizzo (2000) and Arge et al. (2004). The solar wind density is obtained using an empirical relationship between the density and velocity (Merkin et al. 2016b). To derive plasma temperature, we assume total pressure balance at GAMERA's inner boundary according to the equation where K B is the Boltzmann constant. B, n, and T are the magnetic field, solar wind number density, and temperature, respectively. The right-hand side of Equation (1) represents the vicinity of the heliospheric current sheet where we have neglected the magnetic field pressure. For n slow we picked the maximum number density on the grid that occurs at the position of the current sheet, yielding n slow = 1389 cm −3 . We then assumed T slow to be 3.6 × 10 5 K, which is reasonable for the slow wind. While GAMERA can utilize arbitrary hexahedral grids, in this work, we use a uniform spherical grid spanning 21.5 R s r 215 R s (here, R s is the solar radius and r is the heliocentric distance) in the radial dimension, 0.1 π θ 0.9 π in the polar dimension, and 0 f 2π in the azimuthal dimension. The grid used in this study had N r × N θ × N f = 1024 × 512 × 1024 cells (where N is the number of cells in each direction), leading to the cell size of Δr × Δθ × Δf = 0.2 R s × 0.3°× 0.35°. The simulation at this resolution is the main focus of this study, and we refer to this case as the "high-resolution case" hereafter. We note that due to the use of high-order reconstruction and an aggressive limiter in GAMERA, the above grid cell size is close to the physical scale size of solar wind structures the code is able to resolve, which is not the case for lower-order schemes (Zhang et al. 2019).
An inertial Sun-centered coordinate system is used such that the z-axis coincides with the solar rotation axis, the x-axis points from Sun to Earth at the beginning of the simulation, and the y-axis completes the right-handed coordinate system (see Merkin et al. 2011, for more details). In the inertial frame, the solar wind is purely radial and thus the off-radial components of the solar wind velocities are set to zero in the initial and boundary conditions of the simulation. The coronal boundary condition from the upper boundary of the WSA model is applied to GAMERA's inner boundary at 21.5 R s where the solar wind is already supersonic and super Alfvénic. The application of the inner boundary condition results in the solar wind blowing out the initial conditions specified on the grid. The simulation is then run long enough such that it relaxes to a solution that is in a steady state in the rotating frame . The high-resolution simulation was run on 128 nodes each consisting of 36 cores on the Cheyenne supercomputer at the US National Center for Atmospheric Research (NCAR). GAMERA uses hybrid parallelism, i.e., it uses Message Passing Interface (MPI) for grid decomposition and threading within each MPI rank using OpenMP. The simulation required about 11 hours of wall-clock time (i.e., ∼50,000 core hours).

Results
This section presents the results of the inner heliosphere simulations using GAMERA. Figure 1 shows latitude−longitude maps of the solar wind parameters and the radial component of the interplanetary magnetic field (IMF) in the second radial grid slice just above the inner boundary of the simulation (22 R s ). The radial velocity and magnetic field represent the coronal boundary values applied to the simulation from the WSA model, while the density and temperature were derived as described in Section 2. The map represents Carrington rotation 2210, centered on 2018 November 9, 12:32 UT, and corresponds to the first Parker Solar Probe (PSP) solar encounter. This time was chosen to compare the simulation results with the PSP observations that will be presented elsewhere. The heliospheric current sheet is defined as the B r = 0 isosurface, where B r is the radial magnetic field. Figure 2 shows the same variables as Figure 1 at 214 R s (∼1 au) and demonstrates how the solar wind and IMF features propagate and evolve from the inner boundary at 0.1 au (Figure 1) to the outer boundary of the simulation at 1 au. We observe that the simulation captures both large and mesoscale structures, e.g., folds in the heliospheric current sheet (Merkin et al. 2011(Merkin et al. , 2013, as well as field and plasma compression regions within stream interactions, developed as the solar wind propagates to 1 au. Figure 3 compares the high-resolution simulation described above with a lower-resolution simulation that has 1/4 as many cells in each dimension (256 × 128 × 256). The figure shows that the higher-resolution simulation not only captures the same features with increased sharpness but also introduces an emergent mesoscale structure that does not exist in the lowerresolution simulation. It should be noted that in both cases the same WSA solution with the angular resolution of 2.5°× 2.5°i s used whereby the WSA map is interpolated to the GAMERA grid to drive the corresponding GAMERA simulations. Since the resolution of the WSA map is even coarser than the lowerresolution GAMERA simulation, the mesoscale structure visible in the right panel of Figure 3 develops self-consistently in the MHD domain. However, in cases that we have analyzed so far, all structure seen in the MHD domain is ultimately seeded by the WSA boundary conditions. Larger-scale structures imposed by WSA, such as stream interactions, develop with the radial distance as the solar wind propagates away from the inner boundary, and mesoscale structure emerges from these interactions self-consistently, commensurate with the resolution of the MHD simulation. Figure 4 shows the radial and angular (nonradial) components of the velocity, as well as the radial component of the flow vorticity (i.e., where V θ and V f are polar (θ) and azimuthal (f) components of the velocity, respectively) at 1 au. These nonradial velocities and the corresponding flow vorticity are due to local deflections away from the radial direction within the stream interaction regions, generated at the interface between the fast and slow solar wind (for more details, see Pizzo 1978;Gosling & Pizzo 1999). In other words, the nonradial flows do not originate in the corona but develop in the MHD domain in response to the imposed coronal boundary condition that only exhibits a radial solar wind flow albeit with strong shears. Some of these shears appear in the vicinity of the heliospheric current sheet, while others develop at pseudostreamers. Figure 5 shows the different components of the solar wind velocity at different distances from the Sun. Note that the nonradial flows are negligible just above the inner boundary (22 R s ), and they grow as the solar wind propagates to larger heliocentric distances. The green lines in the last two panels show a stream interaction region that develops a predominantly azimuthal flow shear and will be the focus of the remainder of this paper. The relatively small flow away from the shear in the polar direction (middle column of Figure 5) allows us to consider the layer approximately as sheared only in the azimuthal direction and apply the linear KH stability analysis by Miura & Pritchett (1982), who considered an idealized onedimensional shear layer.

Analyzing the KH Waves
Compressibility (Miura & Pritchett 1982) and the magnetic field component parallel to the flow (Parker 1964) can stabilize the growth of KH waves. Study of the linear dispersion relation (Miura & Pritchett 1982) showed that in the parallel case (i.e., V 0 ∥ B; where V 0 and B are the sheared flow velocity and magnetic field vectors, respectively), both the Alfvén Mach number and sound Mach number are relevant to the stability of the MHD problem. The combined effects of the field line tension and compressibility require that the amplitude of the velocity shear be in the range 2V a < V 0 < 2V s for instability to grow, where V a and V s are Alfvén and sound speeds, respectively. Hence, plasma 2 ) is greater than one in this case. By contrast, only the fast magnetosonic mode is destabilized for the transverse case (k ∥ V 0 ⊥ B 0 ). As magnetosonic Mach 2 is the fast magnetosonic speed) increases, compressibility becomes important and the growth rate decreases. The shear layer becomes completely stable if M f > 2 (see more details in Figures 3 and 4 of the Miura & Pritchett 1982 paper).
Here, we use this analysis to estimate the possibility of KH wave generation in radial and azimuthal velocity shears in the solar wind. We note that the velocity shears in our simulation are neither purely parallel nor purely perpendicular to the magnetic field. Therefore, we use the dispersion relation analysis by Miura & Pritchett (1982) in a simplistic fashion, i.e., we consider whether a simulated shear would be unstable if it was a purely radial or a purely azimuthal shear. In the latter case, one can think of it as an azimuthal shear in the frame of the solar wind moving radially outward. This analysis is only used to roughly estimate the stability properties of the simulated shear from analytical theory; however, the ultimate answer on stability can only be provided by a numerical simulation, given that it is sufficiently resolved. We return to this issue below. Figure 6 represents an example of a velocity shear across a stream interaction region at 1 au. Figure 6(a) shows a velocity shear in the radial direction with an amplitude of ∼100 km s −1 . The approximate average values of characteristic velocities across the velocity shear are V a ∼ 22 km s −1 , V s ∼ 15 km s −1 , and V f ∼ 31 km s −1 (Figure 6(c)). The Alfvénic Mach number of this shear is M a ∼ 4.5 > 2, meaning that the radial shear is not stabilized by the magnetic tension. However, the sonic Mach number is M s ∼ 6.6 > 2, which leads to the stabilization of the radial shear by the compressibility. Note that this radial shear is relatively weak; for instance, Figure 5(a) clearly exhibits many stream interaction regions that are much more intense. Since the velocity jump across the shear is even larger in these regions, they can be expected to be stabilized by compressibility as well. We conclude that it is unlikely that linear KH waves driven by radial velocity shears could grow inside stream interaction regions. Figure 6(b) shows the change of the azimuthal velocity across the same shear with the amplitude of ∼25 km s −1 and the thickness essentially at the grid scale, i.e., Δ ∼ 1 R s (the cell size across the shear at 1 au is 1 au × dθ = 215 R s × 0.8π/512 ∼ 1 R s ). Based on the Miura & Pritchett (1982) dispersion relation analysis, such a shear, assuming for simplicity that it is purely perpendicular, would not be stabilized by compressibility since M f ∼ 25/22 < 2. The fastest growing mode would occur at kΔ ∼ 0.9 with its maximum growth rate equal to g = 0.12 max V 0 /Δ (see more details in Miura & Pritchett 1982, Figure 3). This leads to an estimated e-folding time and distance T = 1/γ ∼ 2.7 days and D ∼ 0.62 au (assuming a solar wind radial speed of 400 km s −1 ), respectively. Therefore, this azimuthal shear could become linearly unstable within 1 au. However, Figure 5 shows that the shear is not horizontally (azimuthally) aligned at smaller distances, but becomes progressively more so by the time it reaches 1 au. Therefore, if it becomes linearly unstable, it would likely happen beyond 1 au.
We mentioned above a number of caveats regarding the use of the dispersion analysis by Miura & Pritchett (1982). Among them are the following. (1) We picked a stream interaction region that aligns nearly azimuthally by 1 au, which minimizes the outflows in the direction perpendicular to the shear.
(2) We assumed the azimuthal shear to be purely perpendicular, i.e., we neglected the stabilizing role of the azimuthal component of the IMF. The parallel and transverse configurations only represent the limiting cases for the relative orientation of the velocity shear and magnetic field. Miura & Pritchett (1982) showed that the effect of the field line tension, compressibility, and the angle between flow speed and magnetic field (θ BV ) set the condition for the KH instability. The growth rate decreases as θ BV varies from 90°to 0°(see Miura & Pritchett 1982, Figure 7). The angle between the velocity shear and magnetic field at 1 au is about 45°, and the growth rate would be higher for the purely perpendicular case and lower for the purely parallel case.
Thus, as noted above, the estimates of the linear growth rate and e-folding distance we derived should only be used as a qualitative guide. In a complex and structured solar wind background such as simulated here, only numerical simulations will be able to address the question of shear layer stability. Our results above indicate that even if the azimuthal shear layer could be unstable in reality, it is still not resolved by our simulation, even given its high resolving power. This suggests that further simulations at even higher resolution are necessary. Indeed, because the growth rate is normalized to the shear width (Miura & Pritchett 1982), we expect that a thinner shear layer produced by a higher-resolution simulation would yield a proportionally higher growth rate. For example, a ten-fold increase of the resolution would lead to a drop of the e-folding distance to ∼0.06 au or only ∼13 R s . Therefore, there would be ample room for the KH waves to develop even within 1 au.
To confirm that our simulation remains stable we show in Figure 7 the radial velocity at 1 au at two different moments of time. Since the simulation is driven by an inner boundary condition that remains the same in the rotating frame, we expect the two maps to be identical unless there has been some dynamic evolution in the simulation domain as the solar wind propagated to 1 au. In Figure 7(b), we shift the solution back over time (i.e., f − ωΔt, where Δt = 50 and ω = 2π/25.38 days) to match Figure 7(a). The lack of noticeable differences between the two panels in Figure 7 shows that even our highresolution simulation remained stable. Therefore, even higher resolution simulations will be needed to confirm or refute the hypothesis that off-radial shears can become linearly unstable to KH waves in the solar wind.

Conclusions
Velocity flow shears forming inside stream interaction regions (including corotating interaction regions) can be considered a possible source of in situ generation of turbulence (Bavassano & Bruno 1989;Smith et al. 2011). For example, Bavassano & Bruno (1989) analyzed a velocity shear at 0.3 au and showed an increase in turbulence amplitudes associated with its passage. However, Borovsky & Denton (2010) analyzed different CIRs observed by the ACE spacecraft at 1 au and found no evidence of increased turbulence in the vicinity of velocity shears and thus concluded that the shears were not capable of driving turbulence. Therefore, the fundamental question of the extent to which turbulent fluctuations can be generated in situ in the solar wind, as opposed to propagating outward from the corona, still remains open. High-resolution global simulations of the inner heliosphere that are entering the inertial range of solar wind turbulence have now become feasible. They make it possible to investigate KH stability of velocity shears associated with  Figure 4 with the black rectangle. The green lines in the last two plots represent the stream interaction region and the predominantly azimuthal flow shear considered in the rest of the paper and analyzed further in Figure 6. large-scale stream interaction regions and thus start addressing the question above.
In this paper, we used the GAMERA MHD code to simulate the global inner heliosphere and analyzed the linear stability of KH waves in solar wind velocity shears. The simulation was performed at an exceedingly high resolution, resolving structure down to 0.2 R s (1.5 × 10 5 km) scale in the radial direction. The simulation results show the formation of mesoscale solar wind structures in stream interaction regions. The mesoscale structures are characterized by the off-radial flows that develop with increasing distance from the Sun and lead to the generation of radial flow vorticity. These structures are seeded by the larger-scale features of the WSA solution at the inner boundary of GAMERA and evolve self-consistently as the solar wind propagates away from the Sun. The structures are primarily located near the heliospheric current sheet and pseudostreamers, where the solar wind radial speed is highly variable, allowing interactions of plasma flows with high and low speeds.
Using the linear KHI dispersion relation by Miura & Pritchett (1982) as a guide, we first showed that, although, as expected, the most intense shears form at gradients of radial velocity, they are most likely stabilized by the solar wind plasma compressibility. We further demonstrated that shears in the azimuthal direction, which develop and intensify as the solar wind propagates away from the Sun, could become unstable to KH waves in the inner heliosphere, based on the shear properties at 1 au. However, the heliospheric magnetic field turns more azimuthal with the distance from the Sun as well, which would have a stabilizing effect on the KH waves generated in azimuthal shears. We hypothesize that the interplay of these effects ultimately determines whether KH waves can be generated in situ and in what range of radial distances it happens.
Despite its unprecedentedly high resolution, the simulation presented here still barely entered the inertial range of solar wind turbulence. We found that this resolution was still insufficient to resolve KH waves even if the azimuthal shear layers had been unstable according to linear theory. The simulation we presented used 0.5 billion cells, which is not particularly large (e.g., Ferrand et al. 2020). In addition, we used a simple spherical grid without leveraging GAMERA's capability to work on arbitrary hexahedral grids, where higher resolution can be adapted to the region of interest (see the terrestrial magnetosphere simulation by Sorathia et al. 2020 that leverages this capability). We conclude that even higher resolution simulations of the solar wind, capturing global to mesoscales and reaching deeper into the inertial range of solar wind turbulence, are now feasible. Pursuing these capabilities will be necessary to reproduce in situ generation of KHI in the solar wind. This, in turn, will help us better understand the physics of solar wind turbulence observed by PSP, Solar