The Evolution of Compressible Solar Wind Turbulence in the Inner Heliosphere: PSP, THEMIS, and MAVEN Observations

Turbulence is a unique nonlinear phenomenon in fluid and plasma flows that allows for the transfer of energy between different scales. Turbulence plays a major role in controlling the dynamical features of many astrophysical plasmas, such as accretion disks, star formation, solar wind heating, and energy transport in planetary magnetospheres (e.g., Balbus & Hawley 1998; Sorriso-Valvo et al. 1999; Schekochihin et al. 2009; Federrath & Klessen 2012; Hennebelle & Falgarone 2012; Tao et al. 2015; Kritsuk et al. 2017; Huang et al. 2020b, 2020c). Thanks to the availability of in situ measurements from various orbiting spacecraft, the solar wind provides a unique opportunity to investigate plasma turbulence (Osman et al. 2011; Howes et al. 2012; Alexandrova et al. 2013; Kiyani et al. 2015; Hadid et al. 2018; Andrés et al. 2019; Huang & Sahraoui 2019; Chen et al. 2020; Sahraoui et al. 2020; Huang et al. 2020b). A long-established challenge in the solar-wind community is the so-called heating problem. It is manifested by the fact that the solar wind's proton temperature decreases slowly as a function of the radial distance from the Sun, in comparison to the prediction of the adiabatic expansion model of the solar wind (Marsch et al. 1982; Vasquez et al. 2007; Pine et al. 2020). While several scenarios have been proposed to explain those observations (see, Marsch 1991; Matthaeus et al. 1999), the main candidate is certainly local heating of the solar wind plasma via turbulent cascade (Bruno & Carbone 2005; Matthaeus & Velli 2011). In this picture, the energy that is injected at the largest scales in the solar wind will cascade within the inertial range until it reaches the dissipation scales, where it is eventually converted into thermal heat of the plasma particles (see, Kiyani et al. 2015; Sahraoui et al. 2020). This framework has led to several investigations to estimate the energy cascade rate in the solar wind at different scales and different heliocentric distances using theoretical, numerical, and observational strategies.

The first exact theoretical relation (or law) for incompressible hydrodynamic (HD) turbulence was derived from the von Kármán–Howarth dynamical equation (von Kármán & Howarth 1938) and represents one of the very few exact results in turbulence theories (Frisch 1995). Under the assumption of homogeneity and full isotropy, the so-called 4/5 law (Kolmogorov 1941a, 1941b) predicts a linear scaling for the longitudinal, third-order structure function of the velocity field with the distance between points. This exact relation (valid only in the inertial range) gives an expression for the energy dissipation or cascade rate, ε, as a function of the structure functions of the turbulent fields (see, e.g., Monin & Yaglom 1975, and references therein). Galtier & Banerjee (2011) and Banerjee & Galtier (2014) generalized this exact result to compressible HD turbulence using isothermal and polytropic approximations, respectively. The authors found the presence of a new term that acts in the inertial range as a source (or a sink) for the mean energy cascade rate, in contrast with incompressible HD turbulence, where only flux terms act to transfer energy in the inertial range. Numerical results of supersonic isothermal HD turbulence have shown that these new source terms are smaller than the flux terms in the inertial range (Kritsuk et al. 2013).

The first generalization of these exact relations to a magnetized plasma was made by Politano & Pouquet (1998a, 1998b) using an incompressible magnetohydrodynamics (MHD) model. The validity of this exact result has been subjected to several numerical tests using direct numerical simulations (DNSs) of MHD turbulence (see, e.g., Boldyrev et al. 2009; Mininni & Pouquet 2009; Wan et al. 2010). Moreover, the exact law has been used to estimate the energy cascade rate (Sorriso-Valvo et al. 2007; Sahraoui 2008; Coburn et al. 2015) and the magnetic and kinetic Reynolds numbers (Weygand et al. 2007) in solar wind turbulence, and in large-scale modeling of the solar wind (Matthaeus et al. 1999; MacBride et al. 2008). Banerjee & Galtier (2013) derived an exact law for a two-point correlation function of the fields for isothermal compressible MHD turbulence, which was expressed in terms of flux or source terms. Recently, Andrés & Sahraoui (2017) revisited that work and provided a new derivation using classical plasma variables, i.e., the plasma density and velocity field and the compressible Alfvén speed. The new expression reported by Andrés & Sahraoui (2017) showed four types of terms that are involved in nonlinear cascade: hybrid and β-dependent terms, in addition to well-known flux and source terms (see, Andrés et al. 2019; Ferrand et al. 2020). It is this latter formulation that we shall use in the present study.

