The Transport and Evolution of MHD Turbulence throughout the Heliosphere: Models and Observations

Abstract

A detailed study of solar wind turbulence throughout the heliosphere in both the upwind and downwind directions is presented. We use an incompressible magnetohydrodynamic (MHD) turbulence model that includes the effects of electrons, the separation of turbulence energy into proton and electron heating, the electron heat flux, and Coulomb collisions between protons and electrons. We derive expressions for the turbulence cascade rate corresponding to the energy in forward and backward propagating modes, the fluctuating kinetic and magnetic energy, the normalized cross-helicity, and the normalized residual energy, and calculate the turbulence cascade rate from 0.17 to 75 au in the upwind and downwind directions. Finally, we use the turbulence transport models to derive cosmic ray (CR) parallel and perpendicular mean free paths (mfps) in the upwind and downwind heliocentric directions. We find that turbulence in the upwind and downwind directions is different, in part because of the asymmetric distribution of new born pickup ions in the two directions, which results in the CR mfps being different in the two directions. This is important for models that describe the modulation of cosmic rays by the solar wind.

1. Introduction

Turbulence in solar wind plays an important role in many aspects of plasma behavior: (i) the observed nonadiabatic radial profile of the solar wind temperature [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23]; (ii) the scattering of the solar energetic particles [24,25]; (iii) the heating of the coronal plasma to millions of degrees Kelvin and the acceleration of the solar wind from a subsonic to a supersonic speed [26,27,28,29,30,31,32,33,34,35,36,37,38,39,40], (iv) the transport of cosmic rays throughout the heliosphere [22,41,42,43,44,45], and (v) the increase in the solar wind proton and electron entropy [20,46,47].

Since the early space age observations of solar wind turbulence [48,49,50,51,52,53,54,55,56,57], interplanetary fluctuations have been described either as evidence of a turbulent solar wind [48] or as the superposition of forward propagating Alfvén waves [50]. These early results are contradictory. The turbulent characteristics of the solar wind are based on the idea that the solar wind plasma and magnetic field fluctuations exhibit a power-law like spectra as predicted by a Kolmogorov [58] or Kraichnan [59] phenomenology. However, the strongly Alfvénic character of the solar wind fluctuations observed by Belcher and Davis [50] suggests the existence of pure Alfvén waves, and does not, therefore, support nonlinear interactions between counter-propagating Alfvén waves. In the latter case, the observed unidirectional propagation of Alfvén waves that exhibit a Kolmogorov-like power law [60,61] has been ascribed to the interaction with dominant quasi-2D turbulence [62]. Coleman [48] suggested that sheared solar wind flows that are generated by the difference in fast and slow solar wind speeds, generates low-frequency MHD turbulence in the inner heliosphere. Tu et al. [63] proposed a model to explain the radial evolution of the power spectrum of the Alfvénic fluctuations observed by $Helios$ 1 and 2 in the inner heliosphere between 0.3 and 1 au [51]. In their model, they assumed that most of the fluctuations correspond to forward propagating incompressible (Alfvénic) wave modes, with fewer of the fluctuations corresponding to backward propagating wave modes. This indicates that there is a weak interaction between counter-propagating modes, which is consistent with Coleman [48] in the sense that the turbulent shear source generates backward propagating wave modes [17,64].

The modeling of the transport of incompressible interplanetary fluctuations began with the Wentzel–Kramers–Brillouin (WKB) theory. WKB theory describes the transport of linear, non-interacting MHD Alfvén waves in an inhomogeneous background [65,66,67,68,69]. The description incorporates convection, wave propagation, and expansion effects, but does not consider the nonlinear interaction between forward and backward propagating Alfvén waves. Forward and backward propagating Alfvén waves interact on the nonlinear timescales ${\tau }_{nl}^{\pm }\sim {\lambda }^{\pm }/{\langle {z^{\mp }}^2\rangle }^{1/2}$, where ${\lambda }^{\pm }$ are correlation lengths corresponding to the Elsässer energies $\langle {z^{\pm }}^2\rangle$. Similarly, the propagation of Alfvén waves introduces the Alfvénic timescales ${\tau }_A^{\pm }\sim {\lambda }_{\left|\right|}^{\pm }/\left|{\mathbf{V}}_{A0}\right|$, where ${\lambda }_{\left|\right|}^{\pm }$ are the correlation lengths along the direction of the magnetic field, and $|{\mathbf{V}}_{A0}|$ is the magnitude of the Alfvén velocity. If ${\tau }_{nl}^{\pm }>{\tau }_A^{\pm }$, the counter-propagating modes experience many interactions before their energy is transferred into smaller scales [59,70]. This is known as a weak turbulence. If ${\tau }_{nl}^{\pm }<{\tau }_A^{\pm }$, the Alfvénic fluctuations undergo multiple decays before they interact. This is known as a strong turbulence. When the two timescales are comparable, i.e., ${\tau }_{nl}^{\pm }\sim {\tau }_A^{\pm }$, this leads to the scaling of $k_{\left|\right|}\propto k_{\perp }^{2/3}$, and is known as a “critical balanced” state [71].

The outline of the paper is as follows. In Section 1, we present an overview of incompressible MHD turbulence transport models, and in Section 2, we discuss 3D solar wind turbulence models. In Section 3, we discuss 1D steady-state turbulence models, including the incorporation of electrons, and present comparisons between the theoretical results with those measured by $PSP$, $Helios$ 2, $Voyager$ 2, $New$$Horizons$$Solar$$Wind$$around$$Pluto$ ($NH$$SWAP$) in the upwind direction, and with $Pioneer$ 10 in the downwind direction. In Section 4, we derive the turbulence cascade rate for the energy in forward and backward propagating modes, the total turbulent energy, the cross-helicity, the residual energy, the turbulent kinetic energy, and the turbulent magnetic energy, and discuss the results. Section 5 discusses the cosmic rays diffusion tensor and its relation to the turbulence transport models. Finally, further discussion and conclusions are found in Section 6.

2. Incompressible MHD Turbulence Models

The fluctuating Elsässer variables ${\mathbf{z}}^{\pm }=\mathbf{u}\pm \mathbf{b}/\sqrt{{\mu }_0\rho }$ [72] (where $\mathbf{u}$ is the fluctuating solar wind speed, $\mathbf{b}$ the fluctuating the magnetic field, $\rho$ the proton mass density, and ${\mu }_0$ the magnetic permeability) define the forward and backward propagating modes [63,73,74,75,76,77,78]. The incompressible MHD transport equations can be expressed in terms of the Elsässer variables as [73,74,75],

$$\begin{array}{cc}\hfill & \frac{\partial {\mathbf{z}}^{\pm }}{\partial t}+(\mathbf{U}\mp {\mathbf{V}}_A)·\nabla {\mathbf{z}}^{\pm }+\frac{1}{2}\nabla ·\left(\frac{\mathbf{U}}{2}\pm {\mathbf{V}}_A\right){\mathbf{z}}^{\pm }+{\mathbf{z}}^{\mp }·[\nabla \mathbf{U}\pm \frac{\nabla \mathbf{B}}{\sqrt{{\mu }_0\rho }}\hfill \\ \hfill & -\frac{1}{2}I\nabla ·\left(\frac{\mathbf{U}}{2}\pm {\mathbf{V}}_A\right)]={\mathbf{NL}}_{\pm }+S^{\pm },\hfill \end{array}$$

(1)

where I is the identity matrix, ${\mathbf{NL}}_{\pm }$ the nonlinear terms, $\mathbf{U}$ the solar wind speed, ${\mathbf{V}}_A$ the Alfvén velocity, and $S^{\pm }$ are sources of turbulence. When constructing ensemble-averaged moments of the Elsässer variables using (1), such as $\langle {\mathbf{z}}^{+}·{\mathbf{z}}^{+}\rangle$ or $\langle {\mathbf{z}}^{-}·{\mathbf{z}}^{-}\rangle$ etc., terms inside the squared brackets are introduced like $\langle z_i^{+}{z}_j^{-}\rangle$, which is an MHD analog of the hydrodynamic Reynolds stress. Due to the presence of the $\langle z_i^{+}{z}_j^{-}\rangle$ terms, the backward and forward propagating modes interact through the small-scale fluctuations, and large-scale solar wind speed and magnetic field.

The advantage of the transport Equation (1) is that it includes the nonlinear interaction between ${\mathbf{z}}^{+}$ and ${\mathbf{z}}^{-}$, and introduces the cross-helicity and the residual energy (definitions given below), which are measures of the Alfvénicity of the solar wind turbulence. The turbulence transport model equations are, therefore, useful to study the radial evolution of MHD turbulence, Alfvénicity, and applications, such as coronal heating, and the solar wind heating throughout the heliosphere.

Starting from Equation (1), Zank et al. [76] developed a simple model comprising two coupled turbulence transport equations for the turbulent magnetic energy density $E_b$ (=$\langle b^2\rangle /{\mu }_0\rho$, where $\langle b^2\rangle$ is the fluctuating magnetic energy) and the correlation length of the magnetic field fluctuations ${\lambda }_b$. The Zank et al. [76] model was compared successfully with turbulent magnetic energy observations measured by $Voyager$ 1, 2 and $Pioneer$ 11 between 1 and 40 au.

Matthaeus et al. [4], Smith et al. [5], and Smith et al. [6] included the dissipation of magnetic energy [76] in the proton temperature equation to study proton heating throughout the heliosphere [3,8,9]. Matthaeus et al. [4], Smith et al. [5], and Smith et al. [6] found that the radial profile of the solar wind temperature is similar to that observed by $Pioneer$ 11 and $Voyager$ 2, including the observed temperature increase in the distant heliosphere. Of particular note was an observed temperature increase [1] in the solar wind temperature beyond ∼25–30 au. This was explained by Williams et al. [3] as being due to the dissipation of PUI-driven waves in the outer heliosphere. Adhikari et al. [16] extended the Zank et al. [76] model by using time-dependent boundary conditions and time-dependent turbulent sources. They compared their model results with $Voyager$ 2 measurements, and found periodic changes in the magnetic energy density, correlation length, and the solar wind temperature, especially beyond $\sim 10$ au.

In developing an incompressible MHD turbulence model, Matthaeus et al. [79] proposed that the radial evolution of the normalized cross-helicity vary [80], but also assumed that the normalized residual energy (the energy difference between the fluctuations in the velocity and magnetic field) is constant. Later, Breech et al. [10] extended the turbulence transport model to study the evolution of turbulence energy, normalized cross-helicity, and the turbulent correlation length from the inner to outer heliosphere by assuming that the normalized residual energy is constant, and the solar wind is super-Alfvénic, i.e., $\left|\mathbf{U}\right|\gg |{\mathbf{V}}_A|$. Therefore, the Breech et al. [10] model cannot be used for sub-Alfvénic flows, i.e., $\left|\mathbf{U}\right|\ll |{\mathbf{V}}_A|$, where the Alfvén velocity is larger than the solar wind speed.

The incompressible MHD turbulence model was further developed by Zank et al. [77] to include the residual energy, various correlation lengths, and the Alfvén velocity, making it suitable for sub-Alfvénic flows. Zank et al [77] used Equation (1) to develop a six coupled equations turbulence model that describes the transport of the total energy density, cross-helicity, residual energy, and correlation lengths corresponding to forward and backward propagating modes and the residual energy. To model the dissipative terms for the various transport equations, Zank et al. [77] used a single-point closure related to that of von Kármán and Howarth [81], and further developed [10,76,82,83,84]. The six coupled turbulence transport equations of Zank et al. [77] describe the transport of turbulence in inhomogeneous sub-Alfvénic and super-Alfvénic flows. The Zank et al. [77] turbulence equations are quasi-linear in the spatial evolution operators, and nonlinear in the dissipation terms.

