Nature and Scalings of Density Fluctuations of Compressible MHD Turbulence with Applications to the Solar Wind

The solar wind is a magnetized and turbulent plasma. Its turbulence is often dominated by Alfv\'enic fluctuations and often deemed as nearly incompressible far away from the Sun, as shown by in-situ measurements near 1AU. However, for solar wind closer to the Sun, the plasma $\beta$ decreases (often lower than unity) while the turbulent Mach number $M_t$ increases (can approach unity, e.g., transonic fluctuations). These conditions could produce significantly more compressible effects, characterized by enhanced density fluctuations, as seen by several space missions. In this paper, a series of 3D MHD simulations of turbulence are carried out to understand the properties of compressible turbulence, particularly the generation of density fluctuations. We find that, over a broad range of parameter space in plasma $\beta$, cross helicity and polytropic index, the turbulent density fluctuations scale linearly as a function of $M_t$, with the scaling coefficients showing weak dependence on parameters. Furthermore, through detailed spatio-temporal analysis, we show that the density fluctuations are dominated by low-frequency nonlinear structures, rather than compressible MHD eigen-waves. These results could be important for understanding how compressible turbulence contributes to solar wind heating near the Sun.


INTRODUCTION
The solar wind is a high-speed plasma flow emitted from the solar surface, as first predicted by Parker (1958). It was soon discovered that the solar wind is turbulent, containing fluctuations in a vast range of spatial and temporal scales. The free energy source for solar wind turbulence is typically at large scales, associated with activities near the Sun. (Energy can also be injected by instabilities at kinetic scales driven by beams or temperature anisotropies.) These large-scale perturbations will interact nonlinearly, causing cascade of energy to smaller scales and forming power-law spectra of magnetic, velocity and density fluctuations, which are believed to contribute to the solar wind heating.
In the classic picture of solar wind turbulence (see a recent review by Bruno & Carbone 2013), the fluctuations are treated as Alfvénic and incompressible, i.e., the divergence of velocity fluctuations is zero. Theories have been developed to deal with nearly incompressible turbulence in the solar wind and they have explained many of the key features of solar wind turbulence (Montgomery et al. 1987;Matthaeus & Brown 1988;Matthaeus et al. 1991;Zank & Matthaeus 1992;Zank et al. 2017). However, there exists a non-negligible component of compressible turbulence in the solar wind throughout the heliosphere, characterized by enhanced density fluctuations. Furthermore, there have been several theoretical and numerical studies on compressible MHD turbulence in astrophysical plasmas as well. Radio wave scintillation observations reveal a nearly Kolmogorov spectrum of density fluctuations in the ionized interstellar medium, suggestive of compressible MHD turbulence. Past theory and simulations have suggested that the density fluctuations are due to the slow mode and the entropy mode (Lithwick & Goldreich 2001). Many numerical simulations have been used to understand the nature of compressible MHD turbulence as well (Cho & Lazarian 2003;Kowal et al. MHD turbulence simulations, particularly with relatively high M t , to investigate how various turbulence quantities correlate among themselves.
The paper is organized as follows. Sec.2 describes the numerical model used in this study. Simulations results and analyses are presented in detail in Sec.3. We summarize main results in Sec.4, with discussion of unresolved issues and possible future directions.

NUMERICAL SIMULATION
To study compressible turbulence and density fluctuations in the context of the solar wind in the inner heliosphere (< 1 AU), we carry out a series 3D compressible MHD simulations. We select plasma parameters that represent typical solar wind conditions within 1 AU. We use the high-performance code ATHENA++ (Stone et al. 2020) to solve the ideal compressible MHD equations: with the polytropic equation of state pρ −γ = const, where γ is the polytropic index. Since we are interested in local properties of turbulence driven by large-scale free-energy input, our model does not include global structure of the solar wind. Instead, we assume a system with uniform background magnetic field and periodic boundary conditions. The simulation domain is an elongated box with L x = 8π, L y = L z = 2π and the background magnetic field is B 0 e x . Periodic boundary conditions are applied in all three directions. The Alfvén crossing time is For most of our runs we use 512 cells to resolve each dimension. Starting from a quiet uniform plasma (ρ 0 = 1, B 0 = 1 , and zero fluctuations), we perturb the MHD system by adding driving terms in both momentum and magnetic field equations (Eq. 2 and Eq. 3). We maintain the divergence free condition for the magnetic field by keeping ∇ · f b = 0 to machine precision (Gan et al. 2022). For all of the runs, we also keep ∇ · f v = 0 to model incompressible driving at large scales. f v and f b are applied at large scales where their wave number k ≤ k inj = 3. By varying the relative strength and phase of f v and f b , we are able to vary the Alfvén ratio r A ≡ |δv| 2 /|δB| 2 and cross helicity σ c ≡ (|z + | 2 − |z − | 2 )/(|z + | 2 + |z − | 2 ) (where z ± ≡ δv ± δB are Elsässer variables) of the injected fluctuations. To some degree, they replicate the observed solar wind conditions. By keeping the relative phase and strength of f v and f b and changing the magnitude of both quantities proportionally, we can also vary the driving strength and achieve different levels of turbulence. Typically, an adiabatic equation of state is adopted with γ = 5/3 to simulate MHD turbulence in the solar wind. But we will vary γ and study its effect on compressible turbulence and density fluctuations, as inspired by recent observations from PSP (Nicolaou et al. 2020). γ > 5/3 implies heating of the system by an external energy source (not included in the model) and γ < 5/3 implies cooling.
We should emphasize that, even though we have tried to keep the driving process ∇ · f v = 0, this does not mean that ∇ · (ρv) ≈ 0 during the evolution. For example, because we have used f b , the "added" Lorentz force J × B due to injection is typically non-zero. This leads to momentum variations that are not guaranteed to be ∇ · (ρv) = 0.
Although ∇ · f v = 0 is satisfied through the injection process, the density variations can be excited through the joint injection process of f v and f b . We will discuss this point in more details in the next section.

RESULTS
The key parameters and turbulent quantities of our MHD simulations are summarized in Table 1. Note that, in addition to these listed runs, we have also performed a few runs with higher resolution (1024 3 ) and many lower resolution (256 3 ) runs to map out the parameter space (see below), as well as to check the consistency of the results. In this section, we will first describe the evolution of turbulent quantities under different driving conditions. Second, we present the scaling of density fluctuations as a function of turbulent Mach number M t and examine how the scaling varies with key plasma/turbulence parameters. Third, we investigate the nature of the density fluctuations and their relation to compressible MHD waves and nonlinear structures. Table 1. Main input parameters (column 2-4) and turbulent quantities (column 5-9) for 3D MHD simulations. The system size is 8π × 2π × 2π with a uniform magnetic field B0 = 1 in the x direction. Turbulent quantities are measured in the late stage of each simulation when the system is in a quasi-steady state.

Evolution of Turbulent Quantities
Our nominal case is Run 1 with an initial β 0 = 1 and an adiabatic equation of state (i.e., γ = 5/3). The system is driven by Alfvénic incompressible perturbations at large scales (with k inj = 3) and it has a final cross helicity σ c = 0.93. Figure 1 (a) shows the time history of averaged velocity field, magnetic field and density fluctuations in Run 1, with where the brackets indicate the root-mean-square value over the whole simulation domain. These quantities are measures of the overall level of turbulence. With a continuous energy injection, all three fluctuation quantities increase gradually as a function of time. After an initial growth stage (0 < t < 60), the turbulence reaches a quasi-steady state with a saturated level of fluctuations. Between t = 70 and t = 90, we check the spectra of δv, δρ and δB and confirm that they remain nearly unchanged. Since we use an adiabatic equation of state, the total energy of the system is increasing with the continuous injection. After the turbulence is saturated, most of the input energy is converted into the internal energy of plasma, causing plasma heating in the later stage of the simulation (β increases from 1.0 at t = 0 to ∼ 1.03 at t = 90). Overall, the turbulence and plasma parameters are roughly controlled by the injection.
For Run 1, M t reaches 0.27 in the quasi-steady state and δρ reaches 0.13. By varying the strength of driving forces (f b and f v ), we can obtain turbulence with different levels of M t . In Run 2 (3), which uses the same parameters as Run 1 except stronger (weaker) driving forces, we obtain a larger (smaller) M t = 0.37 (M t = 0.19) and stronger (weaker) density fluctuations δρ = 0.19 (δρ = 0.09) (see Table 1 for details).
Note that although we keep the injected fluctuations Alfvénic with δv = δB, the plasma has freedom to respond to the injection and total δv is slightly larger than δB throughout the simulation. This can also be seen in Figure  1 (b), where r A first increases from 1.2 to 1.6 then reduces to 1.1 (r A = 1 for Alfvénic fluctuations) at the later stage. Similarly, we intend to have an imbalanced turbulence with σ c = 0.9, the resulting σ c varies from ∼ 0.8 at the beginning to ∼ 0.9 at the end of the simulation.
To contrast with Run 1, we carry out Run 4 with the same parameters except a lower β 0 = 0.2. This represents the solar wind condition much closer to the Sun. As shown in Figure 1 (c), the fluctuations saturate at t ∼ 40, with a similar level of velocity fluctuation (δv ∼ 0.24) to Run 1. It is notable that the density fluctuation generated in Run 4 is much higher than Run 1, with δρ ∼ 0.33 (vs. 0.13 in Run 1). This is the characteristic response of low-β plasma. Because we hold M A = δv/v A the same for both runs, Run 4, being lower β, achieves much higher turbulence March number M t with same δv. Another feature of Run 4, different from Run 1, is that the Alfvén ratio remains higher than unity (∼ 1.5 at the later stage, i.e. more energy goes into the velocity fluctuation than the magnetic fluctuation) throughout the simulation, even though the driving forces are Alfvénic. Similarly, the cross helicity remains around 0.7, lower than the value of Run 1. These quantities can be used to study the possible relationship between Alfvénicity and compressibility in the solar wind too (Bruno & Carbone 2013). As we emphasized earlier, even though the driving ensures that ∇ · f v = 0 and ∇ · f b = 0, the evolution of both magnetic field (thus Lorentz force) and momentum does not guarantee that ∇ · (ρv) = 0. Using Run 1 as an example, Figure 2 depicts the evolution of RMS values of various terms in the continuity equation (panel a) along with their probability density functions (PDFs) in panels b-d. (These quantities are calculated using the cell-based values.) It can be seen that, even though ∇ · v ≈ 0 when integrated over volume (i.e., the averaged value is near zero), there are large variances for ∇ · v and ∇ · (δρv), both of which contribute to large ∇ · (ρv) variations. It is also interesting to notice that the ∇ · v term has major contributions throughout the evolution (being dominant in the beginning), the ∇ · (δρ v) term also becomes important as time goes on, even becoming the bigger term after t > 70.

Scaling of Density Fluctuations
As shown above, the density variations depend strongly on the turbulent Mach number M t and show some dependence on plasma β and other parameters. To systematically study density fluctuations under various conditions relevant to the solar wind, we explore a large parameter space with different β, σ c , γ and M t . Before we present our main results, let us first examine the spectral properties of density fluctuations in the nominal case. Figure 3 shows a surface contour of δρ at t = 80 in Run 1. Clearly, δρ is anisotropic with much smoother variation along B 0 . Although the root-mean-square value of δρ is 0.13, there are localized structures with density enhancement or depletion with |ρ − ρ 0 | > 0.4. The power spectra of δρ as a function of k ⊥ or k are shown in Figure 4(a). First, we observe that most of the power is contained in the perpendicular direction. It is same for the power spectra of magnetic field and velocity fluctuations, as shown in Figure 4(a), which is due to the anisotropic cascade of MHD turbulence in magnetized plasmas (Shebalin et al. 1983;Goldreich & Sridhar 1995). Second the perpendicular spectra have two power law segments with a break point near k = 100, below which it follows a Kolmogorov law with k −5/3 , a signature of the inertial range. The roll-off of the spectra above the break point is due to numerical dissipation at small scales. Third, fluctuations at k ≤ 3 are mainly due to driving applied in the simulation. In the growing phase of turbulence (t < 60), the total divergence ∇ · (ρv) is dominated by the linear term ρ0∇ · v. In the saturation phase, the nonlinear term ∇ · (δρ v) is also quite strong.
In order to study how the density fluctuations will correlate with various quantities, we decide to apply a bandpass to the density power spectra, excluding fluctuations in both the injection range and the dissipation range. This is because the injection region dominates the spectral power, yet its dynamics is subject to driving forces. We can calculate density fluctuations in the inertial range as, Similarly, we can calculate magnetic field and velocity field fluctuations in the inertial range as δB and δv, respectively. Note that the injection wave number k inj is set by the driving forces (k inj = 3 in all our simulations), and the roll-over wave number k ro is set by the numerical scheme for solving the MHD equations and the resolution of the simulation. With other parameters being the same, Run 7 (1024 3 cells) has k ro ∼ 100, Run 1 (512 3 cells) has k ro ∼ 40, and Run 8 (256 3 cells) has k ro ∼ 25. But all three runs yield similar fluctuation levels in the inertial range, δB ∼ 0.13, δv ∼ 0.14 and δρ ∼ 0.08. They are understandably smaller than the total fluctuations which include the injection range (cf. Table 1).  . Power spectra of (a) total velocity field, magnetic field and density field fluctuations and (b) magnetic and velocity field components at t = 80 for Run 1 . In the inertial range (3 < k < 40), both the magnetic field and velocity field power spectra follow approximately the Kolmogorov law with k −5/3 ⊥ in the perpendicular direction.
Next, we study the relation of density fluctuation δρ with turbulent Mach number M t (≡ δv/c s ). Due to the limitation of computing resources, we use a large number of low-resolution (256 3 ) simulations to cover the parameter space with β = 0.2, 0.5, 1.0, 2.0, γ = 1.67, 2.70, σ c = 0.0, 0.3, 0.9 and a range of M t . The range of each parameter is inspired by observations. For example, PSP observations reported by Nicolaou et al. (2020) suggest that the solar wind protons follow an equation state with a wide range of γ, and the averaged γ increases from 1.7 at 1AU to 2.7 in the inner heliosphere. High σ c and low β are typical for the solar wind close to the Sun (e.g. Kasper et al. 2021), and σ c ∼ 0 and β ∼ 1 are typical for the solar wind at 1 AU. Figure 5 summarizes results from these low-resolution runs with different combinations of these key parameters. With other parameters (β, σ c and γ) fixed, we find δρ scales nearly linearly as a function of M t , i.e., The circles are data points extracted from simulations, and the dashed lines are linear fits to each set of runs. Note that although we have some control over σ c , it fluctuates in the simulations. So the numbers for σ c are not exact, but averaged values. β 0 is the initial value of β at t = 0, and β increases slightly in simulation due to heating by injection, as discussed above. Several features can be summarized from these results: First, the coefficient C has relatively weak dependence on β when using M t as shown by comparing different rows in Figure 5 (left panel). If we plot the scaling of density as a function of Alfvén Mach number (right panel), there will be strong β dependence, which is not surprising as M A ≡ M t γβ/2. Second, the dependence of C on both γ and σ c is relatively weak. More balanced turbulence (σ c = 0) yields slightly higher density fluctuations than the imbalanced turbulence (σ = 0.9). This trend is more prominent in runs with higher β values. As we approach the Sun, generally β decreases, γ increases and σ c increases. Changes of density fluctuation level can be complicated due to different effects of these parameters. Interestingly, for the same range of M t , density fluctuations near the Sun might be of similar magnitudes as those near 1 AU.

Nature of Density Fluctuations
Lastly, we examine the nature of density fluctuations in our simulations. Traditionally, density fluctuations are associated with compressible waves. Among the three linear MHD modes, Alfvén mode is incompressible, so only fast and slow modes can contribute to density fluctuations. Using linear mode decomposition in the spatial domain, fluctuations can be decomposed into the three MHD modes (Cho & Lazarian 2002). This decomposition is typically done on the velocity field, because the velocity fields from the three MHD modes form bases of a complete orthogonal system. We have performed such decomposition using frames near the end of Run 1. We find that the mode decomposition method gives about 75% the velocity fluctuation power in the Alfvén mode, about 20% in the slow mode and less than 5% in the fast mode. The spatial linear mode decomposition approach, however, does not take into account the temporal variations. For example, nonlinear structures, which do not belong to any linear mode, can still be decomposed spatially into three MHD modes. Therefore, as pointed out by Gan et al. (2022), this decomposition may overestimate the contribution from MHD waves.
To better analyze fluctuations in the simulations, we employ 4D FFT analysis (3D for space and 1D for time), in which both spatial and temporal information are used to calculate the Fourier power in the ω − k space (Gan et al. 2022). Using this approach, the power along the dispersion surfaces of MHD waves (numerically we sum over power within 10% of theoretical frequencies) will be regarded as "waves", and fluctuations elsewhere will be called non-wave structures. Figure 6 shows such an analysis of the density fluctuations for Run 1, where we use data from the whole simulation domain in a time window 50 < t < 80 (we exclude the initial growing phase). There is some power in two compressible modes: fast and slow modes. However, it is clear that the majority of the energy resides in the very low frequency structures, whose frequency is at or near zero and wave number is perpendicular to the background magnetic field (referred as "low-frequency structure" hereafter). This is more clear if we make a 2D cut in the ω − k ⊥ plane with k = 0 (panel b). Summing up the power in the low-frequency structure, we find that it comprises nearly 91% of the total power in the simulation box. MHD waves, on the other hand, contribute only to a small fraction, less than 9%. Furthermore, the majority of wave power resides in the slow mode (8%) and the power of fast mode is negligible (< 1%), which can also be seen in Panel (c). Similar results were reported by Gan et al. (2022). Finally, Panel (d) shows the anisotropy of the density fluctuation, consistent with Figure 4 (a).
Using 4D FFT analysis of the density fluctuations in the low-β Run 4, we obtain qualitatively similar results. As in Figure 7(a) and (c), the dispersion surface/curve of the slow mode (ω = k c s = k v A βγ/2) is lower with a lower β. Still, the nonlinear structures contribute to 97% of the total power, slow mode contributes 3%, and the fast mode is negligible.

CONCLUSION AND DISCUSSION
In this paper, using a large set of 3D MHD simulations, we investigate the development and properties of compressible turbulence driven by continuous large-scale perturbations. We have found: 1. Although the driving (at large scales) is maintained incompressible throughout the simulations, compressible fluctuations of varying magnitudes can still be generated via driving and they can be monitored via the RMS values of ρ, ∇·v and ∇·(δρv). In principle, compressible fluctuations can be generated via wave-wave interaction or mode conversion in an MHD system, and it is likely that some of these processes are occurring (particularly for low β simulations), but we believe that the density fluctuations in our simulations are primarily due to driving, which is enforced in both the momentum and magnetic fields. The turbulent Mach number in our simulations ranges from a few percent to the order of unity (including the driving range). After removing fluctuations in the driving range, the turbulent Mach number drops down to be less than 0.3.

2.
Simulations show a linear scaling of density fluctuation in the inertial range as a function of turbulence Mach number δρ = CM t , where the coefficient C depends weakly on σ c and γ. 3. Density fluctuations are dominantly caused by nonlinear structures with very low frequency and large k ⊥ . The nonlinear structures contribute to more than 90% of the fluctuation based on 4D FFT analysis whereas the actual compressible MHD waves contribute to a few percent in total power. These nonlinear structures could be mistakenly attributed to MHD modes when using the mode decomposition method based on the spatial variations alone. The contribution from the fast modes is negligible.
We have obtained a linear scaling of the density fluctuation as a function of turbulent Mach number based on our 3D MHD simulations. This is different from the prediction of nearly-incompressible MHD theory (Matthaeus & Brown 1988;Zank & Matthaeus 1992), i.e., δn ∝ M 2 t . The discrepancy may be caused by several reasons. First of all, NI-MHD theory assumes M t 1, whereas in our simulations M t ranges between 0.1 − 0.9 (including the injection range). Second, NI-MHD theory attributes the density fluctuations to pseudo-sound which is a zero-frequency structure with ∇ · v ≈ 0. As shown in our simulations, there are non-negligible ∇ · v fluctuations and the frequency of these nonlinear structures is low but not zero. The scaling of density fluctuation was also compared with NI-MHD theory using solar wind data (Matthaeus et al. 1991;Klein et al. 1993;Tu & Marsch 1994;Bavassano et al. 1995), but the spreads in both M t and δρ data points were too large to draw any firm conclusions. Note that, in-situ observation is 1D sampling through the solar wind turbulence, and the statistical properties are sensitive to sampling conditions such as the angle between the sampling line and the background magnetic field. They may also depend on plasma and turbulence parameters such as σ c and γ, as suggested by the simulations. Our linear scaling can be tested by observations if observational data can be grouped further by these parameters. Furthermore, how these results could vary depending on the driving processes needs to be further investigated (Gan et al. 2022).
Finally, we have identified that the majority of the density fluctuations comes from nonlinear structures that do not follow the dispersion relation of linear MHD waves, based on 4D Fourier analysis. However, the exact nature of the structures and their generation mechanism are not fully understood. How these structures are related to the pressure-balance-structures (Vasquez & Hollweg 1999), critically-balanced structures (Goldreich & Sridhar 1995) or psuedo-sound structures (Matthaeus et al. 1990) remains to be investigated. We leave these topics to future investigation.