From an observational viewpoint, using a reduced form of the exact relation for compressible MHD turbulence and in situ measurements from the Time History of Events and Macroscale Interactions during Substorms (THEMIS) spacecraft (Auster et al. 2009), Banerjee et al. (2016) and Hadid et al. (2017) studied the role of compressibility in the energy cascade of solar wind turbulence. The authors found a more prominent role of density fluctuations in amplifying the energy cascade rate in slow, rather than fast, solar wind. Another interesting feature that has been evidenced in the Earth's magnetosheath is that density fluctuations reinforce the anisotropy of the energy cascade rate with respect to the local magnetic field (Hadid et al. 2018) and increase the cascade rate as it enters into sub-ion scales (Andrés et al. 2019). Recently, using Parker Solar Probe (PSP) observations during its first encounter with the Sun, Bandyopadhyay et al. (2020) computed the incompressible energy transfer rate between 55 and 35 solar radii. Their findings showed that the incompressible energy cascade rate obtained near the first perihelion was about 100 times higher than the average value at 1 au (e.g., Sorriso-Valvo et al. 2007; Hadid et al. 2017). Moreover, Andrés et al. (2020) computed the first estimation of the incompressible energy cascade rate at MHD scales in the plasma upstream of the Martian bow shock (at ∼1.38–1.67 au). Using Mars Atmosphere and Volatile EvolutioN (MAVEN) observations, the authors found that the nonlinear cascade of energy at MHD scales is slightly amplified when proton cyclotron waves are present in the plasma. However, in all these recent cases, only the incompressible cascade rate was estimated.

The main goal of the present paper is to generalize and extend previous observational studies using a more complete theory of turbulence to investigate the compressible energy cascade rate, ε_C, at MHD scales at different heliocentric distances. In particular, we use three data sets: observations from PSP at ∼0.2–0.4 au, THEMIS around ∼1 au, and MAVEN at ∼1.5–1.7 au. The paper is organized as follows: in Section 2 we describe the compressible MHD set of equations and recall briefly the main steps to derive the exact laws for fully developed compressible MHD turbulence (and its incompressible limit). In Section 3, we present the three observational data sets and the selection criteria used in the present work. Finally, in Section 4 we discuss our main observational results and their physical implications on solar wind turbulence.

2.1. Compressible MHD Equations

The three-dimensional (3D) compressible MHD equations are the continuity equation for the mass density, ρ, the momentum equation for the velocity field, u , in which the Lorentz force is included, the induction equation for the magnetic field, B , and the differential Gauss' law. These equations can be written as (see, e.g., Marsch & Mangeney 1987; Andrés et al. 2017),

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho }{\partial t}=-{\boldsymbol{\nabla }}\cdot (\rho {\boldsymbol{u}}),\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\partial {\boldsymbol{u}}}{\partial t} & = & -{\boldsymbol{u}}\cdot {\boldsymbol{\nabla }}{\boldsymbol{u}}+{{\boldsymbol{u}}}_{{\rm{A}}}\cdot {\boldsymbol{\nabla }}{{\boldsymbol{u}}}_{{\rm{A}}}-\displaystyle \frac{1}{\rho }{\boldsymbol{\nabla }}(P+{P}_{M})\\ & & -{{\boldsymbol{u}}}_{{\rm{A}}}({\boldsymbol{\nabla }}\cdot {{\boldsymbol{u}}}_{{\rm{A}}})+{{\boldsymbol{f}}}_{k}+{{\boldsymbol{d}}}_{k},\end{array}\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {{\boldsymbol{u}}}_{{\rm{A}}}}{\partial t}=-{\boldsymbol{u}}\cdot {\boldsymbol{\nabla }}{{\boldsymbol{u}}}_{{\rm{A}}}+{{\boldsymbol{u}}}_{{\rm{A}}}\cdot {\boldsymbol{\nabla }}{\boldsymbol{u}}-\displaystyle \frac{{{\boldsymbol{u}}}_{{\rm{A}}}}{2}({\boldsymbol{\nabla }}\cdot {\boldsymbol{u}})+{{\boldsymbol{f}}}_{m}+{{\boldsymbol{d}}}_{m},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{{\boldsymbol{u}}}_{{\rm{A}}}\cdot {\boldsymbol{\nabla }}\rho +2\rho ({\boldsymbol{\nabla }}\cdot {{\boldsymbol{u}}}_{{\rm{A}}})=0,\end{eqnarray} \tag{ 4 }$$

where we have defined the compressible Alfvén velocity as ${{\boldsymbol{u}}}_{{\rm{A}}}\equiv {\boldsymbol{B}}/\sqrt{4\pi \rho }$. In this manner, both field variables, u and u _A, are expressed in speed units. For the sake of simplicity, we assume that the plasma follows an isothermal equation of state, i.e., $P={c}_{s}^{2}\rho$, where c_s is the sound speed (temperature-dependent), which allows us to close the hierarchy of the fluid equations (no energy equation is further needed), and ${P}_{M}\,\equiv \rho {u}_{{\rm{A}}}^{2}/2$ is the magnetic pressure. Finally, f _k,m are, respectively, a mechanical forcing and the curl of the electromotive large-scale forcing, and d _k,m are, respectively, the small-scale kinetic and magnetic dissipation terms.

2.2. Exact Relation in Compressible MHD Turbulence

Using Equations (1–4), and following the usual assumptions for fully developed homogeneous turbulence (i.e., infinite kinetic and magnetic Reynolds numbers and a steady state with a balance between forcing and dissipation; Galtier & Banerjee 2011; Banerjee & Kritsuk 2018), an exact relation for compressible MHD turbulence can be obtained as,