Adhikari et al. [17] studied the evolution of turbulence in the super-Alfvénic solar wind flow from 0.3 to 75 au using the Zank et al. [77] turbulence model and $Voyager$ 2 measurements. They compared their model results with $Helios$ 2, $Ulysses$, and $Voyager$ 2 measurements, and found that (i) shear driving is responsible for the in situ generation of backward propagating modes; (ii) the correlation lengths corresponding to forward and backward propagating modes increase with distance, and are approximately equal beyond ∼30 au, and (iii) the fluctuating kinetic and magnetic energies are approximately equal beyond ∼30 au.

Adhikari et al. [20] developed a conservative formulation of the coupled solar wind and incompressible MHD turbulence transport model equations using the Zank et al. [77] turbulence model. They showed that not only the pickup ions, but also stream-shear leads to a decrease in the solar wind speed. They showed that the sum of the solar wind flow (kinetic plus enthalpy) energy and turbulent (magnetic) energy is constant, illustrating that kinetic solar wind energy is transferred into turbulent energy via stream-shear and pickup ion isotropization, which then in turn heats the solar wind via the dissipation of turbulence.

Consider the following moments of the Elsässer variables,

$$\begin{array}{cc}\hfill & E_T\equiv \frac{\langle {\mathbf{z}}^{+}·{\mathbf{z}}^{+}\rangle +\langle {\mathbf{z}}^{-}·{\mathbf{z}}^{-}\rangle }{2}=\frac{\langle {z^{+}}^2\rangle +\langle {z^{-}}^2\rangle }{2}=\langle u^2\rangle +〈\frac{b^2}{{\mu }_0\rho }〉;\hfill \\ \hfill & E_C\equiv \frac{\langle {\mathbf{z}}^{+}·{\mathbf{z}}^{+}\rangle -\langle {\mathbf{z}}^{-}·{\mathbf{z}}^{-}\rangle }{2}=\frac{\langle {z^{+}}^2\rangle -\langle {z^{-}}^2\rangle }{2}=2〈\mathbf{u}·\frac{\mathbf{b}}{\sqrt{{\mu }_0\rho }}〉;\hfill \\ \hfill & E_D\equiv \langle {\mathbf{z}}^{+}·{\mathbf{z}}^{-}\rangle =\langle u^2\rangle -〈\frac{b^2}{{\mu }_0\rho }〉,\hfill \end{array}$$

(2)

where $E_T$ is twice the total energy of the fluctuations per unit mass (or the sum of the square of the fluctuating plasma and Alfvénic velocity), $E_C$ is the cross-helicity of the fluctuations (or a measure of the obliquity of the fluctuating velocity and magnetic field vectors), and $E_D$ is twice the difference between the fluctuating kinetic and magnetic energy per unit mass. The variable $\langle {z^{\pm }}^2\rangle$ is the energy in forward/backward propagating modes. Similarly, the following relations hold [77,85],

$$\begin{array}{cc}\hfill & r_A\equiv \frac{\langle u^2\rangle }{E_b}=\frac{E_T+E_D}{E_T-E_D}=\frac{1+{\sigma }_D}{1-{\sigma }_D};\langle {z^{+}}^2\rangle =E_T+E_C;\langle {z^{-}}^2\rangle =E_T-E_C;\hfill \\ \hfill & {\sigma }_c=\frac{E_C}{E_T};{\sigma }_D=\frac{E_D}{E_T};\langle u^2\rangle =\frac{E_T+E_D}{2};E_b=〈\frac{b^2}{{\mu }_0\rho }〉=\frac{E_T-E_D}{2};\hfill \\ \hfill & {\lambda }_u=\frac{\langle {z^{+}}^2\rangle {\lambda }^{+}+\langle {z^{-}}^2\rangle {\lambda }^{-}+E_D{\lambda }_D}{4\langle u^2\rangle };{\lambda }_b=\frac{\langle {z^{+}}^2\rangle {\lambda }^{+}+\langle {z^{-}}^2\rangle {\lambda }^{-}-E_D{\lambda }_D}{4E_b},\hfill \end{array}$$

(3)

where $r_A$ is the Alfvén ratio, ${\sigma }_c$ the normalized cross-helicity, ${\lambda }^{\pm }$ the correlation length corresponding to forward/backward propagating modes, ${\lambda }_D$ the correlation length corresponding to the residual energy, and ${\lambda }_{u,b}$ the correlation lengths corresponding to the velocity/magnetic field fluctuations.

Similarly, we can also introduce the following relations,

$$\begin{array}{cc}\hfill & L^{\pm }\equiv \int \langle {\mathbf{z}}^{\pm }·{{\mathbf{z}}^{\pm }}^{\prime }\rangle dy\equiv \langle {z^{\pm }}^2\rangle {\lambda }^{\pm }=(E_T\pm E_C){\lambda }^{\pm };\hfill \\ \hfill & L_D\equiv \int \langle {\mathbf{z}}^{+}·{{\mathbf{z}}^{-}}^{\prime }+{{\mathbf{z}}^{+}}^{\prime }·{\mathbf{z}}^{-}\rangle dy\equiv E_D{\lambda }_D.\hfill \end{array}$$

(4)

where, $L^{\pm }$ and $L_D$ denote the correlation functions corresponding to forward/backward propagating modes, and the residual energy. The prime denotes the spatial lagged Elsässer variables separated by the distance y.

The detailed derivation of the six coupled turbulence transport equations is presented in Zank et al. [77], and the equations are given by,

$$\begin{array}{cc}\hfill \frac{\partial E_T}{\partial t}+& \mathbf{U}·\nabla E_T+\frac{1}{2}\nabla ·\mathbf{U}{E}_T-{\mathbf{V}}_A·\nabla E_C+\nabla ·{\mathbf{V}}_A{E}_C+\left(2a-\frac{1}{2}\right)\nabla ·\mathbf{U}{E}_D\hfill \\ \hfill & -2a{n}_i{n}_j\frac{\partial U_j}{\partial x_i}{E}_D=-\frac{(E_T+E_C){|E_T-E_C|}^{1/2}}{{\lambda }^{+}}-\frac{(E_T-E_C){|E_T+E_C|}^{1/2}}{{\lambda }^{-}}\hfill \\ \hfill & +\langle {\mathbf{z}}^{+}·{\mathbf{S}}^{+}\rangle +\langle {\mathbf{z}}^{-}·{\mathbf{S}}^{-}\rangle ;\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill \frac{\partial E_C}{\partial t}& +\mathbf{U}·\nabla E_C+\frac{1}{2}\nabla ·\mathbf{U}{E}_C-{\mathbf{V}}_A·\nabla E_T+\nabla ·{\mathbf{V}}_A{E}_T+(2b-1)\nabla ·{\mathbf{V}}_A{E}_D\hfill \\ \hfill & -2b{n}_i{n}_j\frac{\partial V_{Aj}}{\partial x_i}{E}_D+\frac{b}{\rho }\left({\mathbf{V}}_{\mathbf{A}}·\nabla \rho -n_i{n}_j{V}_{Aj}\frac{\partial \rho }{\partial x_i}\right)E_D=-\frac{(E_T+E_C){|E_T-E_C|}^{1/2}}{{\lambda }^{+}}\hfill \\ \hfill & +\frac{(E_T-E_C){|E_T+E_C|}^{1/2}}{{\lambda }^{-}}+\langle {\mathbf{z}}^{+}·{\mathbf{S}}^{+}\rangle -\langle {\mathbf{z}}^{-}·{\mathbf{S}}^{-}\rangle ;\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill \frac{\partial E_D}{\partial t}& +\mathbf{U}·\nabla E_D+\frac{1}{2}\nabla ·\mathbf{U}{E}_D+\left(2a-\frac{1}{2}\right)\nabla ·\mathbf{U}{E}_T-(2b-1)\nabla ·{\mathbf{V}}_A{E}_C\hfill \\ \hfill & +\frac{1}{\sqrt{E_T^2-E_C^2}}(E_C{\mathbf{V}}_A·\nabla E_T-E_T{\mathbf{V}}_A·\nabla E_C)-\frac{b}{\rho }\left[{\mathbf{V}}_A·\nabla \rho -n_i{n}_j{V}_{Aj}\frac{\partial \rho }{\partial x_i}\right]E_C\hfill \\ \hfill & -2n_i{n}_j\left(a\frac{\partial U_j}{\partial x_i}{E}_T-b\frac{\partial V_{Aj}}{\partial x_i}{E}_C\right)=-E_D\left(\frac{|E_T-E_C{|}^{1/2}}{{\lambda }^{+}}+\frac{|E_T+E_C{|}^{1/2}}{{\lambda }^{-}}\right)\hfill \\ \hfill & +\langle {\mathbf{z}}^{-}·{\mathbf{S}}^{+}\rangle +\langle {\mathbf{z}}^{+}·{\mathbf{S}}^{-}\rangle ;\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill & \frac{\partial L^{\pm }}{\partial t}+(\mathbf{U}\mp {\mathbf{V}}_A)·\nabla L^{\pm }+\nabla ·\left(\frac{\mathbf{U}}{2}\pm {\mathbf{V}}_A\right)L^{\pm }+\left(a-\frac{1}{4}\right)\nabla ·\mathbf{U}{L}^D\pm \left(b-\frac{1}{2}\right)\hfill \\ \hfill & \times \nabla ·{\mathbf{V}}_A{L}^D\pm \frac{b}{2\rho }{\mathbf{V}}_A·\nabla \rho L^D-n_i{n}_j\frac{\partial }{\partial x_i}(a{U}_j\pm b{V}_{Aj})L^D\mp \frac{b}{2\rho }{n}_i{n}_j{V}_{Aj}\frac{\partial \rho }{\partial x_j}{L}^D=0;\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill & \frac{\partial L_D}{\partial t}+\mathbf{U}·\nabla L^D+\frac{1}{2}\nabla ·\mathbf{U}{L}^D+{\left(\frac{L^{+}}{L^{-}}\right)}^{1/2}{\mathbf{V}}_A·\nabla L^{-}-{\left(\frac{L^{-}}{L^{+}}\right)}^{1/2}{\mathbf{V}}_A·\nabla L^{+}\hfill \\ \hfill & \left[\left(2a-\frac{1}{2}\right)\nabla ·\mathbf{U}-2a{n}_i{n}_j\frac{\partial U_j}{\partial x_i}\right](L^{+}+L^{-})-\left[(2b-1)\nabla ·{\mathbf{V}}_A-2a{n}_i{n}_j\frac{\partial V_{Aj}}{\partial x_i}\right]\hfill \\ \hfill & \times (L^{+}-L^{-})-\frac{b}{\rho }\left({\mathbf{V}}_A·\nabla \rho -n_i{n}_j{V}_{Aj}\frac{\partial \rho }{\partial x_j}\right)(L^{+}-L^{-})=0,\hfill \end{array}$$

(9)

where a and b are the structural similarity parameters related to the large-scale velocity $\mathbf{U}$, and the Alfvén velocity ${\mathbf{V}}_A$, respectively. The parameters a and b can be chosen freely, and choosing $a=1/2$ and $b=1/2$ or $a=1/3$ and $b=1/3$ produces 2D or 3D mixed tensors [83]. In Equations (5)–(7), the first three terms on the lhs are similar to the standard WKB model, while the remaining terms describe the indirect mixing and coupling of forward and backward propagating modes. The first two terms on the rhs are the nonlinear terms, and the latter two terms indicate turbulence sources. The nonlinear dissipative term for the total turbulent energy is derived by following a Kolmogorov phenomenology, and has been verified by numerical simulations [86]. The nonlinear terms can also be expressed using instead an Iroshnikov–Kriachnan phenomenology e.g., [13].

3. 3D Solar Wind–Incompressible MHD Turbulence Model

A global turbulence model has been developed by coupling various forms of the 3D turbulence transport model with a global MHD model of the solar wind [14,21,22,45,87,88,89,90,91,92,93]. Usmanov et al. [14] developed a self-consistent 3D solar wind-turbulence model to study the effect of turbulence on large-scale solar wind properties and solar wind heating from 0.3 to 100 au. Usmanov et al. [87] extended their three-dimensional solar wind model to incorporate turbulence transport and pickup ions as a separate fluid. The global MHD model is based on solving the large-scale Reynolds-averaged MHD equations when coupled to pickup ions and the associated generation of turbulence from 0.3 to 100 au. They confirmed well-known effects of pickup ions on the large-scale solar wind, such as the slowing down of the solar wind speed due to the momentum loss from the solar wind protons to pickup ions, and the compressed azimuthal magnetic field in the outer heliopshere [94]. They also found that the decelerated solar wind speed and the pickup ions weaken corotating interaction regions (CIRs), confirming results by Zank et al. [95] and Rice and Zank [96]. Usmanov et al. [88] further extended the 3D MHD solar wind and turbulence model to include eddy viscosity, turbulent resistivity, and turbulent heating. Usmanov et al. [88] considered a system of co-moving solar wind protons, electrons and interstellar pickup ions, each species having a separate energy equation. They presented three-dimensional distributions of the mean-flow and turbulence parameters throughout the heliosphere under the given boundary conditions at the coronal base, finding that the effect of the eddy viscosity, and velocity shear on the mean-flow parameters increases the solar wind proton and electron temperatures. Usmanov et al. [89] further extended their 3D solar wind turbulence model to include electron heat conduction, radiative cooling, Coulomb collisions, Reynolds stresses, eddy viscosity, and turbulent heating of protons and electrons. Chhiber et al. [45] used the 3D solar wind-turbulence model of Usmanov et al. [14] to study the cosmic ray diffusion tensor throughout the heliosphere [41].

In a related study, Kryukov et al. [90] developed a 3D solar wind turbulence model that incorporates pickup ions and solar wind heating over a heliocentric distance 1–80 au. Like the 3D MHD models of Usmanov et al. [14,87,88,89], Kryukov et al. [90], and Chhiber et al. [93] incorporate the back reaction of turbulence on the global MHD flow by including both the turbulent heating term and Reynolds-averaged turbulence terms in the MHD model. However, these papers combined the simpler turbulence model [10,76] with 3D MHD solar wind models. Similarly, they assumed a constant normalized residual energy in the turbulence model.

Wiengarten et al. [21] developed a 3D model in which they coupled turbulence transport equations with the Reynolds-averaged ideal MHD equations, and applied their model to study the propagation of coronal mass ejections (CMEs). Wiengarten et al. [21] introduced the effect of turbulence driven by shear, which was not taken into account in Usmanov et al. [14]. They found that a CME enhances turbulence levels and reduces the cross-helicity, but the overall effect of CMEs is negligible. Wiengarten et al. [22] coupled a two-component turbulence model of Oughton et al. [15] self-consistently with the large-scale MHD equations, and calculated the mean-free path and the drift lengthscale of cosmic rays throughout the heliosphere.

Finally, Shiota et al. [91] coupled a more sophisticated turbulence model, Equations (5)–(9) of Zank et al. [77], with a large-scale 3D MHD solar wind model [97] to study the basic physical changes of turbulent transport in a more complex structured solar wind. The model included the inhomegeneity of the background solar wind and IMF, and sources of turbulence. To obtain a simple structure for the solar wind, Shiota et al. [91] assumed two steady magnetic field configurations, such as an outward monopole magnetic field whose intensity is spherically symmetric at $r=2.5$ R${}_{\odot }$, and a tilted dipole magnetic field at an angle 30${}^{\circ }$ from the rotational axis of the Sun. Shiota et al [91] also considered (i) a spherically symmetric constant solar wind speed of 400 kms${}^{-1}$ or 600 kms${}^{-1}$, and (ii) a bimodal (fast and slow) solar wind flow. In the latter case, the speed is determined by the titled dipole potential field and the Wang–Sheely–Arge 2000 [98,99] model for coronal magnetic fields. Shiota et al [91] found that the non-axisymmetric solar wind speed distribution with a spherically symmetric IMF enhances the energy in backward propagating modes and the normalized residual energy in the interaction region.

4. Incompressible MHD Turbulence Model: Proton and Electron Heating

Most of the turbulence transport models include proton heating [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23], and neglect electron heating. Electrons carry about half of the thermal energy of the plasma, which indicates that the role of electrons should not be neglected. Electron effects includes the decomposition of turbulence energy into proton and electron heating, Coulomb collisions between electrons and protons, and electron heat flux.

The study by Leamon et al. [100] found that about 60% of the turbulent energy heats solar wind protons at 1 au. Using 60% of the turbulence energy to heat solar wind protons and 40% of the turbulence energy to heat solar wind electrons, Breech et al. [11] developed a turbulence model that included electrons [11,47,93,101,102,103]. They found that the solar wind proton and electron temperatures decrease slower than that expected of an adiabatic radial profile, viz. $r^{-4/3}$. Recently, Boldyrev et al. [104] found that according to their kinetic theory, the electron temperature decreases as $r^{-2/5}$. The electron heat flux mediated by wave particle interactions between electrons and whistler waves [105] can also influence the radial profiles of the electron and proton temperatures.

Adhikari et al. [47] developed a model that coupled the solar wind equations with the nearly incompressible MHD turbulence transport equations, including the effects of electrons. Their model produced theoretical solar wind proton and electron temperatures that resemble measurements made by PSP and Helios 2. Similarly, they found that the electron and proton entropy increases by about 3% and 2.55% over the heliocentric distance 45.5–215 R${}_{\odot }$.

In this paper, we use a 1D form of the Equations (5)–(9) to describe the steady-state transport of turbulent energy in forward and backward propagating modes $\langle z^{\pm 2}\rangle (=E_T\pm E_C)$, the residual energy $E_D$, and the corresponding correlation functions $L^{\pm }$ and $L_D$. We neglect the mixing terms in Equations (5)–(9), and assume a spherically symmetric constant solar wind speed U. The 1D steady-state turbulence transport equations can then be written as [17,77],

$$\begin{array}{cc}\hfill & (U\mp V_A)\frac{d\langle {z^{\pm }}^2\rangle }{dr}+\left(\frac{U}{r}\pm \frac{d{V}_A}{dr}\pm \frac{2V_A}{r}\right)\langle {z^{\pm }}^2\rangle +2\left(2a-\frac{1}{2}\right)\frac{U}{r}{E}_D\pm (2b-1)\hfill \\ \hfill & \times \left(\frac{d{V}_A}{dr}+\frac{2V_A}{r}\right)E_D\pm b\frac{V_A}{\rho }\frac{d\rho }{dr}=-2\frac{\langle {z^{\pm }}^2\rangle {\langle {z^{\mp }}^2\rangle }^{1/2}}{{\lambda }^{\pm }}+2C_{sh}^{\pm }\frac{r_0|\Delta U|V_{A0}^2}{r^2}\hfill \\ \hfill & +\frac{f_D{n}_H^{\infty }U{V}_A^{\prime }}{n_{sw}^0{\tau }_{ion}^0}exp\left(\frac{-\zeta L}{rsin\zeta }\right);\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill & U\frac{d{E}_D}{dr}+\frac{U}{r}{E}_D+2\left(2a-\frac{1}{2}\right)\frac{U}{r}\frac{\langle {z^{+}}^2\rangle +\langle {z^{-}}^2\rangle }{2}-((2b-1)\left(\frac{d{V}_A}{dr}+\frac{2V_A}{r}\right)\hfill \\ \hfill & +b\frac{V_A}{\rho }\frac{d\rho }{dr})\frac{\langle {z^{+}}^2\rangle -\langle {z^{-}}^2\rangle }{2}+\frac{V_A}{2}\left(\sqrt{\frac{\langle {z^{+}}^2\rangle }{\langle {z^{-}}^2\rangle }}\frac{d\langle {z^{-}}^2\rangle }{dr}-\sqrt{\frac{\langle {z^{-}}^2\rangle }{\langle {z^{+}}^2\rangle }}\frac{d\langle {z^{+}}^2\rangle }{dr}\right)\hfill \\ \hfill & =-E_D\left[\frac{{\langle {z^{+}}^2\rangle }^{1/2}}{{\lambda }^{-}}+\frac{{\langle {z^{-}}^2\rangle }^{1/2}}{{\lambda }^{+}}\right]+C_{sh}^{E_D}\frac{r_0|\Delta U|V_{A0}^2}{r^2};\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill & (U\mp V_A)\frac{d{L}^{\pm }}{dr}+\left(\frac{U}{r}\pm \frac{d{V}_A}{dr}\pm \frac{2V_A}{r}\right)L^{\pm }+2\left(a-\frac{1}{4}\right)\frac{U}{r}{L}_D\pm \left(b-\frac{1}{2}\right)\hfill \\ \hfill & \left(\frac{d{V}_A}{dr}+\frac{2V_A}{r}\right)L_D\pm b\frac{V_A}{2\rho }\frac{d\rho }{dr}{L}_D=0;\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill & U\frac{d{L}_D}{dr}+\frac{U}{r}{L}_D+V_A\left(\sqrt{\frac{L^{+}}{L^{-}}}\frac{d{L}^{-}}{dr}-\sqrt{\frac{L^{-}}{L^{+}}}\frac{d{L}^{+}}{dr}\right)+2\left(2a-\frac{1}{2}\right)\frac{U}{r}(L^{+}+L^{-})\hfill \\ \hfill & -(2b-1)\left(\frac{d{V}_A}{dr}+\frac{2V_A}{r}\right)(L^{+}-L^{-})-b\frac{V_A}{\rho }\frac{d\rho }{dr}(L^{+}-L^{-})=0,\hfill \end{array}$$

(13)

where we choose $a=1/2$ and $b=1/2$ to invoke two-dimensional turbulence [77,83,91]. In the transport Equations (10)–(11), the second term on the rhs is the turbulent shear source [78], and the third term on the rhs of Equation (10) is the pickup ion source of turbulence [76,78]. Here, we assume that the pickup ion source of turbulence provides equal energy to the forward and backward propagating modes. Note that the pickup ion source of turbulence is not included in the transport equation for the residual energy because the pickup ions generate Alfvénic fluctuations in the outer heliosphere, which leads to zero residual energy. The parameters $C_{sh}^{\pm }$ and $C_{sh}^{E_D}$ denote the strength of the stream-shear source of turbulence for the forward and backward propagating modes and the residual energy, $|\Delta U|$ is the difference between the fast and slow solar wind speeds, and $V_{A0}$ is the Alfvén velocity at a reference point $r_0$. We use $V_A^{\prime }=50$ kms${}^{-1}$. The parameter $f_D$ represents the fraction of the pickup ion source of turbulence that drives the turbulence in the outer heliosphere [9].