$$\begin{eqnarray}&&-2{\varepsilon }_{C}=\displaystyle \frac{1}{2}{{\boldsymbol{\nabla }}}_{{\ell }}\cdot {{\boldsymbol{F}}}_{{\rm{C}}}+{{\rm{S}}}_{{\rm{C}}}+{{\rm{S}}}_{{\rm{H}}}+{{\rm{M}}}_{\beta },\end{eqnarray} \tag{ 5 }$$

where ε_C is the total compressible energy cascade rate and F _C, S_C, S_H, and M_β represent the total compressible flux, source, hybrid, and β-dependent terms, respectively (for a detailed derivation see, Andrés & Sahraoui 2017). These terms are defined as,

$$\begin{eqnarray}\begin{array}{rcl}{{\boldsymbol{F}}}_{{\rm{C}}} & \equiv & \,\langle [(\delta (\rho {\boldsymbol{u}})\cdot \delta {\boldsymbol{u}}+\delta (\rho {{\boldsymbol{u}}}_{{\rm{A}}})\cdot \delta {{\boldsymbol{u}}}_{{\rm{A}}}]\delta {\boldsymbol{u}}\\ & & -[\delta (\rho {\boldsymbol{u}})\cdot \delta {{\boldsymbol{u}}}_{{\rm{A}}}+\delta {\boldsymbol{u}}\cdot \delta (\rho {{\boldsymbol{u}}}_{{\rm{A}}})]\delta {{\boldsymbol{u}}}_{{\rm{A}}}\rangle \\ & & +2\langle \delta e\delta \rho \delta {\boldsymbol{u}}\rangle ,\end{array}\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}\begin{array}{rcl}{{\rm{S}}}_{{\rm{C}}} & \equiv & \,\langle [{R}_{E}^{{\prime} }-\displaystyle \frac{1}{2}({R}_{B}^{{\prime} }+{R}_{B})]({\boldsymbol{\nabla }}\cdot {\boldsymbol{u}})\\ & & +[{R}_{E}-\displaystyle \frac{1}{2}({R}_{B}+{R}_{B}^{{\prime} })]({\boldsymbol{\nabla }}^{\prime} \cdot {\boldsymbol{u}}^{\prime} )\rangle \\ & & +\langle [({R}_{H}-{R}_{H}^{{\prime} })-\bar{\rho }({\boldsymbol{u}}^{\prime} \cdot {{\boldsymbol{u}}}_{{\rm{A}}})]({\boldsymbol{\nabla }}\cdot {{\boldsymbol{u}}}_{{\rm{A}}})\\ & & +[({R}_{H}^{{\prime} }-{R}_{H})-\bar{\rho }({\boldsymbol{u}}\cdot {{\boldsymbol{u}}}_{{\rm{A}}}^{{\prime} })]({\boldsymbol{\nabla }}^{\prime} \cdot {{\boldsymbol{u}}}_{{\rm{A}}}^{{\prime} })\rangle ,\end{array}\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}\begin{array}{rcl}{{\rm{S}}}_{{\rm{H}}} & \equiv & \,\Space{0ex}{0.25ex}{0ex}\langle \left(\displaystyle \frac{{P}_{M}^{{\prime} }-P^{\prime} }{2}-E^{\prime} \right)({\boldsymbol{\nabla }}\cdot {\boldsymbol{u}})\\ & & +\left(\displaystyle \frac{{P}_{M}-P}{2}-E\right)({\boldsymbol{\nabla }}^{\prime} \cdot {\boldsymbol{u}}^{\prime} )\Space{0ex}{0.25ex}{0ex}\rangle \\ & & +\langle H^{\prime} ({\boldsymbol{\nabla }}\cdot {{\boldsymbol{u}}}_{{\rm{A}}})+H({\boldsymbol{\nabla }}^{\prime} \cdot {{\boldsymbol{u}}}_{{\rm{A}}}^{{\prime} })\rangle \\ & & +\displaystyle \frac{1}{2}\Space{0ex}{0.25ex}{0ex}\langle \left(e^{\prime} +{\displaystyle \frac{{u}_{{\rm{A}}}}{2}}^{{\prime} 2}\right)\left[{\boldsymbol{\nabla }}\cdot (\rho {\boldsymbol{u}})\right]\\ & & +\left(e+{\displaystyle \frac{{u}_{{\rm{A}}}}{2}}^{2}\right)\left[{\boldsymbol{\nabla }}^{\prime} \cdot (\rho ^{\prime} {\boldsymbol{u}}^{\prime} )\right]\Space{0ex}{0.25ex}{0ex}\rangle ,\end{array}\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{{\rm{M}}}_{\beta }\equiv -\displaystyle \frac{1}{2}\langle {\beta }^{-1^{\prime} }{\boldsymbol{\nabla }}^{\prime} \cdot (e^{\prime} \rho {\boldsymbol{u}})+{\beta }^{-1}{\boldsymbol{\nabla }}\cdot (e\rho ^{\prime} {\boldsymbol{u}}^{\prime} )\rangle ,\end{eqnarray} \tag{ 9 }$$

where we have defined the total energy (i.e., the free energy) and the density-weighted cross-helicity per unit volume, respectively, as,

$$\begin{eqnarray}&&E({\boldsymbol{x}})\equiv \displaystyle \frac{\rho }{2}({\boldsymbol{u}}\cdot {\boldsymbol{u}}+{{\boldsymbol{u}}}_{{\rm{A}}}\cdot {{\boldsymbol{u}}}_{{\rm{A}}})+\rho e,\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&H({\boldsymbol{x}})\equiv \rho ({\boldsymbol{u}}\cdot {{\boldsymbol{u}}}_{{\rm{A}}}),\end{eqnarray} \tag{ 11 }$$

and their associated two-point correlation functions as,

$$\begin{eqnarray}&&{R}_{E}({\boldsymbol{x}},{\boldsymbol{x}}^{\prime} )\equiv \displaystyle \frac{\rho }{2}({\boldsymbol{u}}\cdot {\boldsymbol{u}}^{\prime} +{{\boldsymbol{u}}}_{{\rm{A}}}\cdot {{\boldsymbol{u}}}_{{\rm{A}}}^{{\prime} })+\rho e^{\prime} ,\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{R}_{H}({\boldsymbol{x}},{\boldsymbol{x}}^{\prime} )\equiv \displaystyle \frac{\rho }{2}({\boldsymbol{u}}\cdot {{\boldsymbol{u}}}_{{\rm{A}}}^{{\prime} }+{{\boldsymbol{u}}}_{{\rm{A}}}\cdot {\boldsymbol{u}}^{\prime} ),\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{R}_{B}({\boldsymbol{x}},{\boldsymbol{x}}^{\prime} )\equiv \displaystyle \frac{\rho }{2}({{\boldsymbol{u}}}_{{\rm{A}}}\cdot {{\boldsymbol{u}}}_{{\rm{A}}}^{{\prime} }).\end{eqnarray} \tag{ 14 }$$

We have also defined the internal compressible energy for an isothermal plasma, i.e., $e\equiv {c}_{s}^{2}\mathrm{ln}(\rho /{\rho }_{0})$, where ρ₀ is a constant (reference) mass density. In all cases, a prime denotes a field evaluation at ${\boldsymbol{x}}^{\prime} ={\boldsymbol{x}}+{\ell }$ (ℓ being the displacement vector) and an angular bracket 〈 · 〉 denotes an ensemble average. It is worth mentioning that the properties of spatial homogeneity imply (assuming ergodicity) that the results of averaging over a large number of realizations can be obtained equally well by averaging over a large region of space for one realization (Batchelor 1953). In analyses of spacecraft data, this generally amounts to time averaging based on the Taylor hypothesis (e.g., Taylor 1938a, 1938b; Huang & Sahraoui 2019; Treumann et al. 2019). We have introduced the usual increments and local mean definitions, i.e., $\delta \alpha \equiv \alpha ^{\prime} -\alpha$ and $\bar{\alpha }\equiv (\alpha ^{\prime} +\alpha )/2$ (where α is any scalar or vector function), respectively. Finally, we recall that the derivation of the exact law (Equation (5)) does not require the assumption of isotropy and it is independent of the dissipation mechanisms acting in the plasma (assuming that the dissipation acts only at the smallest scales in the system) (see also, Galtier & Banerjee 2011; Andrés et al. 2016a, 2016b).

The quantity in Equation (6) is associated with the energy flux, and is the usual term present in the exact law of incompressible turbulence (Politano & Pouquet 1998a). This term is written as a global divergence of products of increments of different variables. It is worth mentioning that the total compressible flux (Equation (6)) is a combination of two terms of different nature, a Yaglom-like term,