The turbulence transport Equations (10)–(13) include the magnetic field in the form of the Alfvén velocity $V_A$. The magnetic field is given by Weber and Davis [106]

$$B=B_a{\left(\frac{r_a}{r}\right)}^2{\left[1+{\left(\frac{{\omega }_{a}r}{U}\right)}^2{\left(1-{\left(\frac{ra}{r}\right)}^2\right)}^2{sin}^2\theta \right]}^{1/2},$$

(14)

where the subscript a represents the reference point $r_a$. We assume the reference point $r_a\approx 10$R${}_{\odot }$, where R${}_{\odot }$ is a solar radius. We use ${\omega }_a=2.9\times {10}^{-6}$ rads${}^{-1}$, $B_a=2.08\times {10}^3$ nT, and $\theta ={90}^{\circ }$.The conservation of mass, i.e, $\rho U{r}^2={\rho }_a{U}_a{r}_a^2$, yields the solar wind density,

$$\rho ={\rho }_a{\left(\frac{r_a}{r}\right)}^2.$$

With the inclusion of electrons, the 1D steady-state transport equations for the solar wind proton and electron pressures become

$$U\frac{d{P}_p}{dr}+2\gamma \frac{U}{r}{P}_p=(\gamma -1)\left[{\nu }_{pe}(P_e-P_p)+f_p{S}_t\right];$$

(15)

$$U\frac{d{P}_e}{dr}+2\gamma \frac{U}{r}{P}_e=(\gamma -1)\left[{\nu }_{ep}(P_p-P_e)-\nabla ·{\mathbf{q}}_e+(1-f_p)S_t\right],$$

(16)

where $S_t$ is the turbulent heating term, $P_p$ is the solar wind proton pressure, $P_e$ is the solar wind electron pressure, $f_p$ is the fraction of the turbulence energy that heats the solar wind protons, and ${\nu }_{pe}$ and ${\nu }_{ep}$ are the proton–electron and electron–proton Coulomb collision rates. The Coulomb collision frequencies are approximately balanced, so $n_e{\nu }_{ep}\approx n_p{\nu }_{pe}$. We also assume that the electron density is approximately equal to the proton density, i.e., $n_e\approx n_p$. The first term inside the squared brackets of (15) and (16) is the proton–electron Coulomb collision term. The Coulomb collision frequency ${\nu }_{pe}$ is given by Cranmer et al. [101],

$${\nu }_{pe}\approx 8.4\times {10}^{-9}\left(\frac{n_e}{2.5{\mathrm{cm}}^{-3}}\right){\left(\frac{T_e}{{10}^5\mathrm{K}}\right)}^{-3/2}{s}^{-1},$$

(17)

which is derived by assuming massive protons, i.e., $m_p\gg m_e$, and the rate of temperature equilibrium as described by Spitzer [107,108,109]. For electron-proton collisions, this collision frequency produces a large mean free path (mfp) of ∼500–1500 au at 1 au [101]. In the case of electron-electron collisions, the mean free path is of the order of ∼0.5 au at 1 au [107,110].

In Equation (16), the term ${\mathbf{q}}_e$ is the electron heat flux [111]. The proton heat flux is neglected because it is approximately zero for the Maxwellian core protons, and hence it is not important, unlike the electron heat flux. Equations (15) and (16) are a model for isotropic electron and proton pressures, meaning that we implicitly assume both protons and electrons are dominated by the Maxwellian core, i.e., the assumed 1D radially symmetric model-Maxwellian distribution means no heat flux at all, and that parallel and perpendicular contributions to the pressure/temperature are smaller and negligible [105]. We use the electron heat flux parallel to the magnetic field [101], which corresponds to the strahls.

The Spitzer and Härm [112] collisional model for the electron heat flux, $q_{\left|\right|}=-{\kappa }_{\left|\right|}{\nabla }_{\left|\right|}{T}_e\left(r\right)$ yields a larger electron temperature at 1 au than that observed [113]. Hollweg [113] introduced a collisionless form of the electron heat flux, which produced an electron temperature equivalent to the observed electron temperature, but nonlocal effects were ignored [114]. The Hollweg [113] model neglected electron-wave/turbulence scattering, which significantly affects electron transport [105]. Here, we use an empirical expression for the electron heat flux $q_{\left|\right|}$ [101] that was obtained by fitting the electron heat flux measured by Helios 2 from 0.3 to 1 au [115]. The fitted result for the electron heat flux is given by,

$$ln\left(\frac{q_{\left|\right|,e}}{q_0}\right)=-0.7037-2.115x-0.2545x^2,$$

(18)

where $x\equiv ln(r/1au)$ and $q_0=0.01$ erg cm${}^{-2}$ s${}^{-1}$. In a spherically symmetric coordinate system, the term $\nabla ·{\mathbf{q}}_e$ can be written as [101],

$$\nabla ·{\mathbf{q}}_e=\frac{1}{r^2}\frac{d}{dr}\left(r^2{q}_{\left|\right|}{cos}^2\varphi \right),$$

(19)

and $\varphi$ is the Parker spiral angle, and is given by

$$tan\varphi =\frac{\mathsf{\Omega }rsin\theta }{U},$$

(20)

and $\mathsf{\Omega }=2.7\times {10}^{-6}$ rad s${}^{-1}$ is the solar rotation frequency. The parameter $\theta$ is the colatitude angle, and we assume $\theta ={90}^{\circ }$.

In Equations (15) and (16), $S_t$ represents a turbulent heating term, which is distributed between electrons and protons. For example, it is thought that 60% of the turbulence energy heats solar wind protons [11,100], and 40% heats solar wind electrons [11]. We assume $f_p=0.6$ throughout the helioshere [11]. However, a few authors also assumed that $f_p=f_p\left(r\right)$ [102,116,117]. The turbulent heating term $S_t$ can be derived from a von Kármán–Taylor phenomenology, and is given by [23,34],

$$\begin{array}{cc}\hfill S_t& =\alpha m_p{n}_p\left[2\frac{\langle {z^{+}}^2\rangle {\langle {z^{-}}^2\rangle }^{1/2}}{{\lambda }^{+}}+2\frac{\langle {z^{-}}^2\rangle {\langle {z^{+}}^2\rangle }^{1/2}}{{\lambda }^{-}}+E_D\left(\frac{\langle {z^{+}}^2\rangle }{{\lambda }^{-}}+\frac{\langle {z^{-}}^2\rangle }{{\lambda }^{+}}\right)\right].\hfill \end{array}$$

(21)

Here, the first term inside the squared brackets is a nonlinear dissipation term corresponding to the energy in forward propagating modes, the second term corresponds to backward propagating modes, and the third term corresponds to the residual energy.

Results: Proton–Electron Heating in the Upwind and Downwind Directions

We solve the coupled Equations (10)–(16) using a Runge–Kutta fourth order method over the distance of the perihelion of the first orbit of the PSP (i.e., 0.17 au) to 75 au in the upwind and downwind directions. The boundary conditions are shown in

Table 1. Boundary values at 0.17 au for the upwind, and the downwind directions. The upwind boundary conditions are obtained from PSP measurements [19], and the downwind boundary conditions are chosen in such a way that the theoretical results closely reflect the measured values and observed trends of the $Pioneer$ 10 data. We assume that the electron density is approximately equal to the proton density, $n_e\approx n_p$.

Parameters	Upwind	Downwind
$\langle {z^{+}}^2\rangle$	${10}^4$ km${}^2$s${}^{-2}$	$7\times {10}^3$ km${}^2$s${}^{-2}$
$\langle {z^{-}}^2\rangle$	1190 km${}^2$s${}^{-2}$	1190 km${}^2$s${}^{-2}$
$E_D$	$-281.2$ km${}^2$s${}^{-2}$	−281.2 km${}^2$s${}^{-2}$
$L_{+}$	$8.93\times {10}^7$ km${}^3$s${}^{-2}$	$2.63\times {10}^9$ km${}^3$s${}^{-2}$
$L_{-}$	$6\times {10}^8$ km${}^3$s${}^{-2}$	$4.99\times {10}^8$ km${}^3$s${}^{-2}$
$L_D$	$-4.22\times {10}^7$ km${}^3$s${}^{-2}$	−$1.69\times {10}^8$ km${}^3$s${}^{-2}$
$P_e$	$8.28\times {10}^{-10}$ kgm${}^{-1}$s${}^{-2}$	$8.28\times {10}^{-10}$ kgm${}^{-1}$s${}^{-2}$
$P_p$	$1.04\times {10}^{-9}$ kgm${}^{-1}$s${}^{-2}$	$1.04\times {10}^{-9}$ kgm${}^{-1}$s${}^{-2}$

, and

Table 2. Model parameters. In the downwind direction, the values for $L=2$ au and $\zeta ={170}^{\circ }$ are obtained by fitting the hydrogen neutrals by $exp(-L\zeta /sin\zeta )$ along the trajectory of $Pioneer$ 10 [118].

Parameters	Upwind	Downwind
$C_{sh}^{+}$	0.1	0.3
$C_{sh}^{-}$	0.1	0.3
$C_{sh}^{E_D}$	0.01	0.9
U	400 km s${}^{-1}$	400 kms${}^{-1}$
$\Delta U$	200 km s${}^{-1}$	200 kms${}^{-1}$
$V_{A0}$	90 km s${}^{-1}$	90 kms${}^{-1}$
$n_0$	300 cm${}^{-3}$	300 cm${}^{-3}$
$\alpha$	0.1	0.1
$f_D$	0.1	0.02
$n_H^{\infty }$	0.1 cm${}^{-3}$	0.1 cm${}^{-3}$
$n_{nw}$	5 cm${}^{-3}$	5 cm${}^{-3}$
${\tau }_{ion}$	${10}^6$ s	${10}^6$ s
L	8 au	2 au
$\zeta$	$0^{\circ }$	${170}^{\circ }$
$r_0$	0.17 au	0.17 au

shows the parameter values used in the model for the upwind and downwind directions. The parameter values are chosen so that the theoretical results are consistent with observations. We find that other choices of the parameters yield results that do not fit the available observations very well. However, a detailed analysis of the possible parameter set and its statistical fit to the available data sets has not yet been undertaken. How we derive the observed turbulent quantities can be found in our series of papers [16,17,18,44,76,91].

Results of the solar wind proton (dashed curve) and electron (solid curve) temperatures in the upwind direction are shown in the left panel of Figure 1. With the inclusion of electrons, the solar wind plasma temperature is moderately affected by turbulence, electron heat flux and the Coulomb collisions between the solar wind electrons and protons. However, the effect of Coulomb collisions is negligible. The theoretical proton temperature decreases monotonically with distance until ∼20 au, and then increases with increasing heliocentric distance. The theoretical proton temperature is in reasonable agreement with the proton temperature measured by $PSP$ Adhikari et al. [19] (red “*” symbol), $Voyager$ 2 Adhikari et al. [18] (red “.” symbol), and $NHSWAP$ Zank et al. [23] (blue “.” symbol). The theoretical electron temperature decreases gradually until ∼1 au, and further decreases slowly between 1∼10 au. The slow decrease in the electron temperature is due to the electron heat flux. Breech et al. [11] found a shelf-like region between 1 and 10 au, which is not found in the top left panel of Figure 1. The theoretical electron temperature is similar to that measured by $PSP$ Moncuquet et al. [119] (green “.” symbol), and $Ulysses$ McComas et al. [120] (magenta “*” symbol). We obtain $Ulysses$ data from the website https://cdaweb.gsfc.nasa.gov/index.html/. The $Ulysses$ measurements correspond to a solar minimum period from 1994 to 1997. We only selected those data in which the position of $Ulysses$ is within a latitude of $\pm {30}^{\circ }$.