$$\begin{eqnarray}\begin{array}{rcl}{{\boldsymbol{F}}}_{1{\rm{C}}} & \equiv & \,\langle [(\delta (\rho {\boldsymbol{u}})\cdot \delta {\boldsymbol{u}}+\delta (\rho {{\boldsymbol{u}}}_{{\rm{A}}})\cdot \delta {{\boldsymbol{u}}}_{{\rm{A}}}]\delta {\boldsymbol{u}}\\ & & -[\delta (\rho {\boldsymbol{u}})\cdot \delta {{\boldsymbol{u}}}_{{\rm{A}}}+\delta {\boldsymbol{u}}\cdot \delta (\rho {{\boldsymbol{u}}}_{{\rm{A}}})]\delta {{\boldsymbol{u}}}_{{\rm{A}}}\rangle ,\end{array}\end{eqnarray} \tag{ 15 }$$

which is the compressible generalization of the incompressible term (see, Politano & Pouquet 1998a, 1998b), and a new purely compressible flux term,

$$\begin{eqnarray}&&{{\boldsymbol{F}}}_{2{\rm{C}}}\equiv \,2\langle \delta \rho \delta e\delta {\boldsymbol{u}}\rangle ,\end{eqnarray} \tag{ 16 }$$

which is a new contribution to the energy cascade rate due to the presence of density fluctuations in the plasma (see, Banerjee & Galtier 2013; Andrés & Sahraoui 2017).

The purely compressible source terms in Equation (7), i.e., those proportional to the divergence of the Alfvén and kinetic velocity fields (and which involve the two-point correlation functions R_E , R_B , and R_H ), may act as a source (or a sink) for the mean energy cascade rate in the inertial range. The hybrid term offers the freedom to be written either as a flux- or as a source-like term. However, when written as a flux-like term it cannot be expressed as the product of increments, as the classical flux in incompressible HD and MHD turbulence (von Kármán & Howarth 1938; Kolmogorov 1941a, 1941b; Chandrasekhar 1951; Politano & Pouquet 1998a, 1998b) or their counterparts in Equation (6). The mixed β-dependent term (transformed into a flux-like term by Banerjee & Galtier (2013) under certain conditions) has no counterpart in compressible HD turbulence (Galtier & Banerjee 2011; Banerjee & Galtier 2014) and cannot, in general, be expressed as a purely flux or source term. It is worth recalling that these three type of terms can be estimated only using multi-spacecraft data and techniques, since they include local vector divergences (Andrés et al. 2019). In particular, because of their complex form, we could not find any scaling argument that would prove that the source terms are smaller than the flux terms. In our data, since they were measured by a single spacecraft, we could not estimate the non-flux terms. However, using Magnetospheric Multiscale Mission (MMS) observational data of the Earth's magnetosheath it was shown that they were indeed much smaller than the flux terms (Andrés et al. 2019). This observational result agrees with findings from HD and MHD simulations (Kritsuk et al.2013; Andrés et al. 2018; Ferrand et al. 2020). In particular, for high resolution (10,048³ grid points) simulations of HD turbulence (Ferrand et al. 2020) it was shown that the source terms can be comparable to the flux terms only in a supersonic regime of turbulence and within some localized regions of space characterized by strong density filaments. Therefore, in the present paper, to estimate the compressible energy cascade rate (Equation (5)), we shall consider only the flux terms in Equations (15) and (16).

Assuming statistical isotropy, we can integrate Equation (5) over a sphere of radius ℓ to obtain a scalar relation for isotropic turbulence. In compact form, Equation (5) can be cast as,

$$\begin{eqnarray}&&-\displaystyle \frac{4}{3}{\varepsilon }_{C}{\ell }={{\rm{F}}}_{1{\rm{C}}}+{{\rm{F}}}_{2{\rm{C}}},\end{eqnarray} \tag{ 17 }$$

where ${{\rm{F}}}_{1{\rm{C}}}+{{\rm{F}}}_{2{\rm{C}}}\equiv ({{\boldsymbol{F}}}_{1{\rm{C}}}+{{\boldsymbol{F}}}_{2{\rm{C}}})\cdot {\hat{{\boldsymbol{V}}}}_{\mathrm{sw}}$ is the flux term projected into the mean plasma flow velocity field, V _sw. The incompressible limit is easily recovered for ρ → ρ₀,

$$\begin{eqnarray}&&-\displaystyle \frac{4}{3}{\varepsilon }_{I}{\ell }={{\rm{F}}}_{I},\end{eqnarray} \tag{ 18 }$$

where F_I is the projection of F _I = ρ₀〈[(δ u )² + (δ B )²]δ u − 2(δ u · δ B )δ B 〉 along the mean plasma flow velocity field. Equation (18) corresponds to the exact relation for fully developed incompressible MHD turbulence (see, Politano & Pouquet 1998a, 1998b). Here B is expressed in velocity units and ε_I is the incompressible energy cascade rate. One can note that F _I depends only on the increments of the magnetic (and velocity) field, although the total magnetic field has been considered in the derivation. Finally, assuming the Taylor hypothesis (i.e., V τ ≡ ℓ, where V is the mean plasma flow speed), Equations (17) and (18) can be expressed as a function of the time lag, τ.

In order to analyze the solar wind turbulence at different heliocentric distances, our data intervals for PSP and THEMIS/MAVEN were divided into a series of samples of equal duration of 30 and 35 minutes, respectively. This particular time duration ensures having at least one correlation time of the turbulent fluctuations for each particular heliocentric distance (see, Hadid et al. 2017; Marquette et al. 2018; Parashar et al. 2020). Moreover, we avoid intervals that contained significant disturbances or large-scale gradients (e.g., coronal mass ejections or interplanetary shocks). We further considered only intervals that did not show large fluctuations of the energy cascade rate over the MHD scales; typically, we retained events with std(ε_I )/mean(∣ε_I ∣) < 0.75.

Table 1 shows a brief description of the data used. Our analysis of the PSP observations (Bale et al. 2016; Fox et al. 2016; Kasper et al. 2016; Bale et al. 2019; Kasper et al. 2019; Case et al. 2020) involved two particular data sets: one covering the period between 2018 November 1 and November 10, which was dominated essentially by slow-wind flows, and the other one between 2018 November 15 and November 21 dominated by high-speed flows. In both sets, spurious data (i.e., high artificial peaks) in the Solar Probe Cup (SPC) moments (see, Kasper et al. 2016) were removed using a linear interpolation (see, Bandyopadhyay et al. 2020; Parashar et al. 2020) and the data set was re-sampled to a 0.873 s time resolution. For the data at 1 au, we used the same data set analyzed by Hadid et al. (2017), i.e., the period 2008–2011, where a large survey of the THEMIS data (Auster et al. 2008; McFadden et al. 2008, 2009) has been reported. This data set was re-sampled to a 3 s time resolution and we only used events that fulfilled our criteria, i.e., covering both fast and slow solar wind. Finally, magnetic field observations from the MAVEN spacecraft (Coburn et al. 2015; Jakosky et al. 2015a, 2015b) were analyzed to discriminate events in the pristine solar wind with and without wave activity (Russell et al. 1990; Mazelle et al. 2004). As discussed in the 1, recently, the incompressible energy cascade rate has been investigated in a Martian plasma environment (Andrés et al. 2020). Here, we used the same data set for the solar-wind observations that did not contain the presence of proton cyclotron waves (also see, Romanelli et al. 2013, 2016; Halekas et al. 2020; Romeo et al. 2021).

As mentioned by Huang et al. (2020a), SPC on board PSP reports proton measurements derived from both moment and nonlinear fitting algorithms. On the one hand, the moment algorithm returns a single, isotropic proton population. On the other hand, the nonlinear fitting algorithm returns a proton core and a proton beam population, where, usually the proton core corresponds to the peak of the solar wind proton velocity distribution function, and the beam corresponds to its "shoulder". Kasper (2002) suggested that the nonlinear fitting technique provides far more information than the moment algorithm, but the moment algorithm is used for its simplicity. In the present work, we use the measurements provided by SPC from the moment fitting algorithm. Also, it is worth mentioning that the local solar wind temperature from the Solar Wind Ion Analyzer (SWIA) instrument (on board MAVEN) and the electrostatic analyzer (ESA) instrument (on board THEMIS) may overestimate the temperature moments by a factor of ∼2 due to the trace presence of alpha particles (Halekas et al. 2017). In fact, the local temperature enters in the compressible exact law only in the compressible term in Equation (17) through the definition of the sound speed in the internal energy, i.e., $e={c}_{s}^{2}\mathrm{log}(\rho /{\rho }_{0})$, where c_s is the isothermal sound speed. To investigate the impact of a possible overestimation of the solar wind temperature in the MAVEN observations, we have computed the compressible cascade component by taking into account half of the measured temperatures. While in some cases the compressible component decreases (by a factor ∼2), the statistical results, trends, and conclusions in this paper are not affected by this potential temperature overestimation.

Figure 1 shows the occurrence rates for all the analyzed events in the three data sets of the number density, and the absolute values of the velocity and Alfvén velocity field. As we expect, while the THEMIS and MAVEN data sets have similar ranges of mean values for the number density, proton, and Alfvén velocity absolute values, the probability density functions (PDFs) of the PSP observations show clear increases in the number density (two orders of magnitude) and the Alfvén velocity.

Zoom In Zoom Out Reset image size

4.1. The Compressible Cascade Rate and MHD Scales

In order to compute the compressible cascade (Equation (17)), we constructed temporal correlation functions of the different turbulent fields at different time lags in the interval [4, 2100] s. Figure 2 shows the compressible flux term in Equation (17) (absolute value) as a function of time lag for three cases chosen randomly from our list. As a reference, we show a linear scaling in dashed lines and the timescales that we considered to compute the mean values in the MHD inertial range in gray color bands. Figure 3 shows the non-signed compressible energy cascade rate, ∣ε_C ∣, as a function of the time lag, τ, for all the events in each mission. As we expect, our findings show a clear amplification of the compressible cascade as we approach the Sun, which can be explained by an increase of the magnetic and velocity field fluctuations closer to the Sun (see Figure 4). In order to quantify statistically the increment in the cascade rates, we consider the mean of the absolute value at the largest MHD scales (τ ∈ [100, 1800] s) as a representative value of each event.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

4.2. Compressible and Incompressible Energy Cascade Rates

Figure 4 shows 〈∣ε_C∣〉 as a function of 〈∣ε_I∣〉 for the PSP (triangles), THEMIS (squares), and MAVEN (circles) observations. The color bar represents the compressibility (in percent) of each event, defined as 〈∣n_i − n₀∣〉/n₀ (where n₀ is the mean number density). The insets correspond to the PDFs of the absolute values of the velocity and magnetic field fluctuations. In particular, the level of density, velocity, and magnetic fluctuations increase as we approach the Sun (e.g., Bruno & Carbone 2005; Matthaeus & Velli 2011). On the other hand, considering that the plasma compressibility increases up to 25% at PSP's perihelion, we observe moderate increases of the total compressible cascades with respect to the incompressible cascades. This is in agreement with previous work regarding the solar wind and the Earth's magnetosheath (Hadid et al. 2017; Andrés et al. 2019; Huang & Sahraoui 2019).

Figures 5 and 6 show the compressible component 〈∣ε_2C∣〉 (related to Equation (16)) as a function of the Yaglom generalization component 〈∣ε_1C∣〉 (related to Equation (15)), normalized to the total incompressible and compressible cascade rates, respectively. When the density fluctuations are small (≤5%), the compressible Yaglom-like component, ε_1C, coincides with the incompressible component, ε_I , while the compressible term, ε_2C, is negligible. However, when the density fluctuations increase in the plasma (up to 25%), the dominant component is given by a competition between ε_1C and ε_2C. We recall that, independent of the level of the density fluctuations, we assumed that the source, hybrid, and β-dependent terms are negligible in the MHD inertial range based on the simulation results of Andrés et al. (2018), and thus are not estimated here.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

4.3. Relation Between the Cascade Rate and Local Temperature

In order to connect the turbulent cascade rate and the local ion temperature (e.g., Livadiotis et al. 2020; Wu et al. 2020), we studied the correlation of the compressible cascade rate with the ion temperatures for slow and fast solar-wind plasma. Figure 7 shows the solar-wind ion temperatures as a function of the absolute value of the total compressible energy cascade rates. The color bar corresponds to the mean solar-wind velocity. The observational results show a clear trend between the increase in the solar-wind temperature and the compressible cascade rates. As the cascade rate increases, we observe a slight increase in the local ion temperature. Note that the temperature increase also correlated with the increase in the solar-wind speed.

Zoom In Zoom Out Reset image size

To the best of our knowledge, we reported the first estimation of the compressible energy transfer rate in the solar wind at ∼0.2–0.4 au and 1.5–1.7 au, using PSP and MAVEN observations, respectively. As we got closer to the Sun, we observed a slight increase in the solar wind's density, velocity, and magnetic fluctuations (Bruno & Carbone 2005). For the most compressible events (when the compressibility is ∼25%), in PSP's first perihelion (on 2018 November 6), we observed that the compressible energy cascade rate is at least 5 orders of magnitude larger than the cascade at ∼1.6 au. Also, we observed a moderate growth in the compressible cascade rate with respect to the incompressible rate. The increases in the compressible (and incompressible) cascade rates when we compared the results between 0.2 au and 1.7 au are mainly due to an increase in the magnetic and the number density fluctuation levels (see, Figures 1 and 4). Previous studies showed that the density fluctuations would increase considerably the cascade rate when the compressibility is larger than 30% or when it enters into sub-ion scales (see, Hadid et al. 2017; Andrés et al. 2018, 2019). Our results are compatible with these findings at the MHD scales studied here. Also, Adhikari et al. (2020b) computed the frequency distributions of solar wind compressibility between 0.17 au and 0.61 au using PSP observations. Their results showed that density fluctuations (normalized to the mean density) are concentrated mainly around 0.15, which is compatible with the density fluctuation levels found in this study. In contrast, recent computations of the magnetic compressibility coefficient, C_B ≡ (δ∣ B ∣/∣δ B ∣)², at 0.17 au have shown a clear decrease toward smaller heliocentric distances (at all frequencies) (see, Chen et al. 2020). In particular, the authors observed that magnetic compressibility levels at PSP's perihelion are an order of magnitude smaller than at 1 au. However, this decrease in the magnetic compressibility shows no overall impact on the density fluctuation level, and consequently on the compressible cascade rate estimations. The expression for ε_C in Equation (5) is a fully nonlinear equation derived for compressible isothermal MHD turbulence, which is valid in the inertial range under the typical assumptions for fully developed turbulence (e.g., Frisch 1995). As such, it naturally involves the dynamics of compressible magnetosonic (i.e., fast/slow) and entropy (i.e., zero temporal frequency) modes. Therefore, if these modes exist in our data sets, they are naturally captured by our formalism. The moderate enhancement in the cascade rate is due only to the relatively low density fluctuations observed (as is generally the case in the solar wind), which means that magnetosonic-like modes would represent only a small fraction of the total solar-wind fluctuations. Which of the fast or slow (or entropy) modes is dominant in our data sets is beyond the scope of the present work, although there seems to be a consensus that the latter carry the major part of the density fluctuations (e.g., Klein et al. 2014; Andrés et al. 2017).

Our findings here for different heliocentric distances confirm that density fluctuations in the solar wind amplify the cascade rate with respect to the incompressible model (Sorriso-Valvo et al. 2007; Andrés et al. 2020; Bandyopadhyay et al. 2020). When the density fluctuations are relatively larger than ∼15%, our results showed that the leading role in the amplification of the cascade rate is due to a competition between the new compressible flux component, 〈∣ε_2C∣〉, and the compressible Yaglom-like generalization, 〈∣ε_1C∣〉, in the exact relation. As the plasma compressibility locally increases through the heliosphere we observed two clear features: on one hand, the Yaglom-like component, ∣ε_1C∣, has a larger spread (larger and/or smaller values) around the incompressible value. While this feature could be due to the relatively poor statistics in the THEMIS- and MAVEN-analyzed data sets, the spread could be due to the large plasma density fluctuation levels present in the PSP data set. On the other hand, as the plasma compressibility increases we observed that the compressible component, ε_2C, becomes of the order of the Yaglom-like component in the inertial range. Andrés et al. (2018) computed all the terms in Equation (5) using DNSs of sub-sonic isothermal compressible MHD turbulence. The authors found that even for large guide field values, the energy cascade rate, ε_C, is mainly due to the flux terms. Also, they observed that ε_2C increases by at least one order of magnitude when the plasma compressibility increased by a factor of 2 in the system. Furthermore, they found that when compressibility increases in the plasma, the compressible component becomes dominant in energy transfer. These numerical finding are in agreement with our observational results in Figures 5 and 6. We speculate that this trend should be observed in more compressible plasma environments, like the Earth's magnetosheath (Hadid et al. 2017; Huang et al. 2017; Andrés et al. 2019; Li et al. 2020).

Our observational results show a correlation between the solar-wind temperature and the compressible energy cascade rate, where the larger the cascade rate, the higher the temperature. Also, a clear increase in the temperature is connected with fast solar-wind speeds. The relation between the cascade rate and the local temperature in the solar wind has been studied previously in the literature (see, Vasquez et al. 2007; MacBride et al. 2008; Marino et al. 2008; Banerjee et al. 2016). As we discussed above, the turbulent cascade rate has been proposed as a solution to the solar wind heating problem (e.g., Richardson et al. 1995; Matthaeus & Velli 2011). Our observational results here support the idea that compressible turbulent cascades may heat the plasma at different heliocentric distances. We emphasize that all results reported here were obtained inside the solar-wind ecliptic plane. A natural extension of this study would be the use of ESA's recently launched Solar Orbiter mission (Müller et al. 2020) data. Its observations should give us access to the polar regions of the solar wind within heliocentric distance never explored before.

Finally, two limitations should be kept in mind when interpreting the results above. The first one concerns the isothermal closure used to derive our theoretical model. Although the actual "thermodynamics" of the solar wind is an open and hotly debated question (because of the collisionless and magnetized nature of its plasma), the isothermal closure may not be always applicable to solar-wind observations (e.g., Totten et al. 1995; Kartalev et al. 2006; Bhattacharjee et al. 2009; Nicolaou et al. 2014; Verscharen et al. 2019; Nicolaou et al. 2020; Adhikari et al. 2020a). For instance, Nicolaou et al. (2020) analyzed solar-wind proton plasma measurements obtained by PSP to derive the effective polytropic indexes of large and small-scale fluctuations. The authors found that the large-scale variations of the solar-wind proton density and temperature follow a polytropic model with a polytropic index of ∼5/3 and the short-scale fluctuations follow a model with even larger polytropic indexes. A natural extension of our isothermal model is the use of an exact polytropic MHD law. This has been achieved recently by Simon & Sahraoui (2021), who extended the existing turbulence exact laws to isentropic flows (i.e., constant entropy) that encompasses isothermal and polytropic closures. The authors showed that the polytropic law does not impact significantly the cascade rate estimate in the solar wind (based only on flux terms) as long as any density fluctuations remain moderate. A more realistic description of the solar wind within the same formalism of exact laws would be to include entropy variations and/or pressure anisotropy due to the presence of a mean magnetic field. These exact relations are currently unavailable.

The second limitation of the present study is the use of the isotropic version of the exact law. Since PSP is a single spacecraft mission, spatial variations of the turbulent fluctuations can be accessed only along a single direction given by the mean flow speed based on the Taylor assumption. Using isotropy along with the Taylor hypothesis means that observations inferred from sampling the direction of the mean flow is representative of the actual 3D properties of the energy cascade rate. This inevitable limitation obviously comes up against the anisotropic nature of turbulence in magnetized plasmas. Nevertheless, "reduced" forms of anisotropy of the cascade rate have been explored observationally by examining the dependence of the incompressible cascade rate upon the angle θ_BV (MacBride et al. 2008). A similar approach has been adopted for the compressible cascade in the solar wind at 1 au (Hadid et al. 2017) and the terrestrial magnetosheath, where density fluctuations were found to enhance its spatial anisotropy (Hadid et al. 2017). In all these observational studies based on the Taylor assumption it is however important to select only intervals where the angle θ_BV is approximately constant to guaranty sampling the same direction of space, and thus provide a more robust estimate of the cascade rate, as we did in the present work. A general 3D estimate of the cascade rate free from the Taylor hypothesis requires using multi-spacecraft data and developing an appropriate method for that. Two attempts were done in that direction using data from the multi-spacecraft missions Cluster (Osman et al. 2011) and MMS (Andrés et al. 2019), but which still involved the Taylor hypothesis.

N.A. acknowledges financial support from the following grants: PICT 2018 1095 and UBACyT 20020190200035BA. N.A., L.H.Z., and F.S. acknowledge financial support from CNRS/CONICET Laboratoire International Associé (LIA) MAGNETO. N.R. is an Assistant Research Scientist at NASA Goddard Space Flight Center, hired by MAVEN Project Scientist Support, and administered by CRESST II and UMBC. We thank the entire PSP, THEMIS, and MAVEN team and instrument leads for data access and support. PSP data are publicly available (https://research.ssl.berkeley.edu/data/psp/). The THEMIS/ARTEMIS data come from the AMDA database (http://amda.cdpp.eu/). MAVEN data are publicly available through the Planetary Data System (https://pds-ppi.igpp.ucla.edu/index.jsp).