In the downwind direction (∼180${}^{\circ }$), $Pioneer$ 10 shows no significant or obvious increase in proton temperature (red “.” symbols in the top right panel of Figure 1), unlike the $Voyager$ 2 measurements in the upwind direction. As shown in Figure 1 in [118], the hydrogen distribution in a narrow conical region in the downwind direction is less than that in the upwind direction, leading to a reduced pickup ion number density and hence source of turbulence yielding a lower heating rate in the downwind direction compared to the upwind direction. The top right panel of Figure 1 shows that the theoretical proton temperature (dashed curve) decreases gradually with increasing heliocentric distance, similar to the proton temperature measured by $Pioneer$ 10. The electron temperature (solid curve) decreases more slowly than the proton temperature within ∼3–4 au. Electron and proton temperatures are almost flat beyond ∼50–60 au.

The bottom left and right panels of Figure 1 show the theoretical and observed proton and electron entropy as a function of heliocentric distance in the upwind and downwind directions, respectively. In the upwind direction, the theoretical and observed electron entropy increases as a function of heliocentric distance, which is about 12.48% higher at 75 au than at 0.17 au. Similarly, the theoretical and observed proton entropy increases monotonically with distance. The proton entropy increases by about 17.24% at 75 au from 0.17 au. In the downwind direction, the electron entropy increases by about 9.11% and the proton entropy by about 12.5% from 0.17 au.

The top left and right panels of Figure 2 show the energy in forward propagating modes as a function of heliocentric distance in the upwind and downwind directions, respectively. In the left panel, the theoretical energy in forward propagating modes is compared with $PSP$ and $Voyager$ 2 measurements. In the right panel, the theoretical energy in forward propagating modes is compared with $Pioneer$ 10 measurements. The theoretical and observed results are in good agreement. Since the $Pioneer$ 10 magnetometer data are not available after 10 au, the observed energy in forward propagating modes is plotted over the distance 1–10 au. In the upwind direction, the theoretical energy in forward propagating modes tends to flatten after ∼10 au, while in the downwind direction, the theoretical energy in forward propagating modes is slightly flat after ∼50 au. Similarly, the theoretical energy in backward propagating modes is in good agreement with $PSP$ and $Voyager$ 2 measurements in the upwind direction (Figure 2, the left panel of the second row) and $Pioneer$ 10 measurements in the downwind direction (Figure 2, the right panel of the second row). However, after ∼10 au and ∼30 au, the pickup ions begin to affect the theoretical energy in backward propagating modes in the upwind and downwind directions, respectively. The theoretical results of the energy in forward and backward propagating modes clearly show that the turbulence energy is different in the upwind and strictly downwind directions.

The theoretical fluctuating magnetic energy in the upwind direction (Figure 2, left panel of the third row) decreases as $r^{-2.33}$, while that in the downwind direction decreases as $r^{-2.56}$ (Figure 2, right panel of the third row). Thus, turbulent magnetic energy in the downwind direction decreases faster than that in the upwind direction. The results show that the theoretical magnetic energy is in good agreement with the measured values of $PSP$ (red “*” symbols) and $Voyager$ 2 (red “.” symbols) in the upwind direction, and the measured values of $Pioneer$ 10 (red “.” symbols) in the downwind direction. Similarly, in the upwind direction, the theoretical turbulent kinetic energy (bottom left panel of Figure 2) decreases gradually with increasing distance until ∼10 au, and then increases with distance. However, in the downwind direction, the theoretical turbulent kinetic energy decreases monotonically until ∼30–40 au, and then increases slightly with distance. The theoretical turbulent kinetic energy is in good agreement with the observed turbulent kinetic energy in the upwind and downwind directions.

The top left and right panels of Figure 3 plot the normalized residual energy as a function of heliocentric distance in the upwind and downwind directions. The theoretical normalized residual energy shows a similar radial profile to that of the observed normalized residual energy in both directions, although the scatter is large. The normalized residual energy in the downwind region does not tend towards equipartition between the fluctuating kinetic and magnetic energy with increasing distance, as in the upwind region. In the upwind direction, the fluctuations are increasingly magnetic dominated over the distance 2–10 au. Only after that does ${\sigma }_D\to 0$. This means that fluctuations in the upwind region beyond ∼10 au are more Alfvénic than that in the downwind region. The Alfvénicity of the solar wind in the outer heliosphere is due to wave excitation by the pickup ions. In the upwind direction (Figure 3, left panel of the second row), the theoretical normalized cross-helicity is high at the left boundary, decreases gradually as a function of heliocentric distance. In the downwind direction (Figure 3, right panel of the second row), the normalized cross-helicity is high at 0.17 au, and decreases gradually with increasing heliocentric distance. The theoretical normalized cross-helicity is in reasonable agreement with the observed normalized cross-helicity in the upwind direction, while the theoretical cross-helicity is greater than the observed cross-helicity in the downwind direction between 1 and 10 au. The bottom left and right panels of Figure 3 show the Alfvén ratio as a function of heliocentric distance in the upwind and downwind directions, respectively. The results show that the theoretical and observed Alfvén ratios are consistent.

The correlation length is an important quantity because it determines the decay rate of turbulence, and hence controls the turbulent heating rate. The top left and right panels of Figure 4 show the correlation length corresponding to forward and backward propagating modes, and the residual energy as a function of heliocentric distance in the upwind and downwind directions. The theoretical correlation length corresponding to forward and backward propagating modes (the red and blue curves on the left panel, respectively) increases gradually with distance until ∼10 au, and then remains approximately constant with increasing heliocentric distance. Since the pickup ions excite Alfvén waves in the outer heliosphere beyond the ionization cavity, this leads to an increase in turbulent energy, which results in a flattening or a decrease of the correlation length through the adjustment of the Kolmogorov dissipation rate. The correlation length corresponding to forward and backward propagating modes (the red and blue curves on the left panel) in the downwind direction increases gradually until ∼40 au, and then remains approximately constant until ∼75 au. The correlation length of the residual energy (the green curves on the left and right panels) increases gradually as a function of heliocentric distance in the upwind and downwind directions.

The middle left and right panels of Figure 4 show the theoretical and observed correlation length of the magnetic field fluctuations ${\lambda }_b$ in the upwind and downwind directions, respectively. The theoretical correlation length ${\lambda }_b$ is in reasonable agreement with the observed ${\lambda }_b$. The bottom left and right panels show the comparison between the theoretical and observed correlation length of the velocity fluctuations ${\lambda }_u$ as a function of heliocentric distance. In the upwind direction (left panel), the theoretical correlation length ${\lambda }_u$ is in good agreement with the observed ${\lambda }_u$ of PSP measurements, but does not agree with $Voyager$ 2 measurements. Similarly, in the downwind direction, the theoretical ${\lambda }_u$ does not agree with $Pioneer$ 10 measurements.

5. Turbulence Cascade in the Inertial Range throughout the Heliosphere

The turbulent cascade rate is directly related to the heating of the solar wind. Several authors have studied the turbulence cascade via observations [121,122,123,124,125,126,127], global MHD simulation [93], turbulence transport theory [34,103,128], and remote sensing observations [129,130,131]. These studies find that the turbulent heating rate is high near the Sun, and decreases with increasing heliocentric distance.

Vasquez et al. [121] used ACE measurements to examine the turbulence cascade as a function of proton temperature and solar wind speed. They argued that the cascade rate cannot be estimated with a Kolmogorov form, and better agreement is found when using an Irosnikov–Kraichnan form. However, Bandyopadhyay et al. [126] using PSP measurements showed that the cascade rate calculated from the Politano–Pouquet third-order law [132,133] and a von-Kármán phenomenology, which follows a Kolmogorov form, are approximately similar. Smith et al [123] studied the cascade rate of the energy in forward propagating modes, the total turbulent energy, and the cross-helicity from the third-order moments of the velocity and magnetic field fluctuations, finding that the cascade rate strongly depends on the correlation between the velocity and magnetic field fluctuations. Recently, Pine et al. [127] used the expression from Vasquez et al. [121], the turbulence transport theory, and wave excitation by newborn interstellar ions to study the turbulence cascade rate using $Voyager$ 1 and 2 magnetometer data. They found that the first two expressions are very consistent, and the latter becomes important beyond 10 au. Their results also support idea that the heating of the solar wind plasma for $R<10$ au is mainly accomplished by turbulence of a solar origin, while the heating for $R>10$ au seems to be almost entirely due to the wave excitation by interstellar PUIs.

Using density modulation indices ${ϵ}_{N_e}$ obtained using angular broadening observations of the Crab Nebula from 1952 to 2013, Sasikumar Raja et al. [131] calculated the solar wind proton heating rate using ${ϵ}_{k_i}\left(r\right)=c_0{\rho }_p{k}_i\left(r\right)\delta v_{l{k}_i}^3\left(r\right)$ erg cm${}^{-3}$s${}^{-1}$ [129,134]. They found that the heating rate varies from $1.01\times {10}^{-8}$ to $\sim 1.58\times {10}^{-14}$ erg cm${}^{-3}$ s${}^{-1}$ over a distance $5-45$ R${}_{\odot }$. Similarly, Ingale [135] also found that the heating rate varies between $3\times {10}^{-8}-{10}^{-15}$ erg cm${}^{-3}$ s${}^{-1}$ over a heliocentric distance $2-174$ R${}_{\odot }$. Using the coupled quasi-2D and NI/slab turbulence transport model equations [78], Adhikari et al. [128] found that the heating rate associated with quasi-2D turbulence is dominant. The heating rate corresponding to quasi-2D turbulence is $1.06\times {10}^{-4}-1.73\times {10}^{-14}$ erg cm${}^{-3}$ s${}^{-1}$ between $1.6-100$ R${}_{\odot }$, and that of NI/slab turbulence is $4.24\times {10}^{-7}-1.11\times {10}^{-14}$ erg cm${}^{-3}$ s${}^{-1}$ between $1.3-100$ R${}_{\odot }$.

Here, we derive the turbulence cascade rate for the energy in forward and backward propagating modes, the total turbulent energy, the cross-helicity, the residual energy, the fluctuating kinetic energy, and the fluctuating magnetic energy from two perspectives: (i) turbulence transport theory Zank et al. [77] (hereafter, called model 1), and (ii) a dimensional analysis between the power spectrum in the energy-containing range and the inertial range Adhikari et al. [136] (hereafter, called model 2). We first discuss the turbulence cascade rate via turbulence transport theory. The turbulence transport Equations (5)–(9) or (10)–(13) are derived without considering the viscous term. However, in deriving the turbulence transport model, the implicit assumption is that the inertial range of the turbulence spectrum is self-similar to the energy input rate in the inertial range, thus balancing the dissipation rate. This indicates that the energy transfer rate in the turbulence model is determined by the nonlinear term, which determines the heating of the solar wind plasma. Although turbulence transport modeling is a theory describing the transport of energy-containing range fluctuations, the decay of these fluctuations follows a Kolmogorov or Iroshnikov–Kraichnan (I–K) phenomenology. Therefore, the turbulence transport model calculates the balanced rate at which energy enters and exits the inertial range. It is found observationally that magnetic field fluctuations exhibit either a Kolmogorov or I–K type of power law in the inertial range [60,61,137,138,139,140,141,142].

The nonlinear term acts as a dissipative source for coronal and solar wind heating. The first term on the rhs of Equations (10) and (11) determines the cascade of the energy in forward and backward propagating modes, and the residual energy. These terms can be expressed as

$$\begin{array}{cc}\hfill & {ϵ}_{\langle {z^{\pm }}^2\rangle }^1=2\alpha \frac{\langle {z^{+}}^2\rangle {\langle {z^{-}}^2\rangle }^{1/2}}{{\lambda }^{\pm }};\hfill \\ \hfill & {ϵ}_{E_D}^1=\alpha E_D\left(\frac{{\langle {z^{-}}^2\rangle }^{1/2}}{{\lambda }^{+}}+\frac{{\langle {z^{+}}^2\rangle }^{1/2}}{{\lambda }^{-}}\right),\hfill \end{array}$$

(22)

where the ${ϵ}_{\langle {z^{\pm }}^2\rangle }$ denote the turbulence cascade rate of energy in forward/backward propagating modes, and ${ϵ}_{E_D}$ the turbulence cascade rate of the residual energy. The parameter $\alpha$ is a von Kármán–Taylor constant, and we choose $\alpha =0.03$[126]. Since the rhs term is positive in the ${ϵ}_{\langle {z^{\pm }}^2\rangle }$ Equation, the energy in forward and backward propagating modes always shows a forward cascade. However, for the turbulent cascade of the residual energy, there occurs a forward cascade for $E_D>0$, and a reverse cascade for $E_D<0$. Similarly, the turbulence cascade corresponding to the turbulent kinetic energy ${ϵ}_{\langle u^2\rangle }$, the turbulent magnetic energy ${ϵ}_{\langle E_b\rangle }$, the total turbulent energy ${ϵ}_{E_T}$, the (normalized) cross-helicity (${ϵ}_{{\sigma }_c}$) ${ϵ}_{E_C}$, and the normalized residual energy ${ϵ}_{{\sigma }_D}$ can be written as

$$\begin{array}{cc}\hfill & {ϵ}_{\langle u^2\rangle }^1=\frac{{ϵ}_{\langle {z^{+}}^2\rangle }^1+{ϵ}_{\langle {z^{-}}^2\rangle }^1+2{ϵ}_{E_D}^1}{4};{ϵ}_{\langle E_b\rangle }^1=\frac{{ϵ}_{\langle {z^{+}}^2\rangle }^1+{ϵ}_{\langle {z^{-}}^2\rangle }^1-2{ϵ}_{E_D}^1}{4};\hfill \\ \hfill & {ϵ}_{E_T}^1=\frac{{ϵ}_{\langle {z^{+}}^2\rangle }^1+{ϵ}_{\langle {z^{-}}^2\rangle }^1}{2};{ϵ}_{E_C}^1=\frac{{ϵ}_{\langle {z^{+}}^2\rangle }^1-{ϵ}_{\langle {z^{-}}^2\rangle }^1}{2};{ϵ}_{{\sigma }_c}^1=\frac{{ϵ}_{E_C}^1}{{ϵ}_{E_T}^1},{ϵ}_{{\sigma }_D}^1=\frac{{ϵ}_{E_D}^1}{{ϵ}_{E_T}^1},\hfill \end{array}$$

(23)

where the superscript “1” refers to model 1 derived from turbulence transport theory. Similar to the residual energy, the cross-helicity also shows a forward cascade when ${ϵ}_{\langle {z^{+}}^2\rangle }>{ϵ}_{\langle {z^{-}}^2\rangle }$, and a reverse cascade when ${ϵ}_{\langle {z^{+}}^2\rangle }<{ϵ}_{\langle {z^{-}}^2\rangle }$.

The second approach is based on the use of a dimensional analysis [136] to derive the turbulence cascade rate. Figure 5 is a schematic diagram, illustrating the power spectral density (PSD) of Elsässer energies $\langle z^{\pm 2}\rangle$, turbulent magnetic energy density $E_b$, and the turbulent kinetic energy $\langle u^2\rangle$ as a function of wavenumber k, where $P\left(k\right)=A{k}^{-1}$ defines the power spectrum in the energy-containing range, and $P\left(k\right)=C_K{ϵ}^{2/3}{k}^{-5/3}$ the power spectrum in the inertial range. The parameter A needs to be determined, $C_K=1.6$ is a Kolmogorov constant, and $ϵ$ is the turbulence cascade rate. The wavenumber $k_b$ separates the energy-containing range and the inertial range. Since the PSD in the energy-containing range is equal to the PSD in the inertial range at $k_b$, we can write

$$\begin{array}{cc}\hfill & A{k}^{-1}{{|}_{k_b}=C_K{ϵ}^{2/3}{k}^{-5/3}|}_{k_b};\hfill \\ \hfill & \Rightarrow ϵ={\left[\frac{A}{C_K}{k}_b^{2/3}\right]}^{3/2}.\hfill \end{array}$$

(24)

Integrating the PSD $P\left(k\right)=A{k}^{-1}$ from $k_{inj}$ to $k_b$ yields

$$\begin{array}{cc}\hfill & E={\int }_{k_{inj}}^{k_b}A{k}^{-1}dk,\hfill \\ \hfill & \Rightarrow A=\frac{E}{log\left(\frac{k_b}{k_{inj}}\right)},\hfill \end{array}$$

(25)

where we assume that $k_{inj}\sim 1.07\times {10}^{-9}$ km${}^{-1}$, which is equivalent to a solar rotation of ∼27 days [136]. The denominator on the rhs of Equation (25) is dimensionless, and the numerator has the dimension of energy, implying that the parameter A has the dimension of energy. From Equations (24) and (25), we obtain the turbulence cascade rate as,

$$ϵ={\left[\frac{E}{C_K}\frac{k_b^{2/3}}{log\left(\frac{k_b}{k_{inj}}\right)}\right]}^{3/2}=\frac{E^{3/2}}{{\left[C_{K}log\left(\frac{1}{k_{inj}\lambda }\right)\right]}^{3/2}\lambda }.$$

(26)

Equation (26) calculates the turbulent cascade rate through the inertial range, which is derived by assuming that the solar wind fluctuations exhibit a Kolmogorov-like power law, and $k_b\equiv {\lambda }^{-1}$, where $\lambda$ is the correlation length corresponding to the turbulent energy E.

We use Equation (26) to calculate the turbulence cascade of the fluctuating magnetic energy density $E_b$, the energy in forward/backward propagating modes $\langle {z^{\pm }}^2\rangle$, and the turbulent kinetic energy $\langle u^2\rangle$. The turbulence cascade rates are given by

$$\begin{array}{cc}\hfill & {ϵ}_{E_b}^2=\frac{E_b^{3/2}}{{\left[C_{K}log\left(\frac{1}{k_{inj}{\lambda }_b}\right)\right]}^{3/2}{\lambda }_b};{ϵ}_{\langle {z^{\pm }}^2\rangle }^2=\frac{{\langle {z^{\pm }}^2\rangle }^{3/2}}{{\left[C_{K}log\left(\frac{1}{k_{inj}{\lambda }^{\pm }}\right)\right]}^{3/2}{\lambda }^{\pm }};\hfill \\ & {ϵ}_{\langle u^2\rangle }^2=\frac{{\langle u^2\rangle }^{3/2}}{{\left[C_{K}log\left(\frac{1}{k_{inj}{\lambda }_u}\right)\right]}^{3/2}{\lambda }_u}.\hfill \end{array}$$

(27)

The turbulence cascade rate of the total turbulent energy, the (normalized) cross-helicity are given by,

$$\begin{array}{cc}\hfill & {ϵ}_{E_T}^2=\frac{{ϵ}_{\langle {z^{+}}^2\rangle }^2+{ϵ}_{\langle {z^{-}}^2\rangle }^2}{2};{ϵ}_{E_C}^2=\frac{{ϵ}_{\langle {z^{+}}^2\rangle }^2-{ϵ}_{\langle {z^{-}}^2\rangle }^2}{2};\hfill \\ \hfill & {ϵ}_{E_D}^2={ϵ}_{\langle u^2\rangle }^2-{ϵ}_{E_b}^2;{ϵ}_{{\sigma }_c}^2=\frac{{ϵ}_{E_C}^2}{{ϵ}_{E_T}^2},\hfill \end{array}$$

(28)

where the superscript “2” refers to model 2 derived from the dimensional analysis.

Results: Turbulence Cascade in the Upwind and Downwind Directions

In this section, we discuss the upwind and downwind turbulence cascade rates from the perihelion of the first orbit of $PSP$ to 75 au. Figure 6 shows the turbulence cascade rate as a function of heliocentric distance in the upwind and downwind directions. In the figure, the solid line denotes the turbulence cascade rate calculated from model 1, i.e., based on turbulence transport theory, and the dashed line the turbulence cascade rate obtained from model 2, i.e., based on dimensional analysis. The red and blue curves represent the turbulence cascade rate in the downwind and upwind directions, respectively. Panels (a–d) of Figure 6 show the turbulence cascade rate of energy in forward propagating modes, energy in backward propagating modes, turbulent magnetic energy density, and the turbulent kinetic energy, respectively. The results show that the turbulence cascade rate in the downwind direction is higher than that in the upwind direction between ∼0.3 and ∼10–20 au. After ∼10–20 au, the turbulence cascade rate in the upwind direction is larger than that in the downwind direction because more turbulence is generated by pickup ions in the upwind direction than in the downwind direction. In the upwind direction, the turbulence cascade rate corresponding to energy in forward/backward propagating modes $\langle {z^{\pm }}^2\rangle$, fluctuating magnetic energy density $E_b$, and fluctuating kinetic energy $\langle u^2\rangle$ decreases with increasing heliocentric distance until $\sim 10$ au, and then increases or flatten with distance. In the downwind direction, the turbulence cascade rates of $\langle {z^{\pm }}^2\rangle$ and $E_b$ decrease as distance increases until ∼40 au, and then flattens or decreases slightly with increasing distance. Models 1 and 2 reproduce similar turbulence cascade rates of $\langle {z^{\pm }}^2\rangle$, and $E_b$ in the upwind and downwind directions. For the turbulence cascade rate of $\langle u^2\rangle$, model 2 produces a larger rate than that of model 1. This is because the turbulence cascade rate in model 2 is inversely proportional to the correlation length, and the correlation length of velocity fluctuations is smaller in the outer heliosphere, yielding a higher turbulence cascade rate ${ϵ}_{\langle u^2\rangle }$.

Panel (e) of Figure 6 shows the turbulence cascade rate of the normalized residual energy ${ϵ}_{{\sigma }_D}$ as a function of heliocentric distance. The ${ϵ}_{{\sigma }_D}$ in the downwind direction (solid and dashed red curves) is larger than that in the upwind direction (solid and dashed blue curves) from 0.17 au to 20 au, and then it reverses. Model 1 produces a negative ${ϵ}_{{\sigma }_D}$ throughout the heliosphere, indicating that the turbulence cascade rate of the residual energy exhibits a reverse cascade within the inertial range. We note that the dimensional analysis is predicated on a particular form of the inertial range which is unlikely to be satisfied by a quantity that cannot be described by I-K theory. There is no theory for the residual energy analogous to Kolmogorov and so there is no reason to assume that the energy into the inertial range is balanced by the dissipation rate. Therefore, we only show the result ${ϵ}_{{\sigma }_D}$ of model 1.

Panel (f) of Figure 6 displays the turbulence cascade rate of normalized cross-helicity ${ϵ}_{{\sigma }_c}$ as a function of heliocentric distance. In the outer heliosphere, the ${ϵ}_{{\sigma }_c}$ in the upwind region is smaller than that in the downwind region because ${\sigma }_c$ is approximately zero beyond ∼20 au in the upwind region, but not in the downwind region (see the middle right panel of Figure 3).

6. Cosmic Ray Diffusion Tensor throughout the Heliosphere

The propagation of cosmic ray (CR) is affected by low-frequency turbulence throughout the heliosphere [41,42,43,44,45,143,144,145,146]. The power in the magnetic fluctuations and the magnetic correlation length are responsible for the scattering of the CR [147].

A detailed understanding of the evolution of magnetic turbulence is required to calculate the parallel, perpendicular and drift components of the CR diffusion tensor. To calculate the perpendicular and drift components, the quantity $\mathsf{\Omega }\tau$ needs to be evaluated, where $\mathsf{\Omega }$ is the particle gyrofrequency and $\tau$ an effective scattering time. It is generally believed that the magnetic field line wandering is the main physical mechanism that controls $\mathsf{\Omega }\tau$. There are two ways to model the magnetic field line wandering: quasi-linear theory (QLT) [147,148,149,150,151], and a nonperturbative approach [152,153,154,155,156].

QLT is a widely accepted theory to describe particle transport in a medium typically comprised of magnetized waves and/or fluctuations. In the case of purely 2D magnetic turbulence, QLT produces an infinitely large parallel diffusion ${\mathsf{\Lambda }}_{\left|\right|}=0$ since the pitch angle Fokker–Planck coefficient $D_{\mu \mu }^{2D}\left(\mu \right)=0$ [150,151]. The basic principle of QLT is to assume that the charged particle gyro-orbits are weakly perturbed by electromagnetic fluctuations. One can evaluate the diffusion coefficient by expanding the terms of the Vlasov or collisionless Boltzmann equation into mean and fluctuating fields [157]. Matthaeus et al. [158] proposed a nonlinear guiding center (NLGC) theory for the perpendicular diffusion of charged particles, including the effect of parallel scattering and dynamical turbulence, finding good agreement between theory and simulation. Inspired by the success of NLGC theory, Shalchi et al. [159] proposed a weakly nonlinear theory (WNLT) for the parallel and perpendicular diffusion of CRs. In WNLT, the nonlinear effect produces a resonance broadening, which allows charged particles to scatter through ${90}^{\circ }$ [159]. The pitch-angle Fokker–Planck coefficient is no longer zero at ${90}^{\circ }$ and, therefore, yields a reasonable parallel diffusion length [160,161]. However, due to the complexity of the extended nonlinear theories for an ab initio turbulence model, we use a simple QLT description to describe the parallel diffusion of energetic charged particles [41].

For the perpendicular diffusion of energetic charged particles, we use the NLGC theory [158]. If turbulence is assumed to be entirely slab [162,163] or if the Kubo number [164] of the turbulence is small [165], the original NLGC theory produces inaccurate results. There are more sophisticated extensions of NLGC theory [166,167,168]. Shalchi [166] introduced an important extensions of NLGC theory, namely the unified nonlinear theory (UNLT), which addresses the weaknesses of NLGC theory. Turbulence responsible for the scattering of charged particles can be described as the superposition of a dominant 2D component and a minority slab component [41,43,44,154,155,156], in which case UNLT theory and NLGC theory give very similar results. If the turbulence models do not contain a dominant 2D component, the UNLT theory should be used instead of NLGC theory. In this work, we use a 3D incompressible MHD transport model to describe $\langle b^2\rangle$ but to mimic the anisotropic character of solar wind turbulence, we decompose the magnetic energy density into a quasi-2D and a slab component in the ratio 80:20, i.e., 2D turbulence energy $\langle b_{2D}^2\rangle$$(=0.8\langle b^2\rangle$), and slab turbulence energy $\langle b_{slab}^2\rangle$$(=0.2\langle b^2\rangle$). Theory and observations predict that the ratio between the 2D and slab turbulence energy is 80:20 [62,103,169,170]. Similarly, we assume $2{\lambda }_{2D}^b={\lambda }_{slab}^b={\lambda }^b$ to calculate the 2D (${\lambda }_{2D}$) and slab (${\lambda }_{slab}$) magnetic correlation length, where ${\lambda }^b$ is the magnetic correlation length. Zhao et al. [43,44] have used a NI MHD turbulence model Zank et al. [78] to investigate the CR diffusion tensor throughout the heliosphere. The WKB model for Alfvén waves was used originally to study the radial and latitudinal dependence of the CR diffusion tensor [171,172].

In CR modulation studies, the drift effect caused by the turbulent magnetic field is often ignored, and a drift coefficient derived under the assumption of weak scattering is adopted [173,174,175]. However, theoretical and simulation results show that the drift coefficient decreases in the presence of turbulence [176,177]. The famous classical scattering limit uses the parallel mfp and gyro-radius to calculate the reduction factor in the drift coefficient [149,178]. Zhao et al. [43] derived an expression for the drift velocity in a strong turbulence, in which the reduction factor is similar to that proposed by Tautz et al. [163], and the first order approach proposed by Engelbrecht et al. [175,179].

Zhao et al. [44] showed observationally that the CR parallel, perpendicular, and radial mean free paths are influenced by solar cycle at 1 au, and throughout the heliopshere [180]. Moloto et al. [181] developed a 3D time-dependent ab initio CR modulation model by using time-dependent large-scale quantities such as the magnetic field, tilt angle, solar wind speed, turbulent magnetic energy and the correlation length of the magnetic field fluctuations, finding reasonable agreement with spacecraft measurements, and reproduced the major salient features of the observed CR intensity temporal profile [175,182,183].

6.1. Cosmic Ray Diffusion Tensor

The CR diffusion tensor

$${\kappa }_{ij}={\kappa }_{\perp }{\delta }_{ij}+\frac{B_i{B}_j}{B^2}({\kappa }_{\left|\right|}-{\kappa }_{\perp })+{ϵ}_{ijk}{\kappa }_A\frac{B_k}{B},$$

(29)

describes the scattering of CRs by fluctuations in the heliospheric magnetic field. The parameters $B_i$, $B_j$, $B_k$ are the magnetic field components, ${\delta }_{ij}$ is the Kronecker delta tensor, and ${ϵ}_{ijk}$ is the alternating or Levi–Civita tensor. Here, ${\kappa }_{\left|\right|}$ and ${\kappa }_{\perp }$ are the diagonal components of the CR diffusion tensor and describe diffusion along and perpendicular to the magnetic field, respectively. The diffusion term ${\kappa }_A$ is an off-diagonal antisymmetric component that describes the particle drift caused by large-scale curvature and gradients in the IMF. The magnetic field is given by Equation (24) [106].

For solar modulation models, the radial diffusion coefficient is given by

$${\kappa }_{rr}\equiv {\kappa }_{\left|\right|}{cos}^2\psi +{\kappa }_{\perp }{sin}^2\psi ,$$

(30)

where $\psi$ is the winding angle between the mean magnetic field and the radial direction,

$$tan\psi =\frac{\mathsf{\Omega }sin\theta }{U}(r-r_a).$$

(31)

Here, $\theta$ is the colatitude with respect to the solar axis rotation. We discuss the diffusion tensor in terms of the diffusion length scales ${\mathsf{\Lambda }}_{\left|\right|,\perp }=3k_{\left|\right|,\perp }/v$, where v is the particle speed [41,43,44,184].

6.2. CR Parallel Mean Free Path

Bieber et al. [185] calculated the mfp considering the superposition of the 2D component and a minority slab component, and found that the mfp of electrons and protons is consistent with the Palmer consensus [186]. However, the restriction of solar wind turbulence to a purely slab model [147], yielded CR mfps that exceeded the Palmer consensus [185]. This suggests that the superposition of a majority 2D component and a minority slab component needs to be considered to calculate the mfp. The geometric structure of turbulence comprising a superposition of a dominant 2D and a minority slab component turbulence proposed by Bieber et al. [185] is consistent with the theoretical prediction of Zank et al. [169,187] and the observational results of Bieber et al. [170]. Theoretical and observational studies show that about 80% of the turbulence energy resides in the dominant 2D turbulence, and about 20% of the turbulence energy resides in the minority slab turbulence [169,170,187]. The parallel mfp based on the standard QLT and assuming magnetostatic turbulence is approximately [41],

$${\mathsf{\Lambda }}_{\left|\right|}=3.1371\frac{B^{5/3}}{\langle {b^2}_{x,slab}\rangle }{\left(\frac{P}{c}\right)}^{1/3}{\lambda }_{slab}^{2/3}\left[1+\frac{7A/9}{(q+1/3)(q+7/3)}\right],$$

(32)

where

$$A={(1+s^2)}^{5/6}-1;q=\frac{5s^2/3}{1+s^2-{(1+s^2)}^{1/6}},\mathrm{and}s=0.746834R_L/{\lambda }_{slab}.$$

The parameter $P(\equiv pc/Ze$, where p is the momentum, c is the speed of light, $Ze$ is the particle charge) is the particle rigidity, B is the magnetic field strength, $R_L=P/Bc$ is the particle Larmor radius, and $\langle {b^2}_{x,slab}\rangle$ is the variance of the x component of the slab fluctuations, and is assumed to be $\langle b_{slab}^2\rangle /2$. Equation (32) is very close to the exact Fokker–Planck form. However, we note that expression (32) may not be fully accurate for very small rigidities when dynamical MHD turbulence is important [41,185]. The fractional term inside the squared brackets is particularly important in the outer heliosphere when the Larmor radius can be equal to or larger than the slab correlation length ${\lambda }_{slab}$. In this case, the ions no longer resonate with the turbulent MHD fluctuations in the inertial range, but resonate with fluctuations in the energy-containing range. As a result, the scaling of ${\mathsf{\Lambda }}_{\left|\right|}$ relative to the rigidity P and the correlation length changes from the inner to outer heliosphere according to the variation of the correlation length with heliocentric distance.

6.3. CR Perpendicular Mean Free Path

There are different methods to calculate the perpendicular mfp ${\mathsf{\Lambda }}_{\perp }$. One method is to use the hard sphere scattering approach based on a relaxation time approximation [188,189,190]. The other method was developed by Bieber and Matthaeus [191] using the Green–Kubo–Taylor formula. This is different from the QLT, and underestimates the diffusion at low energy. We use the NLGC theory for the perpendicular diffusion, which was first proposed by Matthaeus et al. [158]. The NLGC theory assumes that the perpendicular diffusion is controlled by the velocity of gyrocenters along the field line, which allows both the random walk of the field line and the parallel scattering of CRs along the magnetic field line to be incorporated simultaneously. We use the expression for the perpendicular mfp from the Zank et al. [192] analysis [150],

$$\begin{array}{cc}\hfill & \mathrm{for}{\lambda }_{slab}>{\mathsf{\Lambda }}_{\left|\right|}/\sqrt{3},\hfill \\ \hfill & {\mathsf{\Lambda }}_{\perp }={\left[\sqrt{3}\pi a^{2}C\frac{\langle b_{2D}^2\rangle }{B^2}{\lambda }_{2D}\right]}^{2/3}{\mathsf{\Lambda }}_{\left|\right|}^{1/3}{\left[1+\frac{4.33}{\sqrt{3}}\frac{{\left(a^{2}C\right)}^{1/3}}{{\left(\sqrt{3}\pi \right)}^{2/3}}\frac{\langle b_{slab}^2\rangle {\mathsf{\Lambda }}_{\left|\right|}^{2/3}}{{\langle b_{2D}^2\rangle }^{2/3}{\lambda }_{2D}^{2/3}{\left(B^2\right)}^{1/3}}\right]}^{2/3},\hfill \\ \hfill & \mathrm{and}\mathrm{for}{\lambda }_{slab}<{\mathsf{\Lambda }}_{\left|\right|}/\sqrt{3}\hfill \\ \hfill & {\mathsf{\Lambda }}_{\perp }={\left[\sqrt{3}\pi a^{2}C\frac{\langle b_{2D}^2\rangle }{B^2}{\lambda }_{2D}\right]}^{2/3}{\mathsf{\Lambda }}_{\left|\right|}^{1/3}{\left[1+3.091\frac{{\left(a^{2}C\right)}^{1/3}}{{\left(\sqrt{3}\pi \right)}^{2/3}}\frac{\langle b_{slab}^2\rangle }{{\langle b_{2D}^2\rangle }^{2/3}{\left(B^2\right)}^{1/3}}\frac{{\lambda }_{slab}}{{\mathsf{\Lambda }}_{\left|\right|}^{1/3}{\lambda }_{2D}^{2/3}}\right]}^{2/3},\hfill \end{array}$$

(33)

where $a^2$ is the coefficient corresponding to the gyrocenter velocity, which is numerically found to be ∼1/3, and $C=\mathsf{\Gamma }\left(\nu \right)/\left(2\sqrt{\pi }\mathsf{\Gamma }(\nu -1/2)\right)$ is a constant related to the spectral index $2\nu$, where $\mathsf{\Gamma }\left(\nu \right)$ is the gamma function.

The 2D turbulence power spectrum behavior at the lowest wavenumber, i.e., the outer scale, has an important impact on the perpendicular diffusion coefficients of charged particles in the heliosphere [146,193,194], including galactic cosmic ray proton and electron intensities. Engelbrecht et al. [193] calculated the galactic cosmic ray electron intensities using the observed size of magnetic islands [195,196], and different outer scales [136,145], finding that galactic cosmic ray electron differential intensities above kinetic energy 0.1 GeV are closer to the observed electron differential intensities, and lower than that below 0.1 GeV. Similarly, Engelbrecht et al. [194] proposed a new method to calculate the perpendicular diffusion coefficients by using the random ballistic decorrelation interpretation of the NLGC theory proposed by Ruffolo et al. [197]. The results show that the strength and anisotropy of solar energetic particle (SEP) are very sensitive to the pitch-angle dependence of the perpendicular diffusion coefficient.

6.4. Results: CR mfp in the Upwind and Downwind Directions

Figure 7 shows the radial dependence of the parallel mfp ${\mathsf{\Lambda }}_{\left|\right|}$ (blue lines), perpendicular mfp ${\mathsf{\Lambda }}_{\perp }$ (red lines), radial mfp ${\mathsf{\Lambda }}_{rr}$ (cyan lines), and ${\mathsf{\Lambda }}_{\perp }/{\mathsf{\Lambda }}_{\left|\right|}$ (green lines) for a proton with rigidity 445 MV (equivalent to 100 MeV proton) in the upwind direction (the solid line) and the downwind direction (the dashed line). In the upwind and downwind directions, the parallel mfp increases as a function of heliocentric distance. From 0.17 au to ∼0.3 au, ${\mathsf{\Lambda }}_{\left|\right|}$ in the upwind direction is smaller than that in the downwind direction, and then reverses from ∼0.3 au to ∼10 au, after which ${\mathsf{\Lambda }}_{\left|\right|}$ in the downwind direction becomes larger than that in the upwind region. At 75 au, the parallel mfp in the downwind direction is about 2 times greater than that in the upwind direction.

In the upwind and downwind directions, the perpendicular mfp ${\mathsf{\Lambda }}_{\perp }$ increases gradually as a function of heliocentric distance. The perpendicular mfp in the downwind direction is about five times larger than that in the upwind direction throughout the heliosphere. The radial mfp in the upwind direction is less than that in the downwind direction between 0.17 au and 0.3 au, and then reverses between 0.3 au and 10 au, after which ${\mathsf{\Lambda }}_{rr}$ in the downwind region is about three to five times greater than that in the upwind direction.

The results show that the ratio of the perpendicular to parallel mfp ${\mathsf{\Lambda }}_{\perp }/{\mathsf{\Lambda }}_{\left|\right|}$ in the upwind direction increases gradually from 0.17 au to $\sim 1$ au, and then decreases monotonically with increasing heliocentric distance. Similarly, in the upwind direction, the ratio ${\mathsf{\Lambda }}_{\perp }/{\mathsf{\Lambda }}_{\left|\right|}$ increases slowly from 0.17 au to 1 au, and then decreases gradually as a function of heliocentric distance. The result shows that the parallel mfp dominates the perpendicular mfp in the upwind and downwind directions throughout the heliosphere. The radial trends for the parallel, perpendicular, and radial mfps are similar with those of Zhao et al. [43,44].

7. Discussion and Conclusions

We developed a 1D steady-state turbulence model, including the separate dissipation of turbulence energy in proton and electron heating, that incorporated the electron heat flux and Coulomb collisions between solar wind protons and electrons. We solved the coupled turbulence transport model equations from 0.17 to 75 au in the upwind and downwind directions. In the upwind direction, we compared the theoretical results with $PSP$, $Voyager$ 2, and $NH$$SWAP$ measurements, and with $Pioneer$ 10 measurements in the downwind direction. Since the $Pioneer$ 10 magnetometer data are not available beyond 10 au, the observed quantities are plotted from 1 to 10 au in the downwind direction. We found reasonable agreement between the theoretical and observed results.

We derived expressions for the turbulence cascade rates for the energy in forward and backward propagating modes, total turbulent energy, normalized residual energy, normalized cross-helicity, fluctuating kinetic and magnetic energy from two perspectives: (i) a turbulence transport theory Zank et al. [77] (model 1), and (ii) a dimensional analysis between the power spectrum in the energy containing range and the inertial range Adhikari et al. [136] (model 2). Models 1 and 2 yield similar turbulence cascade rates for the energy in forward and backward propagating modes, magnetic energy density, and the normalized cross-helicity in the upwind and downwind directions. Finally, we calculated the CR diffusion tensor in the upwind and downwind directions. Our basic conclusions can be summarized using three categories.

A. Electron and proton heating in the upwind and downwind directions conclusions:

• Turbulence heats solar wind electrons and protons differently in the upwind and downwind directions. The proton and electron temperatures increase beyond ∼20 au in the upwind direction, but not in the downwind direction. The Coulomb collisions between protons and electrons, and the electron heat flux affect the radial profile of electron and proton temperatures but the effect of Coulomb collision is negligible compared to the turbulent heating term.
• The theoretical and observed electron and proton entropy increase as a function of heliocentric distance. In the upwind direction, the electron and proton entropy increases by about 12.48% and 17.24% from 0.17 au to 75 au, respectively. In the downwind direction from 0.17 to 75 au, the electron entropy increases by about 9.11% and the proton entropy by about 12.5%. The entropy in the upwind direction is higher than that in the downwind direction.
• In the upwind direction, the theoretical and observed energy in forward and backward propagating modes decreases gradually until ∼10–20 au, and then slightly increases as distance increases. In the downwind direction, the theoretical energy in forward and backward propagating modes decreases monotonically until ∼40–50 au, and then slightly increases until 75 au.
• The fluctuating magnetic energy density in the upwind direction decreases as $r^{-2.33}$, while that in the downwind direction decreases as $r^{-2.56}$. Similarly, in the upwind direction, the fluctuating kinetic energy decreases gradually until ∼10 au, and then increases as distance increases. However, in the downwind direction, the fluctuating kinetic energy decreases monotonically until ∼40 au, and then increases as a function of heliocentric distance.
• In the upwind direction, the fluctuating kinetic and magnetic energies tend towards equipartition beyond ∼10 au. In the downwind direction, the normalized residual energy increases after ∼30 au, but the fluctuating magnetic and kinetic energy do not balance. The normalized cross-helicity in the upwind direction is approximately zero beyond ∼20 au, but not in the downwind direction.

B. Turbulence cascade rate in the upwind and downwind directions:

• The turbulence cascade rate of energy in forward and backward propagating modes and the fluctuating magnetic energy decrease gradually until ∼10 au and ∼40 au in the upwind and downwind directions, respectively, and then increases slightly with increasing heliocentric distance.
• Over the heliocentric distance ∼0.2–10 au, the downwind turbulence cascade rate is larger than the upwind turbulence cascade rate. However, the turbulence cascade rate in the upwind direction is larger than that in the downwind direction beyond 10 au when pickup ions begin to influence the solar wind.
• The turbulence cascade rate of the normalized residual energy ${ϵ}_{{\sigma }_D}$ obtained from model 1 is negative throughout the heliosphere in the upwind and downwind directions.
• The turbulence cascade rate of the normalized cross-helicity in the downwind direction is larger than that in the upwind direction.

C. CR mfps in the upwind and downwind directions:

• The cosmic ray parallel mfp ${\mathsf{\Lambda }}_{\left|\right|}$ dominates the perpendicular mfp ${\mathsf{\Lambda }}_{\perp }$ throughout the heliosphere in the upwind and downwind directions. The CR parallel mfp in the upwind direction is larger than that in the downwind direction between ∼0.2–10 au. The CR parallel mfp in the downwind direction is larger than that in the upwind direction when pickup ions begin to influence the solar wind.
• The CR perpendicular mfp ${\mathsf{\Lambda }}_{\perp }$ increases monotonically as a function of heliocentric distance in the upwind and downwind directions. In the downwind direction, the CR perpendicular mfp is larger than that in the upwind direction from 0.17 au to 75 au.
• In the upwind direction, the CR radial mfp ${\mathsf{\Lambda }}_{rr}$ increases initially until ∼1–2 au, and then decreases gradually until ∼10 au. When pickup ions begin to affect the solar wind beyond ∼10 au, the ${\mathsf{\Lambda }}_{rr}$ decreases more rapidly. In the downwind direction, ${\mathsf{\Lambda }}_{rr}$ remains approximately constant between ∼0.1–10 au, and then decreases with increasing heliocentric distance.

We studied solar wind turbulence in the upwind and downwind directions in detail by using an incompressible MHD turbulence transport model [77] and spacecraft measurements. The results show that solar wind turbulence characterized by turbulence energies and correlation lengths, is different in the upwind than downwind direction. The different turbulence properties in the upwind and downwind directions yield different turbulence cascade rates and CR diffusion tensors in these two directions. The latter point is of particular importance to models describing cosmic ray modulation since typically the CR diffusion tensor has been assumed to be the same in both the upwind and downwind directions. Thus, besides the difference in the upwind and downwind distances to the heliospheric termination shock e.g., refs. [94,198], the character of the turbulence and hence CR scattering mfps is different.