Theory and Transport of Nearly Incompressible Magnetohydrodynamic Turbulence

The theory of nearly incompressible magnetohydrodynamics (NI MHD) was developed in the early 1990s (Klainerman & Majda 1981, 1982; Montgomery et al. 1987; Matthaeus & Brown 1988; Zank & Matthaeus 1990, 1991, 1992a, 1992b, 1993; Zank et al. 1990; Matthaeus et al. 1991; Bayly et al. 1992; Ghosh & Matthaeus 1992) for homogeneous flows, and subsequently extended to large-scale inhomogeneous flows (Bhattacharjee et al. 1998; Hunana et al. 2006, 2008; Hunana & Zank 2010). NI MHD is a formulation of the MHD equations in a weakly compressible or nearly incompressible regime. Much of the focus in the earlier work was to understand the apparent incompressibility of the solar wind and other plasma environments, and the relationship of density fluctuations to apparently incompressible manifestations of turbulence in the solar wind and interstellar medium. Specifically, NI MHD was developed initially to explain why the observed interstellar electron density fluctuation spectrum (Armstrong & Woo 1980; Armstrong et al. 1981, 1995; Spangler & Armstrong 1990) appears to follow a ${k}^{-5/3}$ Kolmogorov-like spectrum (Montgomery et al. 1987). Besides weak compressibility, NI MHD has been applied to explain a variety of other solar wind observations, including perhaps most importantly the prediction that solar wind turbulence in plasma beta O(1) or $\ll 1$ regions is a superposition of dominant 2D and minority slab components (Zank & Matthaeus 1992b, 1993).

The NI MHD description provides a physically and mathematically consistent coupling of low-frequency incompressible and low-frequency density fluctuations, and therefore provides a natural framework within which to investigate fundamental questions about low-frequency turbulence, particularly its dimensionality, and compressible turbulence in the solar wind. Despite this, NI MHD has not been applied systematically to the study of turbulence in general and in the solar wind specifically. Part of the difficulty in using the NI MHD formulation lies in its somewhat complicated expression in terms of normalized variables, necessary for the systematic derivation of the NI MHD model, but unnecessary otherwise. By means of an Elsässer variable formulation of NI MHD, distinct descriptions of dominant 2D and minority slab turbulence emerge naturally in the plasma beta O(1) or $\ll 1$ regimes, as do the nonlinear couplings between slab and 2D modes. We examine NI MHD in both the homogeneous and inhomogeneous formulations, the latter being appropriate to the large-scale solar wind and solar corona. The Elsässer formulation allows us to (1) address some basic questions regarding the characteristics and spectra of fluctuations in locally homogeneous flows, (2) formulate a transport description for the Elsässer variables and derived moments or ensemble averages for the inhomogeneous solar wind, and (3) describe density fluctuation spectra and their transport in the solar wind.

The paper is organized as follows. In Section 2, we begin with a brief review of NI MHD, after which we express the homogeneous and inhomogeneous NI MHD equations in non-normalized form. Thereafter, in Section 3, we discuss certain properties of the homogeneous NI MHD equations, derive the Elsässer form of the equations, and discuss the nature of the propagating and advected fluctuations and spectra. Section 4 addresses the Elsässer formulation of the inhomogeneous NI MHD equations in the core incompressible limit. Besides deriving the transport equations for the appropriate Elsässer variables, we derive a system of transport equations that describes basic turbulence variables such as the energy density in forward and backward modes, the residual energy, cross-helicity, and the corresponding correlation lengths, following the approach of Zank et al. (1996, 2012a), Matthaeus et al. (1996), Breech et al. (2005, 2008), Oughton et al. (2006, 2011), Adhikari et al. (2015a), and Wiengarten et al. (2015, 2016). Section 5 derives the transport equations for the NI Elsässer variables in an inhomogeneous flow, from which the corresponding NI transport equations for the NI energy density in forward and backward modes, the residual energy, the cross-helicity, and the corresponding correlation lengths are derived. Solutions to the full system of radially symmetric incompressible and NI transport equations for turbulence in the supersonic solar wind are presented in Section 6. The transport of nearly incompressible density fluctuations in the inhomogeneous solar wind is discussed in Sections 4 and 6 as well. The conclusions, presented in Section 7, summarize the main results.

NI MHD describes a magnetofluid that is in a weakly compressible regime. Observations of the solar wind show that density fluctuations have typical amplitudes that deviate no more than 10% from the mean solar wind density (Matthaeus et al. 1991; Bavassano & Bruno 1995; Bruno & Carbone 2013). Numerous observations also show that the solar wind behaves in many respects as an incompressible magnetofluid, with magnetic field and velocity fluctuations being well described by incompressible turbulence (Bruno & Carbone 2013). For example, power law distributions for solar wind fluctuations (energy, density, temperature, magnetic field) with a Kolmogorov exponent ($-5/3$) or occasionally an Iroshnikov–Kraichnan exponent ($-3/2$) are observed frequently. Such observations suggest that compressible MHD converges in some sense to an incompressible state, passing through a nearly incompressible regime in which weakly compressible corrections are related to the final incompressible state. A critical scaling parameter in the compressible MHD equations is the inverse square of the "turbulent Mach number" (i.e., a Mach number defined in terms of the fluctuating velocity and mean sound speed rather than in terms of the bulk flow velocity). The turbulent Mach number is typically very small and therefore introduces a highly singular term into the MHD equations that leads to the generation of high-frequency (magneto)acoustic fluctuations. Klainerman & Majda (1981, 1982) first described the conditions under which solutions of the fully compressible hydrodynamic equations converge to solutions of the incompressible equations in the limit where the turbulent Mach number approaches zero. They showed rigorously that a small Mach number alone is not sufficient for convergence, and that the initial conditions must comprise incompressible initial data with a superposition of small fluctuations that are of the order of the Mach number O(M) for velocity fluctuations and of the order $O({M}^{2})$ for pressure fluctuations (see also Hunana et al. 2006 and Matthaeus & Brown 1988). Ghosh & Matthaeus (1992) used fully compressible numerical simulations to explore the approach to incompressibility in 2D low Mach number hydrodynamic turbulence. By varying the initial conditions and Mach number, they found that departures from incompressibility depended strongly on the assumed initial data. Flows with a prescribed solenoidal velocity and without initial acoustic fluctuations evolved nearly incompressibly for turbulent Mach numbers as large as M = 0.5. Including initial acoustic fluctuations with a solenoidal velocity yielded solutions that began to depart from incompressibility. However, for initial conditions comprising either random initial velocity fluctuations or acoustic waves, the hydrodynamic system evolved far from an incompressible state, even when the turbulent Mach number was chosen to be as small as M = 0.1. Similar results were obtained by Passot & Pouquet (1987), Shaikh & Zank (2006), and Ghosh & Parashar (2015) in the context of MHD.

Evidently, complications are introduced in MHD because the presence of three wave speeds modifies the definition of the turbulent Mach number, and consequently one needs to specify the plasma beta regime that is of interest. This makes the formal analysis of Klainerman & Majda (1981, 1982) more challenging to implement in the presence of a magnetic field, a more general equation of state (including, e.g., heat conduction), and dissipation. A more intuitive physical approach to the derivation of nearly incompressible MHD was developed by Zank & Matthaeus (1990, 1991, 1992a, 1992b, 1993), and Zank et al. (1990) for homogeneous flows. NI MHD predicted a variety of phenomena that were subsequently observed in the solar wind (see, e.g., Zank et al. 1990, Matthaeus et al. 1991, and Klein et al. 1993 for the predictions related to the correlation of various fluctuations, and Zank & Matthaeus 1992b for a discussion about the prediction that solar wind turbulence is a superposition of dominant 2D and slab components). However, observational studies by Tu & Marsch (1994), Bavassano et al. (1995), Bavassano & Bruno (1995), and Klein et al. (1993) indicated that the NI MHD predicted $O({M}^{2})$ scaling for the amplitude of density fluctuations was not often met in the solar wind. Bhattacharjee et al. (1998) suggested that if the inhomogeneity of the background magnetic field in the solar wind was accounted for, then the amplitude of the density fluctuations became O(M) instead. However, Bhattacharjee et al. (1998) did not develop the nearly incompressible corrections, nor did they include the large-scale flow inhomogeneity. Hunana et al. (2006, 2008) and Hunana & Zank (2010) extended the homogeneous NI MHD theory of Zank & Matthaeus (1993) by assuming an inhomogeneous background flow in equilibrium that carries a magnetic field. Hunana et al. (2006) and Hunana & Zank (2010) predicted O(M) scaling of density fluctuations, consistent with that of Bhattacharjee et al. (1998), using a more general treatment.

It is important to realize that the solar wind is observed to be primarily in a plasma beta regime $\beta \sim 1$, and in and near the solar corona, $\beta \ll 1$ (see reviews by, e.g., Goldstein et al. 1995; Tu & Marsch 1995; Zhou et al. 2004; Bruno & Carbone 2013). The regime $\beta \gg 1$ is therefore not generally applicable to the solar wind. In particular, Bruno & Bavassano (1991) and Bruno & Bavassano (1993) have shown that most Alfvénic periods in the solar wind are characterized by low beta values and low magnetic and density compressibility. Furthermore, the power anisotropy differs typically between slow ($\beta \gt 1$) and fast ($\beta \lt 1$) solar wind, with the former being much more isotropic (Klein et al. 1993). However, as discussed by Lighthill (1978) and derived in the nearly incompressible framework by Zank & Matthaeus (1993), the most widely used incompressible MHD description is valid only for $\beta \gg 1$. Formally, this description should not therefore be used to describe solar wind turbulence, but instead incompressible descriptions based on the $\beta \sim 1$ and $\beta \ll 1$ regimes should be used, depending on the assumed local plasma beta.⁶ Zank & Matthaeus (1993; hereafter ZM93) showed that because nearly incompressible theory is based on an expansion of the normalized equations and collecting terms of similar order, it is necessary in the MHD regime to consider three plasma beta β regimes ($\beta =P/({B}^{2}/2{\mu }_{0})$, where P is the thermal pressure, $B=| {\boldsymbol{B}}|$, ${\boldsymbol{B}}$ the magnetic field, and ${\mu }_{0}$ is the magnetic permeability). ZM93 considered $\beta \ll 1$, β ∼ 1, and $\beta \gg 1$ and showed that the description of solar wind turbulence is very different in these three different regimes. The $\beta \gg 1$ regime yields a leading-order incompressible MHD description that is fully 3D, whereas the leading-order description for the $\beta \ll 1$, $\beta \sim 1$ regimes is reduced to two dimensions in the plane perpendicular to the mean magnetic field. Higher-order nearly incompressible fluctuations for $\beta \sim 1$ comprise 3D propagating magnetosonic and sound waves, and Alfvén waves propagating parallel to the magnetic field. When $\beta \ll 1$, nearly incompressible fluctuations propagate parallel to the magnetic field. Thus ZM93 predicted that in such plasma beta regimes, solar wind turbulence is a superposition of 2D and slab turbulence, dominated by the 2D component. Ghosh & Parashar (2015) present fully compressible 3D MHD simulations that suggest a decoupling of turbulent dynamics between the plane perpendicular to the mean magnetic field and the dynamics along the mean field direction. As noted by Ghosh & Parashar (2015), their results are "tantalizingly consistent" with the predictions of a superposition of dominant 2D and slab turbulence by Zank & Matthaeus (1993).

The 2D-slab turbulence model is now used for cosmic ray transport and modulation in the heliosphere (e.g., Bieber et al. 1994, 1996; Florinski et al. 2003; Shalchi et al. 2004; Zank et al. 2004; Shalchi 2009; Shalchi & Dosch 2008) and particle acceleration at quasi-perpendicular interplanetary shock waves (Zank et al. 2006; Dosch & Shalchi 2010), for example.

2.1. NI MHD Equations for a Homogeneous $\beta \sim 1$ Plasma

Consider first the NI MHD equations for a homogeneous plasma in the presence of a large constant magnetic field ${B}_{0}$ such that $\beta \sim 1$. The NI MHD description comprises a core set of incompressible equations given by the normalized Equations (57)–(59) in ZM93. This set of equations is derived by imposing the condition that all fast-scale magnetoacoustic variation has to vanish. The elimination of fast-timescale variation is accomplished through the principle of bounded derivatives (Kreiss 1980). The bounding of the time derivatives yields the incompressible hydrodynamic or MHD equations as the secular conditions that ensure the elimination of all fast-scale variation (Zank & Matthaeus 1990, 1991). The core equations, despite being derived from the fully 3D compressible MHD equations, are spatially 2D in a pane perpendicular to the mean magnetic field. The weakly compressible equations are given by the normalized Equations (69)–(72) in ZM93, and derived as a perturbative expansion from the core incompressible equations to introduce the effects of weak compressibility. The weakly compressible corrections are fully 3D, unlike the core incompressible equations. Besides the normalized variables, the NI equations in ZM93 include the perturbation parameter ε that is related to the turbulent sonic Mach number and the plasma beta. As described in ZM93, the normalized plasma variables ($\tilde{P}$ denotes thermal pressure, $\tilde{u}$ the plasma velocity, $\tilde{\rho }$ the plasma density, and $\tilde{B}$ the magnetic field) can be expressed in terms of the ansatz,

$$\begin{eqnarray}\tilde{P} & = & {\tilde{P}}_{0}+{\varepsilon }^{2}({\tilde{P}}^{\infty }+{\tilde{P}}^{* });\quad \tilde{{\boldsymbol{u}}}={\tilde{{\boldsymbol{u}}}}^{\infty }+\varepsilon {\tilde{{\boldsymbol{u}}}}_{1};\\ \quad \tilde{\rho } & = & 1+{\varepsilon }^{2}{\tilde{\rho }}^{* };\quad \tilde{{\boldsymbol{B}}}={\tilde{{\boldsymbol{B}}}}_{0}+\varepsilon {\tilde{{\boldsymbol{B}}}}^{\infty }+{\varepsilon }^{2}{\tilde{{\boldsymbol{B}}}}^{* }.\end{eqnarray} \tag{ 1 }$$

Here and henceforth, we will conform to the notation introduced in ZM93 of using an "$\infty$" superscript to denote MHD variables that satisfy the core incompressible equations (i.e., the sound speed is "infinite"), and the superscript "*" or subscript "1" will indicate higher-order corrections. The tilde indicates normalized values. The higher-order corrections can be both compressible and incompressible. The ansatz (1) provides a consistent expansion that reveals several key points. First, the incompressible pressure and the NI pressure correction enter at the same order in the square of the turbulent sonic Mach number, as does the fluctuating density. The incompressible magnetic field component ${{\boldsymbol{B}}}^{\infty }$ is of order the turbulent sonic Mach number, and the NI correction ${{\boldsymbol{B}}}^{* }$ enters at the next order. The NI velocity ${{\boldsymbol{u}}}_{1}$ is also of the order of the turbulent sonic Mach number. While the NI perturbation expansion is best expressed in terms of normalized plasma variables and normalized coordinates, some care has to be exercised in rewriting the ZM93 normalized equations in dimensional form.

Consider first the core Equations (57)–(61) in ZM93. ZM93 introduced characteristic turbulent speed (u₀), time (T), and length (L) scales such that ${u}_{0}T/L=1$ and normalizations t/T, ${\boldsymbol{x}}/L$, ${{\boldsymbol{u}}}^{\infty }/{u}_{0}$, ${P}^{\infty }/{p}_{0}$ (note that ${\tilde{P}}_{0}\equiv {P}_{0}/{p}_{0}$), ${\boldsymbol{B}}/{B}_{0}$, and so on. To illustrate the dimensionalization procedure, consider the momentum Equation (58) in ZM93. The obvious dimensional form becomes

$$\begin{eqnarray*}\displaystyle \frac{\partial {{\boldsymbol{u}}}^{\infty }}{\partial t}+{{\boldsymbol{u}}}^{\infty }\cdot {{\rm{\nabla }}}_{\perp }{{\boldsymbol{u}}}^{\infty } & = & -\,\displaystyle \frac{{{Lu}}_{0}}{{{Tp}}_{0}}{{\rm{\nabla }}}_{\perp }{P}^{\infty }+\displaystyle \frac{{{Lu}}_{0}}{{{TB}}_{0}^{2}}{{\boldsymbol{B}}}^{\infty }\\ & & \cdot {{\rm{\nabla }}}_{\perp }{{\boldsymbol{B}}}^{\infty }-\displaystyle \frac{{{Lu}}_{0}}{{{TB}}_{0}^{2}}{{\rm{\nabla }}}_{\perp }\left(\displaystyle \frac{1}{2}{{B}^{\infty }}^{2}\right).\end{eqnarray*}$$

Here, ${{\rm{\nabla }}}_{\perp }$ is a 2D gradient operator in terms of the spatial coordinates orthogonal to the mean magnetic field. The turbulent sonic Mach number ${M}_{s0}\equiv {u}_{0}/{c}_{s0}\ll 1$, ${c}_{s0}^{2}\equiv \gamma {p}_{0}/{\rho }_{0}$, and ${\varepsilon }^{2}\equiv \gamma {M}_{s0}^{2}$, where γ is the adiabatic index of the fluid. Since the plasma beta $\beta \sim 1$, one can equate ${\varepsilon }^{2}={M}_{A0}^{2}$, where ${M}_{A0}={u}_{0}/{V}_{A0}$ is the turbulent Alfvén Mach number, ${V}_{A0}\equiv {{\boldsymbol{B}}}_{0}/\sqrt{{\mu }_{0}{\rho }_{0}}$ is the Alfvén speed, and ${\mu }_{0}$ is the permeability of free space. In front of the ${{\rm{\nabla }}}_{\perp }{P}^{\infty }$ term, we use ${{Lu}}_{0}/{{Tp}}_{0}({\varepsilon }^{2}/{\varepsilon }^{2})={{Lu}}_{0}/{{Tp}}_{0}(\gamma {M}_{s0}^{2}/\gamma {M}_{s0}^{2})={\varepsilon }^{2}/{\rho }_{0}$, and for the magnetic field terms, ${{Lu}}_{0}/{{TB}}_{0}^{2}({\varepsilon }^{2}/{\varepsilon }^{2})={{Lu}}_{0}/{{TB}}_{0}^{2}({M}_{A0}^{2}/{M}_{A0}^{2})={\varepsilon }^{2}/({\mu }_{0}{\rho }_{0})$. On renormalizing ${P}^{\infty }={\varepsilon }^{2}{P}^{\infty }$ and ${{\boldsymbol{B}}}^{\infty }=\varepsilon {{\boldsymbol{B}}}^{\infty }$, consistent with the ordering of ZM93, we can rewrite their core incompressible Equations (57)–(61) in the dimensional form

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{\perp }\cdot {{\boldsymbol{u}}}^{\infty }=0;\qquad {{\rm{\nabla }}}_{\perp }\cdot {{\boldsymbol{B}}}^{\infty }=0;\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}\displaystyle \frac{\partial {{\boldsymbol{u}}}^{\infty }}{\partial t}+{{\boldsymbol{u}}}^{\infty }\cdot {{\rm{\nabla }}}_{\perp }{{\boldsymbol{u}}}^{\infty } & = & -\,\displaystyle \frac{1}{{\rho }_{0}}{{\rm{\nabla }}}_{\perp }\left({P}^{\infty }+\displaystyle \frac{1}{2{\mu }_{0}}{{B}^{\infty }}^{2}\right)\\ & & +\displaystyle \frac{1}{{\mu }_{0}{\rho }_{0}}{{\boldsymbol{B}}}^{\infty }\cdot {{\rm{\nabla }}}_{\perp }{{\boldsymbol{B}}}^{\infty };\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {{\boldsymbol{B}}}^{\infty }}{\partial t}+{{\boldsymbol{u}}}^{\infty }\cdot {{\rm{\nabla }}}_{\perp }{{\boldsymbol{B}}}^{\infty }={{\boldsymbol{B}}}^{\infty }\cdot {{\rm{\nabla }}}_{\perp }{{\boldsymbol{u}}}^{\infty }.\end{eqnarray} \tag{ 4 }$$

A key point about Equations (2)–(4) is that the incompressible equations to which compressible 3D MHD collapses in the NI limit are the 2D MHD equations in a plane perpendicular to the mean magnetic field (i.e., ${{\boldsymbol{u}}}^{\infty }={{\boldsymbol{u}}}^{\infty }(x,y)$ and ${{\boldsymbol{B}}}^{\infty }={{\boldsymbol{B}}}^{\infty }(x,y)$, where $\hat{{\boldsymbol{z}}}$ is defined by the large-scale mean magnetic field ${{\boldsymbol{B}}}_{0}={B}_{0}\hat{{\boldsymbol{z}}}$). Note too the absence of propagation effects associated with Alfvén waves. Hence, the leading-order incompressible description for a $\beta \sim O(1)$ (and $\ll 1$) plasma is properly 2D MHD in a plane perpendicular to the mean magnetic field.

By contrast, the NI corrections that emerge at the next order are fully 3D and can be expressed in dimensional form as

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\rho }^{* }}{\partial t}+{{\boldsymbol{u}}}^{\infty }\cdot {\rm{\nabla }}{\rho }^{* }+{\rho }_{0}{\rm{\nabla }}\cdot {{\boldsymbol{u}}}_{1}=0;\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {{\boldsymbol{u}}}_{1}}{\partial t}+{{\boldsymbol{u}}}^{\infty }\cdot {\rm{\nabla }}{{\boldsymbol{u}}}_{1}+{{\boldsymbol{u}}}_{1}\cdot {\rm{\nabla }}{{\boldsymbol{u}}}^{\infty }\\ &&\,=-\displaystyle \frac{1}{{\rho }_{0}}{\rm{\nabla }}{P}^{* }-\displaystyle \frac{1}{{\mu }_{0}{\rho }_{0}}{\rm{\nabla }}({{\boldsymbol{B}}}_{0}\cdot {{\boldsymbol{B}}}^{* }+{{\boldsymbol{B}}}^{\infty }\cdot {{\boldsymbol{B}}}^{* })\\ &&\qquad +\displaystyle \frac{1}{{\mu }_{0}{\rho }_{0}}({{\boldsymbol{B}}}_{0}+{{\boldsymbol{B}}}^{\infty })\cdot {\rm{\nabla }}{{\boldsymbol{B}}}^{* }+\displaystyle \frac{1}{{\mu }_{0}{\rho }_{0}}{{\boldsymbol{B}}}^{* }\cdot {\rm{\nabla }}{{\boldsymbol{B}}}^{\infty };\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {P}^{* }}{\partial t}+{{\boldsymbol{u}}}^{\infty }\cdot {\rm{\nabla }}{P}^{* }+\gamma {p}_{0}{\rm{\nabla }}\cdot {{\boldsymbol{u}}}_{1}=-\displaystyle \frac{\partial {P}^{\infty }}{\partial t}-{{\boldsymbol{u}}}^{\infty }\cdot {\rm{\nabla }}{P}^{\infty };\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {{\boldsymbol{B}}}^{* }}{\partial t}+{{\boldsymbol{u}}}^{\infty }\cdot {\rm{\nabla }}{{\boldsymbol{B}}}^{* }-{{\boldsymbol{B}}}^{* }\cdot {\rm{\nabla }}{{\boldsymbol{u}}}^{\infty }=({{\boldsymbol{B}}}_{0}+{{\boldsymbol{B}}}^{\infty })\cdot {\rm{\nabla }}{{\boldsymbol{u}}}_{1}\\ &&\quad -({{\boldsymbol{B}}}_{0}+{{\boldsymbol{B}}}^{\infty }){\rm{\nabla }}\cdot {{\boldsymbol{u}}}_{1}-{{\boldsymbol{u}}}_{1}\cdot {\rm{\nabla }}{{\boldsymbol{B}}}^{\infty }.\end{eqnarray} \tag{ 8 }$$

The * variables (${\rho }^{* }$, etc.) describe the higher-order corrections to the core equations and, at least for the homogeneous $\beta \sim O(1)$ NI theory, enter as the square of the turbulent sonic or Alfvénic Mach number for density and pressure corrections or linearly in the velocity. As in Equations (2)–(4), we renormalized the variables in (5)–(8) using ${\rho }^{* }={\varepsilon }^{2}{\rho }^{* }$, ${P}^{* }={\varepsilon }^{2}{P}^{* }$, and ${{\boldsymbol{u}}}_{1}=\varepsilon {{\boldsymbol{u}}}_{1}$. Several noteworthy points about (5)–(8) are immediately apparent. Equations (5)–(8) are not linearized but instead describe the driving of the higher-order corrections by the dynamics of the leading-order incompressible flow variables. The incompressible flow ${{\boldsymbol{u}}}^{\infty }$, in particular, introduces an important passive scalar transport component to the dynamics of the higher-order plasma variables. The variables ${\rho }^{* }$ and ${{\boldsymbol{u}}}_{1}$ should not be regarded as exclusively compressible, since they also possess an incompressible component, which we discuss later. The incompressible flow acts as a source of higher-order fluctuations, including compressible fluctuations, and therefore represents a generalization of the well-known Lighthill mechanism for the generation of sound and also a possible origin of density fluctuations. Finally, since the NI corrections enter at $O({M}_{s0}^{2})$, and the 2D component is dominant, turbulent fluctuations in a $\beta \sim 1$ or $\ll 1$ plasma are primarily 2D rather than slab. In the solar wind, this implies that the partitioning of energy between 2D and slab components is typically in the ratio 80:20 (Zank & Matthaeus 1992b), which has been confirmed observationally (Bieber et al. 1996). This partitioning also implies a strong variance anisotropy in the solar wind fluctuations, since the parallel variance is greatly reduced compared to the perpendicular variance, consistent with that observed (e.g., Belcher & Davis 1971; Parashar et al. 2016).

Before examining the NI MHD $\beta \sim 1$ equations more closely, we consider the Hunana & Zank (2010) (HZ10) extension of Equations (2)–(8) to an inhomogeneous $\beta \sim 1$ MHD flow.

2.2. NI MHD Equations for an Inhomogeneous $\beta \sim 1$ Plasma

Hunana & Zank (2010; HZ10) extend ZM93 to a time-stationary large-scale inhomogeneous background flow such as the solar wind. As discussed previously, a fundamental difference between homogeneous and inhomogeneous NI MHD is that density fluctuations enter at $O({M}_{s0})$ in the latter case. The collapse in dimensionality is also weaker in the inhomogeneous case, being locally 2D. Other than these points, expressing the normalized $\beta \sim 1$ equations of HZ10 in dimensional form follows the same procedure given previously. Let the background density ${\rho }_{\mathrm{SW}}={\rho }_{\mathrm{SW}}({\boldsymbol{x}})$, pressure ${P}_{\mathrm{SW}}={P}_{\mathrm{SW}}({\boldsymbol{x}})$, flow velocity ${\boldsymbol{U}}={\boldsymbol{U}}({\boldsymbol{x}})$, and magnetic field ${{\boldsymbol{B}}}_{\mathrm{SW}}={{\boldsymbol{B}}}_{\mathrm{SW}}({\boldsymbol{x}})$ be prescribed stationary functions of position. The core incompressible Equations (55)–(61) in HZ10 derived from the 3D compressible MHD equations are given in dimensional form by

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot {{\boldsymbol{u}}}^{\infty }=-\displaystyle \frac{1}{{\rho }_{\mathrm{SW}}}{{\boldsymbol{u}}}^{\infty }\cdot {\rm{\nabla }}{\rho }_{\mathrm{SW}};\qquad {\rm{\nabla }}\cdot {{\boldsymbol{B}}}^{\infty }=0;\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{\rho }_{\mathrm{SW}}\left(\displaystyle \frac{\partial {{\boldsymbol{u}}}^{\infty }}{\partial t}+({\boldsymbol{U}}+{{\boldsymbol{u}}}^{\infty })\cdot {\rm{\nabla }}{{\boldsymbol{u}}}^{\infty }\right)=-{\rm{\nabla }}{P}^{\infty }\\ &&\quad +\displaystyle \frac{1}{{\mu }_{0}}({\rm{\nabla }}\times {{\boldsymbol{B}}}^{\infty })\times {{\boldsymbol{B}}}^{\infty }-{\rho }_{\mathrm{SW}}{{\boldsymbol{u}}}^{\infty }\cdot {\rm{\nabla }}{\boldsymbol{U}};\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}\displaystyle \frac{\partial {{\boldsymbol{B}}}^{\infty }}{\partial t}+({\boldsymbol{U}}+{{\boldsymbol{u}}}^{\infty })\cdot {\rm{\nabla }}{{\boldsymbol{B}}}^{\infty } & = & {{\boldsymbol{B}}}^{\infty }\cdot {\rm{\nabla }}{{\boldsymbol{u}}}^{\infty }-{{\boldsymbol{B}}}^{\infty }{\rm{\nabla }}\cdot {{\boldsymbol{u}}}^{\infty }\\ & & +{{\boldsymbol{B}}}^{\infty }\cdot {\rm{\nabla }}{\boldsymbol{U}}-{{\boldsymbol{B}}}^{\infty }{\rm{\nabla }}\cdot {\boldsymbol{U}};\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\rho }^{\infty }}{\partial t}+({\boldsymbol{U}}+{{\boldsymbol{u}}}^{\infty })\cdot {\rm{\nabla }}{\rho }^{\infty }+{\rho }^{\infty }{\rm{\nabla }}\cdot {{\boldsymbol{u}}}^{\infty }=-{\rho }^{\infty }{\rm{\nabla }}\cdot {\boldsymbol{U}}.\end{eqnarray} \tag{ 12 }$$

Several points about Equations (9)–(12) should be made. The "incompressible" density ${\rho }^{\infty }$ enters in normalized form at $O({M}_{s0})$ according to $\tilde{\rho }={\tilde{\rho }}_{\mathrm{SW}}+\varepsilon {\tilde{\rho }}^{\infty }+{\varepsilon }^{2}{\tilde{\rho }}^{* }$, and behaves as a passive scalar driven by the dynamics of the incompressible flow field ${{\boldsymbol{u}}}^{\infty }$ according to (12). Also different from the homogeneous case, the incompressible flow velocity ${{\boldsymbol{u}}}^{\infty }$ is not solenoidal because of the large-scale background density variation ${\rho }_{\mathrm{SW}}$. As discussed at some length in HZ10, the NI correction to the thermal pressure enters only at the second order $O({M}_{s0}^{2})$ and the $O({M}_{s0})$ correction ${P}_{1}=0$—that is, the normalized pressure satisfies $\tilde{P}={\tilde{P}}_{\mathrm{SW}}+{\varepsilon }^{2}({\tilde{P}}^{\infty }+{\tilde{P}}^{* })$. Notice too that the large-scale velocity cannot be translated out of the incompressible Equations (9)–(12). However, as discussed in HZ10, the incompressible inhomogeneous MHD equations are locally homogeneous and collapse to 2D in a coordinate system orthogonal to the local mean magnetic field. In this sense, HZ10 describe the incompressible Equations (9)–(12) as weakly 2D. A final important point about the system (9)–(12), derived from the fully compressible 3D MHD equations, is that wave propagation effects, including Alfvén waves, are completely absent, indicating again that the equations are essentially 2D.

The NI corrections enter at $O({M}_{s0}^{2})$ in the density and pressure and $O({M}_{s0})$ in the velocity. As discussed previously, this yields a dominant 2D fluctuating component and weaker slab component, again in the ratio ∼80:20. We can express the HZ10 Equations (76)–(79) in dimensional form. We have extended the HZ10 equations slightly by assuming that the parameter introduced by HZ10 $\chi =L/R$, where L is a characteristic length scale of the fluctuations as before and R a characteristic scale measuring the variation in the inhomogeneity of the background parameters, is not necessarily of $O(\varepsilon )$. The NI order corrections can then be expressed in dimensional form as

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\rho }^{* }}{\partial t}+({\boldsymbol{U}}+{{\boldsymbol{u}}}^{\infty })\cdot {\rm{\nabla }}{\rho }^{* }+{\rho }_{\mathrm{SW}}{\rm{\nabla }}\cdot {{\boldsymbol{u}}}_{1}+{\rho }^{\infty }{\rm{\nabla }}\cdot {{\boldsymbol{u}}}_{1}\\ &&\quad +{{\boldsymbol{u}}}_{1}\cdot {\rm{\nabla }}{\rho }^{\infty }=-{{\boldsymbol{u}}}_{1}\cdot {\rm{\nabla }}{\rho }_{\mathrm{SW}};\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{\rho }_{\mathrm{SW}}\left(\displaystyle \frac{\partial {{\boldsymbol{u}}}_{1}}{\partial t}+({\boldsymbol{U}}+{{\boldsymbol{u}}}^{\infty })\cdot {\rm{\nabla }}{{\boldsymbol{u}}}_{1}\right)\\ &&\quad +{\rho }_{\mathrm{SW}}{{\boldsymbol{u}}}_{1}\cdot {\rm{\nabla }}({\boldsymbol{U}}+{{\boldsymbol{u}}}^{\infty })=-{\rm{\nabla }}{P}^{* }\\ &&\quad +\displaystyle \frac{1}{{\mu }_{0}}({\rm{\nabla }}\times {{\boldsymbol{B}}}^{* })\times ({{\boldsymbol{B}}}^{\infty }+{{\boldsymbol{B}}}_{\mathrm{SW}})+\displaystyle \frac{1}{{\mu }_{0}}({\rm{\nabla }}\times {{\boldsymbol{B}}}^{\infty })\\ &&\quad \times {{\boldsymbol{B}}}^{* }+\displaystyle \frac{1}{{\mu }_{0}}({\rm{\nabla }}\times {{\boldsymbol{B}}}_{\mathrm{SW}})\times {{\boldsymbol{B}}}^{* }\\ &&\quad -{\rho }^{\infty }({\boldsymbol{U}}+{{\boldsymbol{u}}}^{\infty })\cdot {\rm{\nabla }}{\boldsymbol{U}}\\ &&\quad -{\rho }^{\infty }\left(\displaystyle \frac{\partial {{\boldsymbol{u}}}^{\infty }}{\partial t}+({\boldsymbol{U}}+{{\boldsymbol{u}}}^{\infty })\cdot {\rm{\nabla }}{{\boldsymbol{u}}}^{\infty }\right)\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {{\boldsymbol{B}}}^{* }}{\partial t}+({\boldsymbol{U}}+{{\boldsymbol{u}}}^{\infty })\cdot {\rm{\nabla }}{{\boldsymbol{B}}}^{* }+{{\boldsymbol{B}}}^{* }\cdot {{\boldsymbol{u}}}^{\infty }-{\rm{\nabla }}\\ &&\quad \times ({{\boldsymbol{u}}}_{1}\times {{\boldsymbol{B}}}^{\infty })={{\boldsymbol{B}}}_{\mathrm{SW}}\cdot {\rm{\nabla }}{{\boldsymbol{u}}}_{1}-{{\boldsymbol{B}}}_{\mathrm{SW}}{\rm{\nabla }}\cdot {{\boldsymbol{u}}}_{1}\\ &&\quad -{{\boldsymbol{u}}}_{1}\cdot {\rm{\nabla }}{{\boldsymbol{B}}}_{\mathrm{SW}}+{{\boldsymbol{B}}}^{* }\cdot {\rm{\nabla }}{\boldsymbol{U}}-{{\boldsymbol{B}}}_{\mathrm{SW}}{\rm{\nabla }}\cdot {\boldsymbol{U}}\\ &&\quad +{{\boldsymbol{B}}}^{* }\displaystyle \frac{1}{{\rho }_{\mathrm{SW}}}{{\boldsymbol{u}}}^{\infty }\cdot {\rm{\nabla }}{\rho }_{\mathrm{SW}};\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {P}^{* }}{\partial t}+({\boldsymbol{U}}+{{\boldsymbol{u}}}^{\infty })\cdot {\rm{\nabla }}{P}^{* }+\gamma {P}_{\mathrm{SW}}{\rm{\nabla }}\cdot {{\boldsymbol{u}}}_{1}+{{\boldsymbol{u}}}_{1}\\ &&\quad \cdot {\rm{\nabla }}{P}_{\mathrm{SW}}=-\displaystyle \frac{\partial {P}^{\infty }}{\partial t}-({\boldsymbol{U}}+{{\boldsymbol{u}}}^{\infty })\cdot {\rm{\nabla }}{P}^{\infty }.\end{eqnarray} \tag{ 16 }$$

Many of the same points as made for the homogeneous NI corrections can be made here and are not repeated.

In this section, we consider properties and solutions of the homogeneous NI MHD $\beta \sim 1$ equations. The 2D core incompressible Equations (2)–(4) can be rewritten in terms of the core Elsässer variables:

$$\begin{eqnarray}&&{{\boldsymbol{z}}}^{\infty \pm }={{\boldsymbol{u}}}^{\infty }\pm \displaystyle \frac{{{\boldsymbol{B}}}^{\infty }}{\sqrt{{\mu }_{0}{\rho }_{0}}}.\end{eqnarray} \tag{ 17 }$$

In so doing, we may rewrite the coupled incompressible core Equations (2)–(4) as

$$\begin{eqnarray}&&\displaystyle \frac{\partial {{\boldsymbol{z}}}^{\infty \pm }}{\partial t}+{{\boldsymbol{z}}}^{\infty \mp }\cdot {{\rm{\nabla }}}_{\perp }{{\boldsymbol{z}}}^{\infty \pm }=-\displaystyle \frac{1}{{\rho }_{0}}{{\rm{\nabla }}}_{\perp }\left({P}^{\infty }+\displaystyle \frac{1}{2{\mu }_{0}}{{B}^{\infty }}^{2}\right),\end{eqnarray} \tag{ 18 }$$

where, as a reminder, ${{\rm{\nabla }}}_{\perp }$ refers to a coordinate system orthogonal to the mean magnetic field that is oriented along the $\hat{{\boldsymbol{z}}}$-axis. Despite the presence of a large mean magnetic field ($\beta \sim 1$), the dominant fluctuating component is 2D, and propagation effects are completely absent, as illustrated in (18). Equation (18) indicates that all core modes have zero frequency and the interactions are purely "nonlinear" (i.e., there is no mediation by Alfvénic wave packets; Fyfe & Montgomery 1976; Fyfe et al. 1977).

Equations (18) can be analyzed using a von Karman–Howarth approach (von Karman & Howarth 1938; Matthaeus et al. 1994; Zank et al. 1996, 2012a). We introduce the nonlinear timescales

$$\begin{eqnarray}&&{({\tau }_{\infty }^{\pm })}^{-1}\equiv \displaystyle \frac{{\langle {{z}^{\infty \mp }}^{2}\rangle }^{1/2}}{{\lambda }_{\perp }^{\pm }}\end{eqnarray} \tag{ 19 }$$

(Pouquet et al. 1976; Dobrowolny et al. 1980a, 1980b; Grappin et al. 1982, 1983; Matthaeus & Zhou 1989), and approximate the quadratic nonlinearity in (18) by

$$\begin{eqnarray}&&{{\boldsymbol{NL}}}^{\infty \pm }\equiv {{\boldsymbol{z}}}^{\infty \pm }\displaystyle \frac{{\langle {{z}^{\infty \mp }}^{2}\rangle }^{1/2}}{{\lambda }_{\perp }^{\pm }}.\end{eqnarray} \tag{ 20 }$$

The length scale ${\lambda }_{\perp }^{\pm }$ will correspond essentially to the correlation length in the 2D "forward" and "backward" modes. On introducing the ensemble-averaging operator $\langle \bullet \rangle$ and the total core energy in fluctuations

$$\begin{eqnarray}&&\langle {{z}^{\infty }}^{2}\rangle \equiv \displaystyle \frac{\langle {{z}^{\infty +}}^{2}\rangle +\langle {{z}^{\infty -}}^{2}\rangle }{2},\end{eqnarray} \tag{ 21 }$$

we make the reasonable assumption for 2D turbulence that

$$\begin{eqnarray}&&\displaystyle \frac{{\langle {{z}^{\infty -}}^{2}\rangle }^{1/2}}{{\lambda }_{\perp }^{+}}=\displaystyle \frac{{\langle {{z}^{\infty +}}^{2}\rangle }^{1/2}}{{\lambda }_{\perp }^{-}}\equiv \displaystyle \frac{{\langle {{z}^{\infty }}^{2}\rangle }^{1/2}}{2{\lambda }_{\perp }^{\infty }},\end{eqnarray} \tag{ 22 }$$

which serves as the definition for ${\lambda }_{\perp }^{\infty }$. Equation (18) then yields the ensemble-averaged 1-point energy-containing transport equation for the total energy as

$$\begin{eqnarray}&&\displaystyle \frac{d}{{dt}}\langle {{z}^{\infty }}^{2}\rangle =-\displaystyle \frac{{\langle {{z}^{\infty }}^{2}\rangle }^{3/2}}{{\lambda }_{\perp }^{\infty }}.\end{eqnarray} \tag{ 23 }$$

If we make the usual assumption that the steady energy transfer flux ${\rm{\Pi }}({\boldsymbol{k}})$ (${\boldsymbol{k}}$ is the wave number vector) is proportional to the timescale for the decay of the transfer function correlations ${\tau }_{3}$ (also known as the triple correlation timescale) and that the energy-containing eddies determine the rate of spectral energy transfer in decaying turbulence, we may adopt the phenomenological steady energy transfer rate

$$\begin{eqnarray}&&\epsilon ={\rm{\Pi }}({\boldsymbol{k}})={\tau }_{3}\displaystyle \frac{\langle {{z}^{\infty }}^{2}\rangle }{{\tau }_{\mathrm{nl}}^{2}},\end{eqnarray} \tag{ 24 }$$

where is the dissipation rate (Zhou et al. 2004). The nonlinear dynamical timescale is denoted by ${\tau }_{\mathrm{nl}}$, and from (23) can be expressed as ${\lambda }_{\perp }^{\infty }/{\langle {{z}^{\infty }}^{2}\rangle }^{1/2}$. By assuming isotropy in the 2D plane orthogonal to the mean ${\boldsymbol{B}}$ and writing $\langle {{z}^{\infty }}^{2}\rangle =\int {E}^{\infty }{{dk}}_{\perp }$, it follows at once from (24) that

$$\begin{eqnarray}&&{E}^{\infty }={\epsilon }_{\infty }^{2/3}{k}_{\perp }^{-5/3},\end{eqnarray} \tag{ 25 }$$

where ${k}_{\perp }=| {{\boldsymbol{k}}}_{\perp }| =| ({k}_{x},{k}_{y})|$ is the length of the wave number vector in the 2D plane orthogonal to ${B}_{0}\hat{{\boldsymbol{z}}}$, and ${\epsilon }_{\infty }$ is the dissipation rate for 2D incompressible MHD turbulence. Thus the core energy spectrum expressed in Elsässer variables possesses a Kolmogorov form in terms of the 2D perpendicular wave number k_⊥ (i.e., $\propto {k}_{\perp }^{-5/3}$). Furthermore, as discussed previously, 2D fluctuations represent the dominant component in the NI MHD theory.⁷

Further insight into the nature of the zero-frequency fluctuations admitted by the core 2D MHD Equations (2)–(4) can be gained by considering a class of exact nonlinear solutions. By introducing the vector potential ${{\boldsymbol{A}}}^{\infty }$ through ${{\boldsymbol{B}}}^{\infty }\equiv {\rm{\nabla }}\times {{\boldsymbol{A}}}^{\infty }={\rm{\nabla }}\times {A}^{\infty }\hat{z}$, one can show (Appendix) that exact "vortex" solutions exist, given by

$$\begin{eqnarray}{A}^{\infty }(r)=\left\{\begin{array}{ll}{{CJ}}_{0}\left(\tfrac{k}{\sqrt{{\mu }_{0}\rho }}r\right) & r\lt {r}_{0}\\ 0 & r\geqslant {r}_{0}\end{array}\right..\end{eqnarray} \tag{ 26 }$$

Here J₀ is the Bessel function of order 0, r is the polar coordinate distance with origin at the center of the vortex structure, r₀ is the edge of the vortex, and C and k are constants. The vorticity ${\xi }^{\infty }\equiv {\rm{\nabla }}\times {{\boldsymbol{u}}}^{\infty }$ and current ${\mu }_{0}{{\boldsymbol{J}}}^{\infty }\equiv {\rm{\nabla }}\times {{\boldsymbol{B}}}^{\infty }$ are aligned and parallel to the large-scale magnetic field ${{\boldsymbol{B}}}_{0}$. Since the contours of ${{\boldsymbol{A}}}^{\infty }$ determined by (26) correspond to closed magnetic field lines, the vortex solutions (26) are a subclass of magnetic islands in the 2D plane perpendicular to the mean field ${{\boldsymbol{B}}}_{0}$.

The existence of 2D magnetic islands or vortex structures is therefore an explicit prediction of $\beta \sim 1$ or $\ll 1$ NI MHD. Observations of magnetic islands on a variety of scales in the supersonic solar wind have been reported, ranging from hourly timescales (Khabarova et al. 2015a, 2015b, 2016) to minute timescales (Verkhoglyadova et al. 2003) to ion kinetic scales (Lion et al. 2016). In the last case, the plasma beta of the solar wind environment in which vortex/magnetic island structures were embedded was $\beta \lt 1$ (∼0.5–0.7).

Consider now the minority fluctuating component described by the NI corrections (5)–(8). Before developing an Elsässer description of the NI fluctuations, it is illuminating to consider a "linearized" form of Equations (5)–(8), in which the nonlinear coupling terms between the core and nearly incompressible corrections, such as ${{\boldsymbol{u}}}^{\infty }\cdot {\rm{\nabla }}{{\boldsymbol{u}}}_{1}$ and so on, are neglected. We neglect also the source terms that arise from driving by the low-frequency core plasma variables. In this case, Equations (5)–(8) reduce to a linear system in the NI variables

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\rho }^{* }}{\partial t}+{\rho }_{0}{\rm{\nabla }}\cdot {{\boldsymbol{u}}}_{1}=0;\end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {{\boldsymbol{u}}}_{1}}{\partial t}=-\displaystyle \frac{1}{{\rho }_{0}}{\rm{\nabla }}{P}^{* }-\displaystyle \frac{1}{{\mu }_{0}{\rho }_{0}}({\rm{\nabla }}\times {{\boldsymbol{B}}}^{* })\times {{\boldsymbol{B}}}_{0};\end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {{\boldsymbol{B}}}^{* }}{\partial t}={{\boldsymbol{B}}}_{0}{\rm{\nabla }}\cdot {{\boldsymbol{u}}}_{1}-{{\boldsymbol{B}}}_{0}{\rm{\nabla }}\cdot {{\boldsymbol{u}}}_{1}={\rm{\nabla }}\times ({{\boldsymbol{u}}}_{1}\times {{\boldsymbol{B}}}_{0});\end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {P}^{* }}{\partial t}+\gamma {p}_{0}{\rm{\nabla }}\cdot {{\boldsymbol{u}}}_{1}=0.\end{eqnarray} \tag{ 30 }$$

Notice that, strictly speaking, Equations (27)–(30) do not simply represent a linearization of the compressible MHD equations since terms of similar order such as ${P}^{\infty }$ and ${{\boldsymbol{B}}}^{\infty }$ are excluded. The beautiful classical wave mode analysis of Lighthill (1960) can be followed. On assuming that ${{\boldsymbol{B}}}_{0}={B}_{0}\hat{{\boldsymbol{z}}}$ as before, we obtain the wave equation

$$\begin{eqnarray}&&\displaystyle \frac{{\partial }^{2}{\rm{\Delta }}}{\partial {t}^{2}}={C}_{s0}^{2}{{\rm{\nabla }}}^{2}{\rm{\Delta }}+{V}_{A0}^{2}{{\rm{\nabla }}}^{2}({\rm{\Delta }}-{\rm{\Gamma }}),\end{eqnarray} \tag{ 31 }$$

where ${\rm{\Delta }}\equiv {\rm{\nabla }}\cdot {{\boldsymbol{u}}}_{1}$ and ${\rm{\Gamma }}\equiv \partial {{u}_{1}}_{z}/\partial z$. On introducing the z-component of the NI vorticity ${\rm{\nabla }}\times {{\boldsymbol{u}}}_{1}$, $\xi \equiv -\partial {{u}_{1}}_{x}/\partial y+\partial {{u}_{1}}_{y}/\partial x$, we find that

$$\begin{eqnarray}&&\displaystyle \frac{{\partial }^{2}\xi }{\partial {t}^{2}}={V}_{A0}^{2}\displaystyle \frac{{\partial }^{2}\xi }{\partial {z}^{2}},\end{eqnarray} \tag{ 32 }$$

that is, the z-component of the vorticity propagates along the mean magnetic field at the Alfvén speed V_A0. By contrast, the gradient in ${{\boldsymbol{u}}}_{1}$ along the magnetic field direction Γ is coupled with the compressive component of the NI velocity ${{\boldsymbol{u}}}_{1}$, propagating at the sound speed according to

$$\begin{eqnarray}&&\displaystyle \frac{{\partial }^{2}{\rm{\Gamma }}}{\partial {t}^{2}}={C}_{s0}^{2}\displaystyle \frac{{\partial }^{2}{\rm{\Delta }}}{\partial {z}^{2}}.\end{eqnarray} \tag{ 33 }$$

As pointed out by Lighthill (1960), (31)–(33) are sufficient to completely determine ${{\boldsymbol{u}}}_{1}$. However, for our purposes, it is more useful to introduce the x- and y-components of the vorticity

$$\begin{eqnarray*}&&\eta \equiv \displaystyle \frac{\partial {{u}_{1}}_{z}}{\partial y}-\displaystyle \frac{\partial {{u}_{1}}_{y}}{\partial z};\qquad \zeta \equiv \displaystyle \frac{\partial {{u}_{1}}_{x}}{\partial z}-\displaystyle \frac{\partial {{u}_{1}}_{z}}{\partial x}.\end{eqnarray*}$$

It follows then that

$$\begin{eqnarray}&&\displaystyle \frac{{\partial }^{2}\eta }{\partial {t}^{2}}={V}_{A0}^{2}\left(\displaystyle \frac{{\partial }^{2}\eta }{\partial {z}^{2}}-\displaystyle \frac{{\partial }^{2}{\rm{\Delta }}}{\partial y\partial z}\right);\end{eqnarray} \tag{ 34 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{\partial }^{2}\zeta }{\partial {t}^{2}}={V}_{A0}^{2}\left(\displaystyle \frac{{\partial }^{2}\zeta }{\partial {z}^{2}}-\displaystyle \frac{{\partial }^{2}{\rm{\Delta }}}{\partial x\partial z}\right).\end{eqnarray} \tag{ 35 }$$

Unlike the z-component of the vorticity ξ, the x- and y-vortical components are compressive and can propagate obliquely to the mean magnetic field. The NI fluctuations are therefore a superposition of incompressible and compressible modes. Furthermore, if we restrict our attention to incompressible NI fluctuations satisfying ${\rm{\Delta }}=0$ (i.e., ${\rm{\nabla }}\cdot {{\boldsymbol{u}}}_{1}=0$), then (34) and (35) become

$$\begin{eqnarray}&&\displaystyle \frac{{\partial }^{2}\eta }{\partial {t}^{2}}={V}_{A0}^{2}\displaystyle \frac{{\partial }^{2}\eta }{\partial {z}^{2}};\qquad \displaystyle \frac{{\partial }^{2}\zeta }{\partial {t}^{2}}={V}_{A0}^{2}\displaystyle \frac{{\partial }^{2}\zeta }{\partial {z}^{2}},\end{eqnarray} \tag{ 36 }$$

that is, identical to (32). In other words, incompressible NI modes all propagate at the Alfvén speed in the parallel direction. (Of course, Alfvén waves propagate in all directions with the usual $\cos \theta$ dependence with respect to a propagation angle θ relative to the mean magnetic field, but the group velocity and hence direction of energy propagation is only in the parallel direction.) The minority incompressible fluctuating component therefore comprises an admixture of counter-propagating Alfvén modes (i.e., slab turbulence), which are coupled nonlinearly to the zero-frequency, non-propagating 2D fluctuations described by the core Equations (2)–(4). To see this point more clearly, introduce the NI Elsässer variables

$$\begin{eqnarray}&&{{\boldsymbol{z}}}^{* \pm }\equiv {{\boldsymbol{u}}}_{1}\pm {{\boldsymbol{B}}}^{* }/\sqrt{{\mu }_{0}{\rho }_{0}}.\end{eqnarray} \tag{ 37 }$$

On rewriting Equations (5)–(8), we obtain

$$\begin{eqnarray}&&\displaystyle \frac{\partial {{\boldsymbol{z}}}^{* \pm }}{\partial t}\mp {{\boldsymbol{V}}}_{A0}\cdot {\rm{\nabla }}{{\boldsymbol{z}}}^{* \pm }+{{\boldsymbol{z}}}^{\infty \mp }\cdot {\rm{\nabla }}{{\boldsymbol{z}}}^{* \pm }+{{\boldsymbol{z}}}^{* \mp }\cdot {\rm{\nabla }}{{\boldsymbol{z}}}^{\infty \pm }\\ &&\quad =-\displaystyle \frac{1}{{\rho }_{0}}{\rm{\nabla }}\left({P}^{* }+\displaystyle \frac{1}{{\mu }_{0}}{{\boldsymbol{B}}}_{0}\cdot {{\boldsymbol{B}}}^{* }+\displaystyle \frac{1}{{\mu }_{0}}{{\boldsymbol{B}}}^{\infty }\cdot {{\boldsymbol{B}}}^{* }\right)\\ &&\qquad \mp \left({{\boldsymbol{V}}}_{A0}+\displaystyle \frac{{{\boldsymbol{B}}}^{\infty }}{\sqrt{{\mu }_{0}{\rho }_{0}}}\right){\rm{\nabla }}\cdot {{\boldsymbol{u}}}_{1}.\end{eqnarray} \tag{ 38 }$$

For incompressible modes, that is, vortical modes (32) and (36) propagating at the Alfvén speed, ${\rm{\nabla }}\cdot {{\boldsymbol{u}}}_{1}=0$, and (38) reduce to

$$\begin{eqnarray}&&\displaystyle \frac{\partial {{\boldsymbol{z}}}^{* \pm }}{\partial t}\mp {{\boldsymbol{V}}}_{A0}\cdot {\rm{\nabla }}{{\boldsymbol{z}}}^{* \pm }+{{\boldsymbol{z}}}^{\infty \mp }\cdot {\rm{\nabla }}{{\boldsymbol{z}}}^{* \pm }+{{\boldsymbol{z}}}^{* \mp }\cdot {\rm{\nabla }}{{\boldsymbol{z}}}^{\infty \pm }\\ &&\quad =-\displaystyle \frac{1}{{\rho }_{0}}{\rm{\nabla }}\,\left({P}^{* }+\displaystyle \frac{1}{{\mu }_{0}}{{\boldsymbol{B}}}_{0}\cdot {{\boldsymbol{B}}}^{* }+\displaystyle \frac{1}{{\mu }_{0}}{{\boldsymbol{B}}}^{\infty }\cdot {{\boldsymbol{B}}}^{* }\right).\end{eqnarray} \tag{ 39 }$$

Evidently, the Elsässer modes propagate at the Alfvén speed $\pm {V}_{A0}$ along the mean magnetic field. However, most importantly, the dominant nonlinear interaction is due to the coupling of minority Alfvénic fluctuations ${{\boldsymbol{z}}}^{* \pm }$ with the dominant convected 2D fluctuating component ${{\boldsymbol{z}}}^{\infty \pm }$ (i.e., the NI Elsässer variables respond to the core ${{\boldsymbol{z}}}^{\infty \pm }$ fluctuations in a passive sense). Nonlinear couplings ${{\boldsymbol{z}}}^{* \mp }\cdot {\rm{\nabla }}{{\boldsymbol{z}}}^{* \pm }$ and ${{\boldsymbol{z}}}^{* \pm }\cdot {\rm{\nabla }}{{\boldsymbol{z}}}^{* \mp }$ enter only at the next higher order. The immediate implication is that the nonlinear cascade of energy in the low-frequency inertial range for slab fluctuations in a $\beta \sim 1$ or $\ll 1$ plasma is governed primarily by the dynamical timescale of the dominant advected 2D component. Although analogous to the conclusions of Shebalin et al. (1983), who argued that the nonlinear cascade proceeded via the three-wave coupling of counter-propagating Alfvén waves and a zero-frequency mode, their analysis is appropriate only to the $\beta \gg 1$ regime, is based on a weak turbulence model, and assumed that the mean magnetic field lay within the 2D plane. The very clean separation of the 2D and slab coupling that emerges from Equation (39) provides an explicit demonstration that for a $\beta \sim 1$ or $\ll 1$ plasma, the 2D core nonlinear timescale rather than the Alfvénic timescale dominates and therefore determines the basic spectral characteristics of a minority slab component. Since the spectral timescale is associated with the core or dominant 2D component of the fluctuations, it is not a resonant wave coupling that describes how ${{\boldsymbol{z}}}^{* \pm }$ behaves primarily; instead, ${{\boldsymbol{z}}}^{* \pm }$ behaves more like a passive scalar in response to the ${{\boldsymbol{z}}}^{\infty \pm }$ fluctuating field.

However, despite this discussion, we do not neglect the nonlinear coupling terms provided here and in the following sections. We have noted already that the nonlinear terms ${{{\boldsymbol{z}}}^{* }}^{\mp }\cdot {\rm{\nabla }}{{{\boldsymbol{z}}}^{* }}^{\pm }$ enter only at the next order within the NI expansion (5)–(8) or (13)–(16) (i.e., terms such as ${{\boldsymbol{u}}}_{1}\cdot {\rm{\nabla }}{{\boldsymbol{u}}}_{1}$, ${{\boldsymbol{B}}}^{* }\cdot {\rm{\nabla }}{{\boldsymbol{B}}}^{* }$, ${{\boldsymbol{u}}}_{1}\cdot {\rm{\nabla }}{{\boldsymbol{B}}}^{* }$, etc.). A straightforward extension of ZM93 and HZ10, in which terms such as ${{\boldsymbol{u}}}_{1}\cdot {\rm{\nabla }}{{\boldsymbol{u}}}_{1}$ are retained, shows that these terms contribute higher-order but equivalent terms, as those already contained in (39) and the only structurally new terms introduced are the nonlinear terms ${{{\boldsymbol{z}}}^{* }}^{\mp }\cdot {\rm{\nabla }}{{{\boldsymbol{z}}}^{* }}^{\pm }$. In order to derive the "richest" NI evolution equation, we explicitly include the nonlinear terms in (40) despite it formally being of higher order—that is, the $\beta \sim 1$ NI transport equation is given by

$$\begin{eqnarray}&&\displaystyle \frac{\partial {{\boldsymbol{z}}}^{* \pm }}{\partial t}\mp {{\boldsymbol{V}}}_{A0}\cdot {\rm{\nabla }}{{\boldsymbol{z}}}^{* \pm }+{{\boldsymbol{z}}}^{\infty \mp }\cdot {\rm{\nabla }}{{\boldsymbol{z}}}^{* \pm }\\ &&\qquad +{{\boldsymbol{z}}}^{* \mp }\cdot {\rm{\nabla }}{{\boldsymbol{z}}}^{\infty \pm }+{{\boldsymbol{z}}}^{* \mp }\cdot {\rm{\nabla }}{{\boldsymbol{z}}}^{* \pm }\\ &&\quad =-\displaystyle \frac{1}{{\rho }_{0}}{\rm{\nabla }}\left({P}^{* }+\displaystyle \frac{1}{{\mu }_{0}}{{\boldsymbol{B}}}_{0}\cdot {{\boldsymbol{B}}}^{* }+\displaystyle \frac{1}{{\mu }_{0}}{{\boldsymbol{B}}}^{\infty }\cdot {{\boldsymbol{B}}}^{* }\right).\end{eqnarray} \tag{ 40 }$$

The implications of including the higher-order nonlinearity are interesting and revealed clearly if we again neglect the passive scalar terms in (40). Since ${{\boldsymbol{z}}}^{* \mp }\cdot {\rm{\nabla }}{{\boldsymbol{z}}}^{* \pm }$ enters formally at the next order in the NI expansion, we may express ${{\boldsymbol{z}}}^{* \pm }={{\boldsymbol{z}}}_{1}^{* \pm }+{{\boldsymbol{z}}}_{2}^{* \pm }$, where ${{\boldsymbol{z}}}_{2}^{* \pm }$ is of higher order such that $O({{\boldsymbol{z}}}_{2}^{* \pm })=O({{\boldsymbol{z}}}_{1}^{* \mp }\cdot {\rm{\nabla }}{{\boldsymbol{z}}}_{1}^{* \pm })$. Then, neglecting the passive scalar terms and the RHS of (40), we have

$$\begin{eqnarray}&&\displaystyle \frac{\partial {{\boldsymbol{z}}}_{1}^{* \pm }}{\partial t}\mp {{\boldsymbol{V}}}_{A0}\cdot {\rm{\nabla }}{{\boldsymbol{z}}}_{1}^{* \pm }=0;\end{eqnarray} \tag{ 41 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {{\boldsymbol{z}}}_{2}^{* \pm }}{\partial t}\mp {{\boldsymbol{V}}}_{A0}\cdot {\rm{\nabla }}{{\boldsymbol{z}}}_{2}^{* \pm }=-{{\boldsymbol{z}}}_{1}^{* \mp }\cdot {\rm{\nabla }}{{\boldsymbol{z}}}_{1}^{* \pm }.\end{eqnarray} \tag{ 42 }$$

Equation (41) is a linear wave equation with solutions ${{\boldsymbol{z}}}_{1}^{* \pm }\propto \exp [i{\boldsymbol{k}}\cdot {\boldsymbol{x}}-i{\omega }^{\pm }t]$, provided the frequency ω and wave number ${\boldsymbol{k}}$ are related via the dispersion relation ${\omega }^{\pm }=\mp {{\boldsymbol{V}}}_{A0}\cdot {\boldsymbol{k}}$. Equation (42) is a wave equation too with a source term. On expressing

$$\begin{eqnarray*}&&{{\boldsymbol{z}}}_{1}^{* \pm }=\displaystyle \sum _{k}\,{{\boldsymbol{A}}}_{k}^{\pm }{e}^{i{\boldsymbol{k}}\cdot ({\boldsymbol{x}}\pm {{\boldsymbol{V}}}_{A0}t)},\end{eqnarray*}$$

and substituting in the RHS of (42), it is easily shown that the condition to ensure the absence of secular terms is that the resonance conditions

$$\begin{eqnarray}&&{{\boldsymbol{k}}}^{\prime\prime\prime }={{\boldsymbol{k}}}^{\prime }+{{\boldsymbol{k}}}^{\prime\prime };\qquad {{\boldsymbol{k}}}^{\prime\prime\prime }\cdot {{\boldsymbol{V}}}_{A0}=-{{\boldsymbol{k}}}^{\prime }\cdot {{\boldsymbol{V}}}_{A0}+{{\boldsymbol{k}}}^{\prime\prime }\cdot {{\boldsymbol{V}}}_{A0}\end{eqnarray} \tag{ 43 }$$

hold. Condition (43) is identical to the resonance condition found by Shebalin et al. (1983; see also Montgomery 1989, p. 75) from their perturbation analysis of the incompressible 3D MHD ($\beta \gg 1$) equations. The distinction here from Shebalin et al. (1983) is that the resonance condition holds only for the NI Elsässer variables (i.e., the slab component), and the perturbation ordering arises as a natural consequence of the NI expansion. The implication of (43) is of course that

$$\begin{eqnarray}&&{{\boldsymbol{k}}}^{\prime }\cdot {{\boldsymbol{V}}}_{A0}=0\quad \iff \quad {{\boldsymbol{k}}}^{\prime }\,\perp \,{{\boldsymbol{B}}}_{0},\end{eqnarray} \tag{ 44 }$$

meaning that the nonlinear cascade of energy for slab turbulence proceeds via the coupling of forward and backward propagating Alfvénic fluctuations with 2D zero-frequency modes. This is consistent with the standard picture of the cascade process in MHD turbulence (Shebalin et al. 1983; Montgomery & Matthaeus 1995). Within the NI framework, the nonlinear cascade of energy therefore drives 2D modes, which are then transferred to the dissipation range and heat the plasma. Consequently, the NI MHD 2D core equations acquire a source term corresponding to the nonlinear loss term in the slab NI description, thereby coupling the leading-order 2D and higher-order NI correction descriptions. In deriving the transport equations for NI MHD in an inhomogeneous medium, we include a source of 2D turbulence that is generated by the dissipation of slab turbulence.

Let us now consider the slab turbulence spectra derived from (40). The forward and backward intensities $\langle {{z}^{* \pm }}^{2}\rangle$ satisfy the one-point transport equations

$$\begin{eqnarray}\displaystyle \frac{\partial }{\partial t}\langle {{z}^{* \pm }}^{2}\rangle \mp {{\boldsymbol{V}}}_{A0}\cdot {\rm{\nabla }}\langle {{z}^{* \pm }}^{2}\rangle & \simeq & -\,2\displaystyle \frac{\langle {{z}^{* \pm }}^{2}\rangle {\langle {{z}^{\infty \mp }}^{2}\rangle }^{1/2}}{{\lambda }_{\perp }^{\pm }}\\ & & -2\displaystyle \frac{{\langle {{z}^{* \mp }}^{2}\rangle }^{1/2}{\langle {{z}^{* \pm }}^{2}\rangle }^{1/2}}{{\lambda }_{* }^{\pm }},\end{eqnarray} \tag{ 45 }$$

after invoking orthogonality between the NI and core Elsässer variables (i.e., $\langle {{\boldsymbol{z}}}^{* \pm }\cdot {{\boldsymbol{z}}}^{\infty \pm }\rangle \simeq 0$). Let us consider a slightly simpler problem and estimate the total energy spectrum for

$$\begin{eqnarray*}&&\langle {{z}^{* }}^{2}\rangle \equiv \displaystyle \frac{\langle {{z}^{* +}}^{2}\rangle +\langle {{z}^{* -}}^{2}\rangle }{2}.\end{eqnarray*}$$

We again make use of the Kolmogorov phenomenology in assuming a form related to (24), now given instead by

$$\begin{eqnarray}&&{\epsilon }_{* }={\rm{\Pi }}({\boldsymbol{k}})={\tau }_{3}\displaystyle \frac{\langle {{z}^{* }}^{2}\rangle }{{\tau }_{\mathrm{nl}}^{2}},\end{eqnarray} \tag{ 46 }$$

but we need to proceed a little more carefully because of the passive coupling of the 2D core variable $\langle {{z}^{\infty }}^{2}\rangle$ to $\langle {{z}^{* }}^{2}\rangle$. The nonlinear coupling terms enter at the next higher order, as expressed in Equation (40), and this corresponds to ${\tau }_{\mathrm{nl}}$ in (46).

There is no good reason to suppose that $\langle {{z}^{* }}^{2}\rangle$ is isotropic⁸ , so we define

$$\begin{eqnarray}\langle {{z}^{* }}^{2}\rangle & = & \displaystyle \frac{1}{2\pi }\displaystyle \int {E}^{* }({\boldsymbol{k}}){d}^{3}{\boldsymbol{k}}=\displaystyle \frac{1}{2\pi }\displaystyle \int \displaystyle \int \displaystyle \int {E}^{* }({{\boldsymbol{k}}}_{\perp },{k}_{\parallel })d{{\boldsymbol{k}}}_{\perp }{{dk}}_{\parallel }\\ & = & \displaystyle \int \displaystyle \int {E}^{* }({k}_{\perp },{k}_{\parallel }){k}_{\perp }{{dk}}_{\perp }{{dk}}_{\parallel }\sim {E}^{* }({k}_{\perp },{k}_{\parallel }){k}_{\perp }^{2}{k}_{\parallel },\end{eqnarray} \tag{ 47 }$$

on assuming isotropy in the 2D plane. Consider three cases. First, if we assume ${\tau }_{3}^{-1}={\tau }_{\mathrm{nl}}^{-1}={\langle {{z}^{* }}^{2}\rangle }^{1/2}{k}_{\perp }$, then, from (46), we obtain

$$\begin{eqnarray}&&{E}^{* }({k}_{\perp },{k}_{\parallel }){k}_{\perp }^{2}={\epsilon }_{* }^{2/3}{k}_{\perp }^{-2/3}{k}_{\parallel }^{-1},\end{eqnarray} \tag{ 48 }$$

which is a slight generalization of the Kolmogorov spectrum. Second, if instead we assume that ${\tau }_{3}^{-1}={\tau }_{A}^{-1}={V}_{A}{k}_{\parallel }$ and ${\tau }_{\mathrm{nl}}^{-1}={\langle {{z}^{* }}^{2}\rangle }^{1/2}{k}_{\perp }$, we recover the extended Iroshnikov–Kraichnan (IK) spectrum

$$\begin{eqnarray}&&{E}^{* }({k}_{\perp },{k}_{\parallel }){k}_{\perp }^{2}={({\epsilon }_{* }{V}_{A})}^{1/2}{k}_{\perp }^{-1}{k}_{\parallel }^{-1/2}.\end{eqnarray} \tag{ 49 }$$

Third, we do however need to incorporate the dominant passive scalar interaction of $\langle {{z}^{* }}^{2}\rangle$ with $\langle {{z}^{\infty }}^{2}\rangle$. As before,

$$\begin{eqnarray}\langle {{z}^{\infty }}^{2}\rangle & = & \displaystyle \frac{1}{2\pi }\displaystyle \int {\tilde{E}}^{\infty }({{\boldsymbol{k}}}_{\perp })d{{\boldsymbol{k}}}_{\perp }=\displaystyle \int {\tilde{E}}^{\infty }({k}_{\perp }){k}_{\perp }{{dk}}_{\perp }\\ & & \sim {\tilde{E}}^{\infty }({k}_{\perp }){k}_{\perp }^{2}\equiv {E}^{\infty }({k}_{\perp }){k}_{\perp },\end{eqnarray} \tag{ 50 }$$

on assuming isotropy of $\langle {{z}^{\infty }}^{2}\rangle$ in the 2D plane and relating ${\tilde{E}}^{\infty }({k}_{\perp })$ to ${E}^{\infty }({k}_{\perp })$ given by Equation (25). Following the suggestion of Matthaeus & Zhou (1989), Zhou & Matthaeus (1990a), and Zhou et al. (2004), we introduce the triple correlation time as the sum of the relevant inverse timescales, of which the dominant ones are the passive scalar term ${\tau }_{s}^{-1}\equiv {\langle {{z}^{\infty }}^{2}\rangle }^{1/2}{k}_{\perp }={\epsilon }_{\infty }^{1/3}{k}_{\perp }^{-5/6}{k}_{\perp }^{3/2}$ and the Alfvén term ${\tau }_{A}^{-1}={V}_{A}{k}_{\parallel }$—that is,

$$\begin{eqnarray}&&{\tau }_{3}^{-1}={V}_{A}{k}_{\parallel }+{\epsilon }_{\infty }^{1/3}{k}_{\perp }^{2/3}={V}_{A}{k}_{\parallel }\left(1+{\left(\displaystyle \frac{{\epsilon }_{\infty }}{{V}_{A}^{3}}\right)}^{1/3}{k}_{\perp }^{2/3}{k}_{\parallel }^{-1}\right).\end{eqnarray} \tag{ 51 }$$

Notice that

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{\epsilon }_{\infty }}{{V}_{A}^{3}}\right)}^{1/3}{k}_{\perp }^{2/3}{k}_{\parallel }^{-1}\end{eqnarray} \tag{ 52 }$$

is the critical balance parameter identified by Goldreich & Sridhar (1995), although here with a quite different physical interpretation. In our case, if the term (52) is of order 1 ("critical balance"), then the energy flux in wave number space is a consequence of a balance of Alfvén wave sweeping and passive scalar convection by leading-order 2D ${{\boldsymbol{z}}}^{\infty }$ fluctuations.

By exploiting the phenomenological steady-state energy transfer rate Equation (46), we obtain

$$\begin{eqnarray}{E}^{* }({k}_{\perp },{k}_{\parallel }){k}_{\perp }^{2} & = & {({\epsilon }_{* }{V}_{A})}^{1/2}{k}_{\perp }^{-1}{k}_{\parallel }^{-1/2}\\ & & \times {\left(1+{\left(\displaystyle \frac{{\epsilon }_{\infty }}{{V}_{A}^{3}}\right)}^{1/3}{k}_{\perp }^{2/3}{k}_{\parallel }^{-1}\right)}^{1/2}.\end{eqnarray} \tag{ 53 }$$

Observe that if we assume ${\epsilon }_{\infty }^{1/3}{V}_{A}^{-1}{k}_{\perp }^{2/3}{k}_{\parallel }^{-1}\sim {\rm{const.}}$ or $\ll 1$, we immediately recover the IK spectrum (49) from (53). Conversely, if ${\epsilon }_{\infty }^{1/3}{V}_{A}^{-1}{k}_{\perp }^{2/3}{k}_{\parallel }^{-1}\gg 1$, then ${E}^{* }({k}_{\perp },{k}_{\parallel }){k}_{\perp }^{2}\,={\epsilon }_{* }^{1/2}{\epsilon }_{\infty }^{1/6}{k}_{\perp }^{-2/3}{k}_{\parallel }^{-1}$, and we recover the Kolmogorov spectrum from (53).

Since NI MHD is a superposition of the core 2D MHD equations and the NI corrections, the total turbulent energy spectrum may be expressed as $\langle {z}^{2}\rangle \equiv \langle {{z}^{\infty }}^{2}\rangle +\langle {{z}^{* }}^{2}\rangle$. On using $\langle {{z}^{\infty }}^{2}\rangle =\int {\tilde{E}}^{\infty }({{\boldsymbol{k}}}_{\perp })\delta ({k}_{\parallel })d{{\boldsymbol{k}}}_{\perp }{{dk}}_{\parallel }$, we obtain the total spectrum (core plus NI) as

$$\begin{eqnarray}E({k}_{\perp },{k}_{\parallel }){k}_{\perp }^{2} & = & {C}_{\infty }{\epsilon }_{\infty }^{2/3}{k}_{\perp }^{-2/3}{k}_{\parallel }^{-1}+{C}_{* }{({\epsilon }_{* }{V}_{A})}^{1/2}{k}_{\perp }^{-1}{k}_{\parallel }^{-1/2}\\ & & \times {\left(1+{\left(\displaystyle \frac{{\epsilon }_{\infty }}{{V}_{A}^{3}}\right)}^{1/3}{k}_{\perp }^{2/3}{k}_{\parallel }^{-1}\right)}^{1/2},\end{eqnarray} \tag{ 54 }$$

where C_∞ and C_* are constants reflecting the ratio of 2D to slab turbulence.

Illustrated in Figure 1 is a set of figures showing the complete spectrum $E({k}_{\perp },{k}_{\parallel }){k}_{\perp }^{2}$, Equation (54), comprising a superposition of the majority 2D plus the minority slab component. Here the ratio ${C}_{\infty }{\rm{:}}{C}_{* }={\rm{80:20}}$ is assumed. Figures 1(a) and (b) show the 3D spectrum (54) as a function of wave numbers k_⊥ and k_∥, illustrating that the spectrum is neither isotropic nor possesses a simple single power law description in either k_⊥ and k_∥. The latter point is exhibited particularly clearly in Figures 1(c)–(h), which show cuts of $E({k}_{\perp },{k}_{\parallel }){k}_{\perp }^{2}$ along fixed values of k_⊥ or k_∥. The wave numbers k_⊥ and k_∥ are normalized to the parameter ${k}_{A}\equiv {\epsilon }_{\infty }/{V}_{A}^{3}$. Figures 1(c)–(h) each show three curves, the red curves representing the total spectrum along a fixed value of either k_⊥ or k_∥, and the blue and green curves representing the contribution by the dominant 2D and minority slab components, respectively. The majority 2D component always contributes a Kolmogorov power law spectrum ${k}_{\perp }^{-5/3}$, whereas the minority component contributes a "broken power law" distribution along both cuts. The break in the spectrum is determined by the parameter k_A and the values of k_⊥ or k_∥ relative to it. The total spectrum, being the sum of the dominant 2D and minority slab components, is therefore not a power law necessarily but possesses a somewhat complex "broken power law" structure.

Zoom In Zoom Out Reset image size

Broken power law spectra are observed in the solar wind, primarily in density fluctuation and fluctuating magnetic field spectra (e.g., Telloni et al. 2009; Grappin et al. 1990). However, all the density spectra are measured in frequency and not in 2D wave number space, and the spectra illustrated in Figure 1 are for the Elsässer energy $\langle {z}^{2}\rangle$. Therefore a connection between the results discussed here and solar wind observations has yet to be properly established.

Before concluding this section, we can draw some interesting inferences about the NI density spectrum. In the incompressible limit ${\rm{\nabla }}\cdot {{\boldsymbol{u}}}_{1}=0$, the homogeneous density transport Equation (5) reduces to the passive scalar form

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\rho }^{* }}{\partial t}+{{\boldsymbol{u}}}^{\infty }\cdot {\rm{\nabla }}{\rho }^{* }=0.\end{eqnarray} \tag{ 55 }$$

The density fluctuations are purely entropic, and the passive scalar behavior is due to convection by the dominant core 2D incompressible velocity fluctuations. This is in marked contrast to the conjecture by Lithwick & Goldreich (2001) that entropic density modes are "passively mixed by the cascade of shear Alfvén waves"—see also the related discussion by Chandran et al. (2009). The transport equation for the density fluctuation variance $\langle {\rho }^{* }\rangle$ is given by

$$\begin{eqnarray}&&\displaystyle \frac{d}{{dt}}\langle {{\rho }^{* }}^{2}\rangle =-\langle {{\boldsymbol{u}}}^{\infty }\cdot {\rm{\nabla }}{{\rho }^{* }}^{2}\rangle \simeq -\displaystyle \frac{{\langle {{u}^{\infty }}^{2}\rangle }^{1/2}\langle {{\rho }^{* }}^{2}\rangle }{{{\ell }}_{u}},\end{eqnarray} \tag{ 56 }$$

where we have introduced a correlation length ℓ_u appropriate to the core 2D incompressible velocity fluctuations ${{\boldsymbol{u}}}^{\infty }$. Since $\langle {{u}^{\infty }}^{2}\rangle =(\langle {{u}^{\infty +}}^{2}\rangle +\langle {{u}^{\infty -}}^{2}\rangle )/2$ and the residual energy ${E}_{D}^{\infty }\equiv \langle {{\boldsymbol{z}}}^{\infty +}\cdot {{\boldsymbol{z}}}^{\infty -}\rangle =\langle {{u}^{\infty }}^{2}\rangle -\langle {{B}^{\infty }}^{2}/({\mu }_{0}{\rho }_{0})\rangle$, we have

$$\begin{eqnarray}&&\langle {{u}^{\infty }}^{2}\rangle =\displaystyle \frac{\langle {{z}^{\infty }}^{2}\rangle +{E}_{D}^{\infty }}{2}.\end{eqnarray} \tag{ 57 }$$

Consider three cases. First, ${E}_{D}^{\infty }=0$ implies $\langle {{u}^{\infty }}^{2}\rangle \,=\langle {{B}^{\infty }}^{2}/({\mu }_{0}{\rho }_{0})\rangle$ and hence $\langle {{u}^{\infty }}^{2}\rangle =\langle {{z}^{\infty }}^{2}\rangle /2$. Second, the turbulent kinetic energy dominates (i.e., ${E}_{D}^{\infty }\simeq \langle {{u}^{\infty }}^{2}\rangle$ and $\langle {{B}^{\infty }}^{2}/({\mu }_{0}{\rho }_{0})\rangle =0$), which implies $\langle {{u}^{\infty }}^{2}\rangle =\langle {{z}^{\infty }}^{2}\rangle /2$. Third and finally, if the magnetic energy dominates, then $\langle {{u}^{\infty }}^{2}\rangle \simeq 0$, in which case there is no passive scalar convection of the entropic density fluctuation. This is a very interesting result because, as shown by Adhikari et al. (2015a), solar wind turbulence beyond ∼2 au becomes increasingly dominated by the magnetic energy, and by ∼10 au, the normalized residual energy ${\sigma }_{D}=-1$ (see also Zank et al. 1996 and references therein). Thus, within ∼2 au, density turbulence in the solar wind evolves dynamically as a passive scalar, being mixed by the dominant 2D turbulent velocity component. Beyond 2 au, the density turbulence "freezes" into a non-evolving statistical state, until about 10 au, when pickup ion driven turbulence (Zank et al. 1996; Adhikari et al. 2014, 2015a) begins to generate turbulent velocity fluctuations in the outer heliosphere. In this region, the density turbulence "thaws" and begins to evolve dynamically again. This is discussed in more detail later when we extend these results to the inhomogeneous solar wind.

For ${E}_{D}^{\infty }=0$ or kinetic energy dominated MHD turbulence, Equation (56) can be approximated as

$$\begin{eqnarray}&&\displaystyle \frac{d}{{dt}}\langle {{\rho }^{* }}^{2}\rangle \simeq -\displaystyle \frac{{\langle {{z}^{\infty }}^{2}\rangle }^{1/2}\langle {{\rho }^{* }}^{2}\rangle }{2{{\ell }}_{u}}=-\displaystyle \frac{{\epsilon }_{\rho }}{2},\end{eqnarray} \tag{ 58 }$$

where ${\epsilon }_{\rho }$ is the density dissipation rate. For 2D isotropic density turbulence, $\langle {{\rho }^{* }}^{2}\rangle =\int {E}_{\rho }{{dk}}_{\perp }\sim {E}_{\rho }({k}_{\perp }){k}_{\perp }$. On using ${E}^{\infty }({k}_{\perp })={\epsilon }_{\infty }^{2/3}{k}_{\perp }^{-5/3}$, it follows that

$$\begin{eqnarray}&&{E}_{\rho }({k}_{\perp })={\epsilon }_{\rho }{\epsilon }_{\infty }^{-1/3}{k}_{\perp }^{-5/3}.\end{eqnarray} \tag{ 59 }$$

Alternatively, if we do not assume isotropy, we may, as before, express

$$\begin{eqnarray}&&\langle {{\rho }^{* }}^{2}\rangle =\int \int {E}_{\rho }({k}_{\perp },{k}_{\parallel }){k}_{\perp }{{dk}}_{\perp }{{dk}}_{\parallel }\simeq {E}_{\rho }({k}_{\perp },{k}_{\parallel }){k}_{\perp }^{2}{k}_{\parallel }.\end{eqnarray} \tag{ 60 }$$

On using $\langle {{z}^{\infty }}^{2}\rangle \simeq {\tilde{E}}^{\infty }({k}_{\perp }){k}_{\perp }^{2}\,=\,{E}^{\infty }({k}_{\perp }){k}_{\perp }$ as done previously, together with ${\epsilon }_{\rho }={\langle {{z}^{\infty }}^{2}\rangle }^{1/2}{k}_{\perp }{E}_{\rho }({k}_{\perp },{k}_{\parallel }){k}_{\perp }^{2}{k}_{\parallel }$, we obtain

$$\begin{eqnarray}&&{E}_{\rho }({k}_{\perp },{k}_{\parallel }){k}_{\perp }^{2}={\epsilon }_{\rho }{\epsilon }_{\infty }^{-1/3}{k}_{\perp }^{-2/3}{k}_{\parallel }^{-1}.\end{eqnarray} \tag{ 61 }$$

Observations of density spectra in the solar wind are not easily interpreted. In particular, the results presented previously refer to entropic fluctuations advected in the background solar wind flow, and are not associated with compressive modes that certainly exist in the solar wind (i.e., from small-scale compressive waves to large-scale interplanetary shock waves). Hnat et al. (2003, 2005), in their analysis of density fluctuations in the solar wind, caution that compressibility is very likely manifested in solar wind density spectra, despite incompressible MHD describing many features of solar wind fluctuations. Telloni et al. (2009), Marsch & Tu (1990b), Kellogg & Horbury (2005), and Yao et al. (2011) present observations of solar wind power density spectra that are steeper at small wave numbers and then flatten at higher. Sometimes this appears to be correlated with fluctuations in the magnetic field strength as well, but this has not been definitively established, and indeed Kellogg & Horbury (2005) state that "it appears that density fluctuations are relatively independent of magnetic fluctuations. They are even uncorrelated with fluctuations in the magnitude of B, the usual indicator of compressible fluctuations...." However, Bruno et al. (2014) find that the observed high-frequency flattening of the density spectra is less evident with increasing heliocentric distance: at high frequencies, the spectra steepen as the solar wind expands. Bruno et al. (2014) suggest that this may be related to the parametric decay instability (in which a large amplitude outwardly propagating Alfvén wave decays into an inward Alfvén wave at smaller ${\boldsymbol{k}}$ and a compressive wave; Goldstein 1978), since it is less effective at larger distances. If intermittent events are filtered from the density fluctuations, Bruno et al. (2014) find that the high-frequency density spectra exhibit approximately the same spectral slope. Moreover, it has been found that the density intermittent events are temporally well correlated with the magnetic intermittent events (D. Telloni 2016, personal communication), implying a possible association with magnetosonic waves. Thus the results of Bruno et al. (2014) suggest that the removal of density fluctuations associated with compressive fluctuations (done using a local intermittency measure technique) yields spectra that might be more consistent with passive scalar evolution. Clearly, the interpretation of density spectra needs to be revisited in light of the NI theory presented here.

Consider now the Elsässer formulation of the inhomogeneous core NI MHD Equations (9)–(12). We define the fluctuating Elsässer variables as before—that is,

$$\begin{eqnarray}&&{{\boldsymbol{z}}}^{\infty \pm }\equiv {{\boldsymbol{u}}}^{\infty }\pm \displaystyle \frac{{{\boldsymbol{B}}}^{\infty }}{\sqrt{{\mu }_{0}\rho ({\boldsymbol{x}})}}\equiv {{\boldsymbol{u}}}^{\infty }\pm {{\boldsymbol{v}}}_{A}^{\infty },\end{eqnarray} \tag{ 62 }$$

where now the background density $\rho =\rho ({\boldsymbol{x}})$ is taken to be time-stationary and spatially varying. It is straightforward to combine Equations (9)–(11) in terms of the inhomogeneous Elsässer variables ${{\boldsymbol{z}}}^{\infty \pm }$,

$$\begin{eqnarray}&&\displaystyle \frac{\partial {{{\boldsymbol{z}}}^{\infty }}^{\pm }}{\partial t}+{\boldsymbol{U}}\cdot {\rm{\nabla }}{{{\boldsymbol{z}}}^{\infty }}^{\pm }+{{{\boldsymbol{z}}}^{\infty }}^{\mp }\cdot {\rm{\nabla }}{{{\boldsymbol{z}}}^{\infty }}^{\pm }\\ &&\quad +{{{\boldsymbol{z}}}^{\infty }}^{\mp }\cdot {\rm{\nabla }}{\boldsymbol{U}}+\displaystyle \frac{{{{\boldsymbol{z}}}^{\infty }}^{\pm }-{{{\boldsymbol{z}}}^{\infty }}^{\mp }}{4}{\rm{\nabla }}\cdot {\boldsymbol{U}}\\ &&\quad -\displaystyle \frac{{{{\boldsymbol{z}}}^{\infty }}^{\pm }-{{{\boldsymbol{z}}}^{\infty }}^{\mp }}{4}{{{\boldsymbol{z}}}^{\infty }}^{\mp }\displaystyle \frac{1}{\rho }\cdot {\rm{\nabla }}\rho =-\displaystyle \frac{1}{\rho }{\rm{\nabla }}\left({P}^{\infty }+\displaystyle \frac{{{B}^{\infty }}^{2}}{2{\mu }_{0}}\right).\end{eqnarray} \tag{ 63 }$$

Here, we used the identity ${{\boldsymbol{v}}}_{A}^{\infty }=({{\boldsymbol{z}}}^{\infty +}-{{\boldsymbol{z}}}^{\infty -})/2$. Expression (63) is very similar in some respects to the transport equation for the fluctuating Elsässer variables derived by Zhou & Matthaeus (1990b), and Marsch & Tu (1989, 1990a, 1990b). However, that transport equation, which has been used widely in deriving models that describe the evolution of turbulence in the solar wind (Matthaeus et al. 1994, 1999; Zank et al. 1996, 2012a; Smith et al. 2001; Breech et al. 2008; Usmanov et al. 2011; Kryukov et al. 2012; Adhikari et al. 2014, 2015a, 2015b; Usmanov et al. 2014; Shiota et al. 2016), was derived from the standard 3D incompressible MHD equations and is therefore applicable only to a $\beta \gg 1$ plasma. Equation (63) differs from Equation (21) of Zhou & Matthaeus (1990b), and Marsch & Tu (1989), in that (63) describes the convection of locally 2D Elsässer variables and does not contain the Alfvén velocity nor Alfvén wave propagation effects.

Several points worth noting are immediately apparent from (63). The convection terms for the ${{\boldsymbol{z}}}^{\infty \pm }$ Elsässer variables contain only the large-scale background flow velocity ${\boldsymbol{U}}$, and the singular "Alfvén critical point" ${\boldsymbol{U}}=\pm {{\boldsymbol{V}}}_{A0}$ is absent. As discussed by Adhikari (2015), the presence of the Alfvén critical point for flows transitioning from a sub-Alfvénic to super-Alfvénic regime in turbulence transport models based on the $\beta \gg 1$ 3D incompressible MHD description prevents the transmission of turbulent energy through this point (technically, Adhikari 2015 found that the critical point behaves as an attractive node). Equation (63), besides being appropriate to a $\beta \ll 1$ plasma environment such as the solar corona, does not suffer from this problem. Second, unlike the homogeneous case, the coupling of ${{\boldsymbol{z}}}^{\infty \pm }$ fluctuations is via both the nonlinear term (as before) and the large-scale flow field. These latter so-called mixing-expansion-compression-shear (MECS) terms (Marsch & Tu 1989, 1990a, 1990b; Zhou & Matthaeus 1990b; Zank et al. 1996) generate inward/outward (±) Elsässer fluctuations from the other. Recall the argument of Dobrowolny et al. (1980a, 1980b), who showed using dimensional arguments for a homogeneous flow that strong turbulence could never evolve to an asymmetric state such that ${{\boldsymbol{z}}}^{\infty +}\ne {{\boldsymbol{z}}}^{\infty -}\ne 0$, and that an initially asymmetric distribution of ${{\boldsymbol{z}}}^{\infty \pm }$ would evolve toward a state in which one or the other exists only. The coupling via the MECS terms of the fluctuations to the large-scale flow both enables an asymmetric state for the turbulence and prevents the complete elimination of one or the other mode. These statements will be made more precise later when we develop and provide solutions to a turbulence transport model based on (63).

We follow Zank et al. (2012a) in deriving a general set of transport equations from Equation (63) that describe the evolution of various turbulence quantities in an inhomogeneous flow. This follows in the spirit of previously developed energy-containing models that utilize a Kolmogorov or Irishnikov-Kraichnan phenomenology to describe the physics of the inertial range (Matthaeus et al. 1994; Zank et al. 1996; Breech et al. 2008; Ng et al. 2010). Moments of the "forward" and "backward" modes ${{\boldsymbol{z}}}^{\infty \pm }$ are again introduced through

$$\begin{eqnarray*}&&\langle {{{\boldsymbol{z}}}^{\infty \pm }}^{2}\rangle \equiv \langle {{\boldsymbol{z}}}^{\infty \pm }\cdot {{\boldsymbol{z}}}^{\infty \pm }\rangle ;\quad {\rm{and}}\quad {E}_{D}^{\infty }\equiv \langle {{\boldsymbol{z}}}^{\infty +}\cdot {{\boldsymbol{z}}}^{\infty -}\rangle ,\end{eqnarray*}$$

and $\langle \bullet \rangle$ denotes an averaging over the small scales. The quantity ${E}_{D}^{\infty }=\langle {{u}^{\infty }}^{2}\rangle -\langle {{B}^{\infty }}^{2}/({\mu }_{0}\rho )\rangle$ denotes the residual energy and is a measure of the dominance of the kinetic $\langle {{u}^{\infty }}^{2}\rangle$ or magnetic $\langle {{B}^{\infty }}^{2}/({\mu }_{0}\rho )\rangle$ energy density.

In deriving a transport formalism, we make the following assumptions:

• 1.  
  The nonlinear terms ${{\boldsymbol{z}}}^{\infty \mp }\cdot {\rm{\nabla }}{{\boldsymbol{z}}}^{\infty \pm }$ in (63) are modeled as

  $$\begin{eqnarray}&&{{NL}}_{\pm }\equiv -{{\boldsymbol{z}}}^{\infty \pm }\displaystyle \frac{{\langle {{{\boldsymbol{z}}}^{\infty \mp }}^{2}\rangle }^{1/2}}{{\lambda }_{\perp }^{\pm }},\end{eqnarray} \tag{ 64 }$$

  and are related to the Kolmogorov phenomenology as described in Section 3.
• 2.  
  To close moments such as $\langle ({{\boldsymbol{z}}}^{\infty -}\cdot {\rm{\nabla }}{\boldsymbol{U}})\cdot {{\boldsymbol{z}}}^{\infty +}\rangle$, we assume that the turbulence is weakly 2D (HZ10) and approximately isotropic in the 2D plane. As discussed in Batchelor (1953), one can express a covariance ${Q}_{{ij}}(0)=C{\delta }_{{ij}}$ for some constant C, allowing us to write

  $$\begin{eqnarray}&&\langle ({{\boldsymbol{z}}}^{\infty \mp }\cdot {\rm{\nabla }}{\boldsymbol{U}})\cdot {{\boldsymbol{z}}}^{\infty \pm }\rangle =\langle {z}_{i}^{\infty -}{z}_{j}^{\infty +}\rangle \displaystyle \frac{\partial {U}_{j}}{\partial {x}_{i}}\\ &&\quad =a\langle {{\boldsymbol{z}}}^{\infty -}\cdot {{\boldsymbol{z}}}^{\infty +}\rangle {\rm{\nabla }}\cdot {\boldsymbol{U}}={{aE}}_{D}^{\infty }{\rm{\nabla }}\cdot {\boldsymbol{U}}.\end{eqnarray} \tag{ 65 }$$

  We retain the parameter a for completeness, although we use $a=1/2$ for 2D turbulence (Zank et al. 2012a).
• 3.  
  Define an orthonormal vector $\hat{{\boldsymbol{n}}}$ orthogonal to the local large-scale mean magnetic field ${{\boldsymbol{B}}}_{0}$ (i.e., $\hat{{\boldsymbol{n}}}\cdot {{\boldsymbol{B}}}_{0}=0$). We approximate the moments

  $$\begin{eqnarray}&&\langle {{z}^{\infty \pm }}^{2}{{\boldsymbol{z}}}^{\infty \pm }\cdot {\rm{\nabla }}\rho \rangle \simeq {\langle {{z}^{\infty \pm }}^{2}\rangle }^{3/2}\hat{{\boldsymbol{n}}}\cdot {\rm{\nabla }}\rho ;\end{eqnarray} \tag{ 66 }$$

  $$\begin{eqnarray}&&\langle {{\boldsymbol{z}}}^{\infty \pm }\cdot {{\boldsymbol{z}}}^{\infty \mp }{{\boldsymbol{z}}}^{\infty \pm }\cdot {\rm{\nabla }}\rho \rangle \simeq {E}_{D}^{\infty }{\langle {{z}^{\infty \pm }}^{2}\rangle }^{1/2}\hat{{\boldsymbol{n}}}\cdot {\rm{\nabla }}\rho ;\end{eqnarray} \tag{ 67 }$$

  Hence, if $\rho ={\rho }_{\mathrm{SW}}(r)$ and ${B}_{0}={B}_{\mathrm{SW}}\hat{{\boldsymbol{r}}}$, $\hat{{\boldsymbol{n}}}\cdot {\rm{\nabla }}{\rho }_{\mathrm{SW}}=0$. Since ${{\boldsymbol{B}}}_{\mathrm{SW}}$ has the form of the Parker spiral, the term $\hat{{\boldsymbol{n}}}\cdot {\rm{\nabla }}{\rho }_{\mathrm{SW}}$ will become increasingly important with increasing heliocentric distance.

Subject to assumptions (1)–(3), the turbulence transport equations for the 2D core equations of NI MHD in the $\beta \sim 1$ (and $\ll 1$) limit can be derived as

$$\begin{eqnarray}&&\displaystyle \frac{\partial \langle {{z}^{\infty \pm }}^{2}\rangle }{\partial t}+{\boldsymbol{U}}\cdot {\rm{\nabla }}\langle {{z}^{\infty \pm }}^{2}\rangle +\displaystyle \frac{1}{2}\langle {{z}^{\infty \pm }}^{2}\rangle {\rm{\nabla }}\cdot {\boldsymbol{U}}\\ &&\qquad +\left(2a-\displaystyle \frac{1}{2}\right){E}_{D}^{\infty }{\rm{\nabla }}\cdot {\boldsymbol{U}}\\ &&\qquad -\displaystyle \frac{1}{2}(\langle {{z}^{\infty \pm }}^{2}\rangle -{E}_{D}^{\infty }){\langle {{z}^{\infty \pm }}^{2}\rangle }^{1/2}\displaystyle \frac{1}{\rho }\hat{{\boldsymbol{n}}}\cdot {\rm{\nabla }}\rho \\ &&\quad =-2\displaystyle \frac{\langle {{z}^{\infty \pm }}^{2}\rangle {\langle {{z}^{\infty \mp }}^{2}\rangle }^{1/2}}{{\lambda }_{\perp }^{\pm }};\end{eqnarray} \tag{ 68 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {E}_{D}^{\infty }}{\partial t}+{\boldsymbol{U}}\cdot {\rm{\nabla }}{E}_{D}^{\infty }+\displaystyle \frac{1}{2}{E}_{D}^{\infty }{\rm{\nabla }}\cdot {\boldsymbol{U}}+\left(2a-\displaystyle \frac{1}{2}\right){E}_{T}^{\infty }{\rm{\nabla }}\cdot {\boldsymbol{U}}\\ &&\qquad -\displaystyle \frac{1}{4}({E}_{D}^{\infty }-{\langle {{z}^{\infty \pm }}^{2}\rangle }^{1/2}{\langle {{z}^{\infty \mp }}^{2}\rangle }^{1/2})\\ &&\qquad \times ({\langle {{z}^{\infty +}}^{2}\rangle }^{1/2}+{\langle {{z}^{\infty -}}^{2}\rangle }^{1/2})\displaystyle \frac{1}{\rho }\hat{{\boldsymbol{n}}}\cdot {\rm{\nabla }}\rho \\ &&\quad =-{E}_{D}^{\infty }\left(\displaystyle \frac{{\langle {{z}^{\infty +}}^{2}\rangle }^{1/2}}{{\lambda }_{\perp }^{+}}+\displaystyle \frac{{\langle {{z}^{\infty -}}^{2}\rangle }^{1/2}}{{\lambda }_{\perp }^{-}}\right);\end{eqnarray} \tag{ 69 }$$

$$\begin{eqnarray}&&{E}_{T}^{\infty }\equiv \displaystyle \frac{\langle {{z}^{\infty +}}^{2}\rangle +\langle {{z}^{\infty -}}^{2}\rangle }{2}.\end{eqnarray} \tag{ 70 }$$

The first three terms of each equation have a "WKB"-like structure and resemble a "pressure" equation with "adiabatic index" 1/2. The forward and backward Elsässer energies are coupled through the large-scale flow variables ${\boldsymbol{U}}$ and ρ, mediated via the residual energy. The nonlinear terms are expressed though dissipation terms that couple forward and backward ${{\boldsymbol{z}}}^{\infty \pm }$ modes, and the rate is determined in part by the correlation length scales ${\lambda }_{\perp }^{\pm }$. The dissipation terms for ${E}_{D}^{\infty }$ were discussed by Dosch et al. (2013), who concluded that the form used originally in Zank et al. (2012a) could lead to unphysical solutions, and they recommended that the form employed in (69) be used instead.

To derive equations describing the correlation lengths ${\lambda }_{\perp }^{\pm }$, we apply the same procedure used in Zank et al. (2012a) to Equation (63). By averaging over small scales, and letting $r=| {\boldsymbol{r}}|$ be the spatial lag between fluctuations, we can define the correlation lengths ${\lambda }_{\perp }^{\pm }$ by

$$\begin{eqnarray}{L}_{\infty }^{\pm } & \equiv & \displaystyle \int \langle {{\boldsymbol{z}}}^{\infty \pm }\cdot {{\boldsymbol{z}}}^{\infty \pm ^{\prime} }\rangle {dr}\equiv \langle {{z}^{\infty \pm }}^{2}\rangle {\lambda }_{\perp }^{\pm };\\ {L}_{D}^{\infty } & \equiv & \displaystyle \int \langle {{\boldsymbol{z}}}^{\infty +}\cdot {{\boldsymbol{z}}}^{\infty -^{\prime} }+{{\boldsymbol{z}}}^{\infty +^{\prime} }\cdot {{\boldsymbol{z}}}^{\infty -}\rangle {dr}\equiv {E}_{D}^{\infty }{\lambda }_{D}^{\infty },\end{eqnarray} \tag{ 71 }$$

where ${{\boldsymbol{z}}}^{\infty \pm ^{\prime} }\equiv {{\boldsymbol{z}}}^{\infty \pm }({\boldsymbol{x}}+{\boldsymbol{r}})$ denotes the lagged Elsässer variable at a location ${\boldsymbol{r}}$ from ${\boldsymbol{x}}$. ${L}_{\infty }^{\pm }$ and ${L}_{D}^{\infty }$ then simply denote the area under the ensemble-averaged two-point variance, which is then related to ${\lambda }_{\perp }^{\infty }$ and ${\lambda }_{D}^{\infty }$ through the definitions in (71). By taking appropriate moments of (63), exploiting assumptions (1)–(3) given previously, and using Equation (C.6) of Appendix C in Zank et al. (2012a), we obtain

$$\begin{eqnarray}&&\displaystyle \frac{\partial {L}_{\infty }^{\pm }}{\partial t}+{\boldsymbol{U}}\cdot {\rm{\nabla }}{L}_{\infty }^{\pm }+\displaystyle \frac{1}{2}{L}_{\infty }^{\pm }{\rm{\nabla }}\cdot {\boldsymbol{U}}+\left(a-\displaystyle \frac{1}{4}\right){L}_{D}^{\infty }{\rm{\nabla }}\cdot {\boldsymbol{U}}\\ &&\quad +\displaystyle \frac{1}{4}{\langle {{z}^{\infty \pm }}^{2}\rangle }^{1/2}({L}_{D}^{\infty }-2{L}_{\infty }^{\pm })\displaystyle \frac{1}{\rho }\hat{{\boldsymbol{n}}}\cdot {\rm{\nabla }}\rho =0;\end{eqnarray} \tag{ 72 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {L}_{D}^{\infty }}{\partial t}+{\boldsymbol{U}}\cdot {\rm{\nabla }}{L}_{D}^{\infty }+\displaystyle \frac{1}{2}{L}_{D}^{\infty }{\rm{\nabla }}\cdot {\boldsymbol{U}}\\ &&\quad +\left(2a-\displaystyle \frac{1}{2}\right)({L}_{\infty }^{+}+{L}_{\infty }^{-}){\rm{\nabla }}\cdot {\boldsymbol{U}}\\ &&\quad -\displaystyle \frac{1}{4}({L}_{D}^{\infty }({\langle {{z}^{\infty +}}^{2}\rangle }^{1/2}+{\langle {{z}^{\infty -}}^{2}\rangle }^{1/2})-2{L}_{\infty }^{+}{\langle {{z}^{\infty -}}^{2}\rangle }^{1/2}\\ &&\quad -2{L}_{\infty }^{-}{\langle {{z}^{\infty +}}^{2}\rangle }^{1/2})\displaystyle \frac{1}{\rho }\hat{{\boldsymbol{n}}}\cdot {\rm{\nabla }}\rho =0.\end{eqnarray} \tag{ 73 }$$

The system of transport equations governing the evolution of 2D incompressible turbulence in the $\beta \sim 1$ limit is now complete and given by Equations (68)–(73). These equations hold for arbitrary inhomogeneous flows, including sub- and super-Alfvénic regimes, provided $\beta \sim 1$ or less. From this system of equations, one can obtain the evolution of the 2D total energy ${E}_{T}^{\infty }$, the cross-helicity ${E}_{C}^{\infty }$, the Alfvén ratio ${r}_{A}^{\infty }$, the kinetic energy density $\langle {{u}^{\infty }}^{2}\rangle$, and the magnetic energy density $\langle {{B}^{\infty }}^{2}/{\mu }_{0}\rho \rangle$ in the fluctuations, via the definitions

$$\begin{eqnarray}{E}_{T}^{\infty } & \equiv & \displaystyle \frac{\langle {{z}^{\infty +}}^{2}\rangle +\langle {{z}^{\infty -}}^{2}\rangle }{2}=\langle {{u}^{\infty }}^{2}\rangle +\langle {{B}^{\infty }}^{2}/{\mu }_{0}\rho \rangle ;\\ \quad {E}_{C}^{\infty } & \equiv & \displaystyle \frac{\langle {{z}^{\infty +}}^{2}\rangle -\langle {{z}^{\infty -}}^{2}\rangle }{2}=2\langle {{\boldsymbol{u}}}^{\infty }\cdot {{\boldsymbol{B}}}^{\infty }/\sqrt{{\mu }_{0}\rho }\rangle ;\\ {r}_{A}^{\infty } & \equiv & \displaystyle \frac{\langle {{u}^{\infty }}^{2}\rangle }{\langle {{B}^{\infty }}^{2}/{\mu }_{0}\rho \rangle }=\displaystyle \frac{{E}_{T}^{\infty }+{E}_{D}^{\infty }}{{E}_{T}^{\infty }-{E}_{D}^{\infty }};\\ \langle {{u}^{\infty }}^{2}\rangle & = & \displaystyle \frac{{E}_{T}^{\infty }+{E}_{D}^{\infty }}{2};\quad \left\langle \displaystyle \frac{{{B}^{\infty }}^{2}}{{\mu }_{0}\rho }\right\rangle =\displaystyle \frac{{E}_{T}^{\infty }-{E}_{D}^{\infty }}{2}.\end{eqnarray} \tag{ 74 }$$

It is useful to relate the three correlation lengths ${\lambda }_{\perp }^{\pm }$ and ${\lambda }_{D}^{\infty }$ to the velocity and magnetic field fluctuation variance correlation lengths ${{\ell }}_{u}^{\infty }$, ${{\ell }}_{b}^{\infty }$, and ${{\ell }}_{{ub}}^{\infty }$. The three latter correlation lengths are defined analogously to the Elsässer variable correlation lengths—that is,

$$\begin{eqnarray*}{{\ell }}_{u}^{\infty }\langle {{u}^{\infty }}^{2}\rangle & \equiv & \displaystyle \int \langle {{\boldsymbol{u}}}^{\infty }\cdot {{{\boldsymbol{u}}}^{\infty }}^{\prime }\rangle {dr};\\ {{\ell }}_{b}^{\infty }\langle {{B}^{\infty }}^{2}\rangle & \equiv & \displaystyle \int \langle {{\boldsymbol{B}}}^{\infty }\cdot {{{\boldsymbol{B}}}^{\infty }}^{\prime }\rangle {dr};\\ {{\ell }}_{{ub}}^{\infty }\langle {{\boldsymbol{u}}}^{\infty }\cdot {{\boldsymbol{B}}}^{\infty }\rangle & \equiv & \displaystyle \frac{1}{2}\displaystyle \int \langle {{\boldsymbol{u}}}^{\infty }\cdot {{{\boldsymbol{B}}}^{\infty }}^{\prime }+{{{\boldsymbol{u}}}^{\infty }}^{\prime }\cdot {{\boldsymbol{B}}}^{\infty }\rangle {dr}.\end{eqnarray*}$$

Following Dosch et al. (2013), it is straightforward to show that

$$\begin{eqnarray}&&{{\ell }}_{u}^{\infty }=\displaystyle \frac{({E}_{T}^{\infty }+{E}_{C}^{\infty }){\lambda }_{\perp }^{+}+({E}_{T}^{\infty }-{E}_{C}^{\infty }){\lambda }_{\perp }^{-}+{E}_{D}^{\infty }{\lambda }_{D}^{\infty }}{2({E}_{T}^{\infty }+{E}_{D}^{\infty })};\end{eqnarray} \tag{ 75 }$$

$$\begin{eqnarray}&&{{\ell }}_{b}^{\infty }=\displaystyle \frac{({E}_{T}^{\infty }+{E}_{C}^{\infty }){\lambda }_{\perp }^{+}+({E}_{T}^{\infty }-{E}_{C}^{\infty }){\lambda }_{\perp }^{-}-{E}_{D}^{\infty }{\lambda }_{D}^{\infty }}{2({E}_{T}^{\infty }-{E}_{D}^{\infty })};\end{eqnarray} \tag{ 76 }$$

$$\begin{eqnarray}&&{{\ell }}_{{ub}}^{\infty }=\displaystyle \frac{({E}_{T}^{\infty }+{E}_{C}^{\infty }){\lambda }_{\perp }^{+}-({E}_{T}^{\infty }-{E}_{C}^{\infty }){\lambda }_{\perp }^{-}}{2{E}_{C}^{\infty }}.\end{eqnarray} \tag{ 77 }$$

These correlation lengths are more easily measured in the solar wind than the correlation lengths for the Elsässer variables, and we need ${{\ell }}_{u}^{\infty }$ below when investigating the transport of the density fluctuation variance.

In concluding this subsection, we use the NI transport Equation (12) to derive an equation describing the evolution of the density fluctuation variance $\langle {{\rho }^{\infty }}^{2}\rangle$. This closely parallels the analysis of Zank et al. (2012b), except that we now compute $\langle {{u}^{\infty }}^{2}\rangle$ directly from the turbulence transport models provided previously. On making the same assumptions (1)–(3) given previously and using the inverse nonlinear time ${\langle {{u}^{\infty }}^{2}\rangle }^{1/2}/{{\ell }}_{u}$, we obtain

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}\langle {{\rho }^{\infty }}^{2}\rangle +{\boldsymbol{U}}\cdot {\rm{\nabla }}\langle {{\rho }^{\infty }}^{2}\rangle -2\langle {{\rho }^{\infty }}^{2}\rangle {\langle {{u}^{\infty }}^{2}\rangle }^{1/2}\displaystyle \frac{1}{\rho }\hat{{\boldsymbol{n}}}\cdot {\rm{\nabla }}\rho \\ &&\quad +2\langle {{\rho }^{\infty }}^{2}\rangle {\rm{\nabla }}\cdot {\boldsymbol{U}}=-\displaystyle \frac{{\langle {{u}^{\infty }}^{2}\rangle }^{1/2}\langle {{\rho }^{\infty }}^{2}\rangle }{{{\ell }}_{u}}.\end{eqnarray} \tag{ 78 }$$

As discussed in Section 3, the density fluctuations, which are independent of ${{\boldsymbol{B}}}^{\infty }$, are mixed passively by the 2D velocity fluctuations ${{\boldsymbol{u}}}^{\infty }$ while being convected in an inhomogeneous flow. These density fluctuations are not associated with slow mode waves but are purely entropic in character and do not propagate.

An Elsässer formulation of the NI model, Equations (13)–(16), can be derived by introducing the fluctuating NI Elsässer variables

$$\begin{eqnarray}&&{{\boldsymbol{z}}}^{* \pm }\equiv {{\boldsymbol{u}}}_{1}\pm \displaystyle \frac{{{\boldsymbol{B}}}^{* }}{\sqrt{{\mu }_{0}\rho ({\boldsymbol{x}})}}\equiv {{\boldsymbol{u}}}_{1}\pm {{\boldsymbol{v}}}_{A}^{* }.\end{eqnarray} \tag{ 79 }$$

It is a straightforward, if lengthy, exercise to derive a coupled pair of equations in ${{\boldsymbol{z}}}^{* \pm }$ from (13) to (16). In our derivation, we did not restrict our attention to incompressible NI corrections from the outset, but instead included both compressible and incompressible terms. This yields a source term of the form

$$\begin{eqnarray*}&&\displaystyle \frac{{\rho }^{\infty }}{\rho ({\boldsymbol{x}})}\left(\displaystyle \frac{\partial {{\boldsymbol{u}}}^{\infty }}{\partial t}+({\boldsymbol{U}}+{{\boldsymbol{u}}}^{\infty })\cdot {\rm{\nabla }}{{\boldsymbol{u}}}^{\infty }\right)\equiv \displaystyle \frac{{\rho }^{\infty }}{\rho ({\boldsymbol{x}})}\displaystyle \frac{D{{\boldsymbol{u}}}^{\infty }}{{Dt}},\end{eqnarray*}$$

for example. The source term can be understood by assuming local homogeneity and neglecting terms that couple to the core 2D equations. If we then follow our discussion in Section 3 about the NI wave modes, we find that

$$\begin{eqnarray}&&\displaystyle \frac{{\partial }^{2}{{\boldsymbol{u}}}_{1}}{\partial {t}^{2}}-{C}_{s0}^{2}{\rm{\nabla }}({\rm{\nabla }}\cdot {{\boldsymbol{u}}}_{1})-{V}_{A0}^{2}\\ &&\quad \times \left({{\rm{\nabla }}}^{2}{u}_{x1}+\displaystyle \frac{\partial \xi }{\partial y},{{\rm{\nabla }}}^{2}{u}_{y1}-\displaystyle \frac{\partial \xi }{\partial x},0\right)=-\displaystyle \frac{{\partial }^{2}{{\boldsymbol{u}}}^{\infty }}{\partial {t}^{2}}.\end{eqnarray} \tag{ 80 }$$

On taking ${{\boldsymbol{B}}}_{0}={B}_{0}\hat{z}$, we obtain Equation (32) again, showing that the core incompressible solutions do not drive vortical modes. From (80), we find that Equation (31) is modified by a source term—that is,

$$\begin{eqnarray}&&\displaystyle \frac{{\partial }^{2}{\rm{\Delta }}}{\partial {t}^{2}}-{C}_{s0}^{2}{{\rm{\nabla }}}^{2}{\rm{\Delta }}-{V}_{A0}^{2}{{\rm{\nabla }}}^{2}({\rm{\Delta }}-{\rm{\Gamma }})=-\displaystyle \frac{{\partial }^{2}}{\partial {t}^{2}}({\rm{\nabla }}\cdot {{\boldsymbol{u}}}^{\infty }).\end{eqnarray} \tag{ 81 }$$

In terms of the NI Mach number ordering, ${\rm{\nabla }}\cdot {{\boldsymbol{u}}}^{\infty }$ is an $O(\varepsilon )$ term and non-zero, and therefore a weak source of compressible (${\rm{\nabla }}\cdot {{\boldsymbol{u}}}_{1}\ne 0$) fluctuations.

In similar fashion, we find from (80) that (33) acquires a source term

$$\begin{eqnarray}&&\displaystyle \frac{{\partial }^{2}{\rm{\Gamma }}}{\partial {t}^{2}}-{C}_{s0}^{2}\displaystyle \frac{{\partial }^{2}{\rm{\Delta }}}{\partial {z}^{2}}=-\displaystyle \frac{{\partial }^{2}{{\rm{\Gamma }}}^{\infty }}{\partial {t}^{2}},\end{eqnarray} \tag{ 82 }$$

where ${{\rm{\Gamma }}}^{\infty }\equiv \partial {u}_{z}^{\infty }/\partial z$ implies $\partial {{\rm{\Gamma }}}^{\infty }/\partial t=-(1/\rho ){\partial }^{2}{P}^{\infty }/\partial {z}^{2}$ is a source term for Γ. Thus incompressible variations in ${{\boldsymbol{u}}}^{\infty }$ and ${P}^{\infty }$ drive compressible fluctuations, which is essentially the Lighthill mechanism for the generation of sound. For incompressible modes, ${\rm{\Delta }}=0$, and vortical modes are not driven by the core 2D ${{\boldsymbol{u}}}^{\infty }$ modes. Accordingly, if we consider incompressible modes ${\rm{\nabla }}\cdot {{\boldsymbol{u}}}_{1}=0$ only, we can neglect all the source terms in the Elsässer formulation of the NI Equations (13)–(16). This then yields the transport equations for ${{\boldsymbol{z}}}^{* \pm }$,

$$\begin{eqnarray}&&\displaystyle \frac{\partial {{{\boldsymbol{z}}}^{* }}^{\pm }}{\partial t}+({\boldsymbol{U}}\mp {{\boldsymbol{V}}}_{A0})\cdot {\rm{\nabla }}{{{\boldsymbol{z}}}^{* }}^{\pm }+{{{\boldsymbol{z}}}^{\infty }}^{\mp }\cdot {\rm{\nabla }}{{{\boldsymbol{z}}}^{* }}^{\pm }+{{{\boldsymbol{z}}}^{* }}^{\mp }\\ &&\qquad \cdot {\rm{\nabla }}{{{\boldsymbol{z}}}^{\infty }}^{\pm }+{{\boldsymbol{z}}}^{* \mp }\cdot {\rm{\nabla }}{{\boldsymbol{z}}}^{* \pm }+\displaystyle \frac{1}{4}({{{\boldsymbol{z}}}^{* }}^{\pm }-{{{\boldsymbol{z}}}^{* }}^{\mp }){\rm{\nabla }}\cdot {\boldsymbol{U}}\\ &&\qquad +{{{\boldsymbol{z}}}^{* }}^{\mp }\cdot {\rm{\nabla }}({\boldsymbol{U}}\pm {{\boldsymbol{V}}}_{A0})\pm \displaystyle \frac{1}{2}\displaystyle \frac{1}{\rho }{{{\boldsymbol{z}}}^{* }}^{\mp }\cdot {\rm{\nabla }}{\rho }_{\mathrm{sw}}{{\boldsymbol{V}}}_{A0}\\ &&\qquad -\displaystyle \frac{1}{4}({{{\boldsymbol{z}}}^{* }}^{+}-{{{\boldsymbol{z}}}^{* }}^{-})\displaystyle \frac{1}{\rho }{{\boldsymbol{V}}}_{A0}\cdot {\rm{\nabla }}\rho \\ &&\qquad \pm \displaystyle \frac{1}{4}({{{\boldsymbol{z}}}^{\infty }}^{+}-{{{\boldsymbol{z}}}^{\infty }}^{-})\displaystyle \frac{1}{\rho }{{{\boldsymbol{z}}}^{* }}^{\mp }\cdot {\rm{\nabla }}\rho \mp \displaystyle \frac{1}{4}({{{\boldsymbol{z}}}^{* }}^{+}-{{{\boldsymbol{z}}}^{* }}^{-})\displaystyle \frac{1}{\rho }{{{\boldsymbol{z}}}^{\infty }}^{\pm }\\ &&\qquad \cdot {\rm{\nabla }}\rho \mp \displaystyle \frac{1}{2}({{{\boldsymbol{z}}}^{* }}^{+}-{{{\boldsymbol{z}}}^{* }}^{-})\cdot {\rm{\nabla }}{\boldsymbol{U}}\\ &&\quad =-\displaystyle \frac{1}{\rho }{\rm{\nabla }}\left({P}^{* }+\displaystyle \frac{1}{{\mu }_{0}}{{\boldsymbol{B}}}^{* }\cdot {\boldsymbol{B}}+\displaystyle \frac{1}{{\mu }_{0}}{{\boldsymbol{B}}}^{* }\cdot {{\boldsymbol{B}}}^{\infty }\right).\end{eqnarray} \tag{ 83 }$$

As before, the higher-order NI Elsässer variables are coupled with the leading-order core ${{\boldsymbol{z}}}^{\infty \pm }$ variables in a passive scalar sense, so that the primary dissipative terms for ${{\boldsymbol{z}}}^{* \pm }$ are due to mixing rather than nonlinearity, which enters at the next order. However, for the reasons discussed in Section 3 for the homogeneous case, we include the higher-order nonlinear terms ${{\boldsymbol{z}}}^{* \mp }\cdot {\rm{\nabla }}{{\boldsymbol{z}}}^{* \pm }$ to ensure the richest equation for ${{\boldsymbol{z}}}^{* \pm }$. As before, the remaining higher-order terms are higher-order analogues of terms already present at a lower order. Like the core 2D transport equations, the forward and backward modes are coupled through gradients in the large-scale inhomogeneous background plasma fields ${\boldsymbol{U}}$, ${{\boldsymbol{B}}}_{0}$, and ρ. Unlike the 2D case, the Alfvén velocity ${{\boldsymbol{V}}}_{A0}$ now enters the transport equation for the ${{\boldsymbol{z}}}^{* \pm }$, including the Alfvén critical point ${\boldsymbol{U}}=\pm {{\boldsymbol{V}}}_{A0}$.

We can take ensemble averages of (83) to derive a turbulence transport model for the energies $\langle {{z}^{* \pm }}^{2}\rangle$ and the residual energy ${E}_{D}^{* }$. The derivation of the transport equations for $\langle {{z}^{* \pm }}^{2}\rangle$ utilizes the following assumptions:

• 1.  
  The nonlinear terms ${{\boldsymbol{z}}}^{* \mp }\cdot {\rm{\nabla }}{{\boldsymbol{z}}}^{* \pm }$ are modeled as

  $$\begin{eqnarray*}&&{NL}{* }_{\pm }\equiv -{{\boldsymbol{z}}}^{* \pm }\displaystyle \frac{{\langle {z}^{* {\mp }^{2}}\rangle }^{1/2}}{{\lambda }_{* }^{\pm }}.\end{eqnarray*}$$
• 2.  
  Since we assume that the minority component is slab with a particular direction ${\boldsymbol{S}}$ defined by the mean magnetic field ${{\boldsymbol{B}}}_{0}$, we can express a covariance with zero lag as ${Q}_{{ij}}(0)=b({\delta }_{{ij}}-{S}_{i}{S}_{j})$ (Batchelor 1953; see the Appendix in Zank et al. 2012a). This yields simplifications such as $\langle {z}_{i}^{* +}{z}_{j}^{* -}\rangle \simeq b\langle {{\boldsymbol{z}}}^{* +}\cdot {{\boldsymbol{z}}}^{* -}\rangle ({\delta }_{{ij}}-{S}_{i}{S}_{j})$, for example.

We obtain

$$\begin{eqnarray}&&\displaystyle \frac{\partial \langle {{{z}^{* }}^{\pm }}^{2}\rangle }{\partial t}+({\boldsymbol{U}}\mp {{\boldsymbol{V}}}_{A0})\cdot {\rm{\nabla }}\langle {{{z}^{* }}^{\pm }}^{2}\rangle +\displaystyle \frac{1}{2}\langle {{{z}^{* }}^{\pm }}^{2}\rangle {\rm{\nabla }}\cdot {\boldsymbol{U}}\\ &&\qquad +\left(2b-\displaystyle \frac{1}{2}\right){\rm{\nabla }}\cdot {\boldsymbol{U}}{E}_{D}^{* }-2{{bE}}_{D}^{* }{S}_{i}{S}_{j}\displaystyle \frac{\partial {U}_{i}}{\partial {x}_{j}}\\ &&\qquad \,\pm {{bE}}_{D}^{* }\left(2{\rm{\nabla }}\cdot {{\boldsymbol{V}}}_{A0}-2{S}_{i}{S}_{j}\displaystyle \frac{\partial {V}_{{Aj}}}{\partial {x}_{j}}\right.\\ &&\qquad \left.+\displaystyle \frac{1}{\rho }{{\boldsymbol{V}}}_{A0}\cdot {\rm{\nabla }}\rho -{S}_{i}{S}_{j}{V}_{{Ai}}\displaystyle \frac{1}{\rho }\displaystyle \frac{\partial \rho }{\partial {x}_{j}}\right)\mp \displaystyle \frac{1}{2}(\langle {{{z}^{* }}^{\pm }}^{2}\rangle -{E}_{D}^{* })\\ &&\qquad \times \left(\displaystyle \frac{1}{\rho }{{\boldsymbol{V}}}_{A0}\cdot {\rm{\nabla }}\rho \pm \langle {{{z}^{\infty }}^{\pm }}^{2}{\rangle }^{1/2}\displaystyle \frac{1}{\rho }\hat{n}\right.\\ &&\qquad \cdot \left.{\rm{\nabla }}\rho \pm 2b\left({\rm{\nabla }}\cdot {\boldsymbol{U}}-{S}_{i}{S}_{j}\displaystyle \frac{\partial {U}_{i}}{\partial {x}_{j}}\right)\right)\\ &&\quad =-2\displaystyle \frac{{\langle {{{z}^{\infty }}^{\mp }}^{2}\rangle }^{1/2}\langle {{z}^{* \pm }}^{2}\rangle }{{\lambda }_{\perp }^{\pm }}-2\displaystyle \frac{{\langle {{z}^{* \mp }}^{2}\rangle }^{1/2}\langle {{z}^{* \pm }}^{2}\rangle }{{\lambda }_{* }^{\pm }}\end{eqnarray} \tag{ 84 }$$

after again assuming that $\langle {{\boldsymbol{z}}}^{* \pm }\cdot {{\boldsymbol{z}}}^{\infty \pm }\rangle \simeq 0$. Notice that the inverse dissipation time is the sum of the inverse mixing time and the inverse nonlinear time, in accordance with the discussion in Section 3. In deriving the transport equation for the residual energy ${E}_{D}^{* }\equiv \langle {{\boldsymbol{z}}}^{* +}\cdot {{\boldsymbol{z}}}^{* -}\rangle$, we have to be careful with the closure that involves the terms

$$\begin{eqnarray}&&\langle {{{\boldsymbol{z}}}^{* }}^{+}\cdot ({{\boldsymbol{V}}}_{A0}\cdot {\rm{\nabla }}{{{\boldsymbol{z}}}^{* }}^{-})\rangle -\langle {{{\boldsymbol{z}}}^{* }}^{-}\cdot ({{\boldsymbol{V}}}_{A0}\cdot {\rm{\nabla }}{{{\boldsymbol{z}}}^{* }}^{+})\rangle \\ &&\quad ={V}_{A0}\left(\langle {{z}_{i}^{* }}^{+}\displaystyle \frac{\partial {{z}_{i}^{* }}^{-}}{\partial z}\rangle -\langle {{z}_{i}^{* }}^{-}\displaystyle \frac{\partial {{z}_{i}^{* }}^{+}}{\partial z}\rangle \right).\end{eqnarray} \tag{ 85 }$$

The term ${{\boldsymbol{V}}}_{A0}\cdot {\rm{\nabla }}$ defines derivatives along the mean magnetic field ${{\boldsymbol{B}}}_{0}$ direction. Unlike ${{\boldsymbol{z}}}^{\infty \pm }$, we cannot assume that ${{\boldsymbol{z}}}^{* \pm }$ is approximately 2D perpendicular to ${{\boldsymbol{B}}}_{0}$. Instead, ${{\boldsymbol{z}}}^{* \pm }$ could correspond to Alfvén waves propagating along ${{\boldsymbol{B}}}_{0}$. Derivatives along $\hat{z}$, say (assuming ${{\boldsymbol{B}}}_{0}={B}_{0}\hat{z}$) $\partial {{\boldsymbol{z}}}^{* \pm }/\partial z$, can therefore have high-frequency/short wavelength variation, making such terms essentially secular, and the terms cannot therefore be eliminated by averaging. Fortunately, if the cross-helicity E_C = 0—that is, $\langle {{z}^{* +}}^{2}\rangle =\langle {{z}^{* -}}^{2}\rangle$—then (85) is $\simeq 0$. If wave propagation is unidirectional, then either ${{\boldsymbol{z}}}^{* +}=0$ or ${{\boldsymbol{z}}}^{* -}=0$, and (85) is zero for either case. We therefore assume that the two secular terms cancel—that is, $\langle {{{\boldsymbol{z}}}^{* }}^{+}\cdot ({{\boldsymbol{V}}}_{A0}\cdot {\rm{\nabla }}{{{\boldsymbol{z}}}^{* }}^{-})\rangle \,-\langle {{{\boldsymbol{z}}}^{* }}^{-}\cdot ({{\boldsymbol{V}}}_{A0}\cdot {\rm{\nabla }}{{{\boldsymbol{z}}}^{* }}^{+})\rangle =0$. Some lengthy algebra then yields the transport equation for ${E}_{D}^{* }$,

$$\begin{eqnarray}&&\displaystyle \frac{\partial {E}_{D}^{* }}{\partial t}+{\boldsymbol{U}}\cdot {\rm{\nabla }}{E}_{D}^{* }+\displaystyle \frac{1}{2}{E}_{D}^{* }{\rm{\nabla }}\cdot {\boldsymbol{U}}\\ &&\qquad +\left(3b-\displaystyle \frac{1}{2}\right){E}_{T}^{* }{\rm{\nabla }}\cdot {\boldsymbol{U}}-3{{bE}}_{T}^{* }{S}_{i}{S}_{j}\displaystyle \frac{\partial {U}_{i}}{\partial {x}_{j}}\\ &&\qquad -{{bE}}_{D}^{* }\left({\rm{\nabla }}\cdot {\boldsymbol{U}}-{S}_{i}{S}_{j}\displaystyle \frac{\partial {U}_{i}}{\partial {x}_{j}}\right)-2{{bE}}_{C}^{* }\\ &&\qquad \times \left({\rm{\nabla }}\cdot {{\boldsymbol{V}}}_{A0}-{S}_{i}{S}_{j}\displaystyle \frac{\partial {V}_{{Ai}}}{\partial {x}_{j}}\right)-\displaystyle \frac{1}{2}{E}_{C}^{* }\displaystyle \frac{1}{\rho }{{\boldsymbol{V}}}_{A0}\cdot {\rm{\nabla }}\rho \\ &&\qquad -\displaystyle \frac{1}{\rho }{{bE}}_{C}^{* }\left({{\boldsymbol{V}}}_{A0}\cdot {\rm{\nabla }}\rho -{S}_{i}{S}_{j}\displaystyle \frac{\partial \rho }{\partial {x}_{j}}{V}_{{Ai}}\right)\\ &&\qquad -\displaystyle \frac{1}{4}\displaystyle \frac{1}{\rho }[{E}_{D}^{* }(\langle {{{z}^{\infty }}^{+}}^{2}{\rangle }^{1/2}+\langle {{{z}^{\infty }}^{-}}^{2}{\rangle }^{1/2})\\ &&\qquad -\langle {{{z}^{* }}^{-}}^{2}\rangle \langle {{{z}^{\infty }}^{+}}^{2}{\rangle }^{1/2}-\langle {{{z}^{* }}^{+}}^{2}\rangle \langle {{{z}^{\infty }}^{-}}^{2}{\rangle }^{1/2}]\hat{n}\cdot {\rm{\nabla }}\rho \\ &&\quad =-{E}_{D}^{* }\left(\displaystyle \frac{\langle {{{z}^{\infty }}^{-}}^{2}{\rangle }^{1/2}}{{\lambda }_{\perp }^{+}}+\displaystyle \frac{\langle {{{z}^{\infty }}^{+}}^{2}{\rangle }^{1/2}}{{\lambda }_{\perp }^{-}}\right)\\ &&\qquad -{E}_{D}^{* }\left(\displaystyle \frac{\langle {{z}^{* +}}^{2}\rangle }{{\lambda }_{* }^{+}}+\displaystyle \frac{\langle {{z}^{* -}}^{2}\rangle }{{\lambda }_{* }^{-}}\right).\end{eqnarray} \tag{ 86 }$$

Finally, we can derive transport equations for ${L}_{* }^{\pm }$ and ${L}_{D}^{* }$, using the definitions,

$$\begin{eqnarray*}&&{L}_{* }^{\pm } = \displaystyle \int {\langle {{{\boldsymbol{z}}}^{*\pm }}\cdot {{{\boldsymbol{z}}}^{* \pm \prime}} \rangle {dr}\equiv \langle {{z}^{* \pm }}\cdot {{z}^{* \pm }}\rangle {\lambda }_{* }^{\pm }};\\ &&{L}_{D}^{* } = \displaystyle \int {\langle {{{\boldsymbol{z}}}^{* \pm }}\cdot {{{\boldsymbol{z}}}^{* \mp \prime}} +{{{\boldsymbol{z}}}^{* \pm \prime}} \cdot {{{\boldsymbol{z}}}^{* \mp }}\rangle {dr}\equiv {E}_{D}^{* }{\lambda }_{D}^{* }},\end{eqnarray*}$$

which define the respective correlation lengths. The transport equations can be expressed as

$$\begin{eqnarray}&&\displaystyle \frac{\partial {L}_{* }^{\pm }}{\partial t}+({\boldsymbol{U}}\mp {{\boldsymbol{V}}}_{A0})\cdot {\rm{\nabla }}{L}_{* }^{\pm }+\displaystyle \frac{1}{2}{L}_{* }^{\pm }{\rm{\nabla }}\cdot {\boldsymbol{U}}\\ &&\quad +\left(b-\displaystyle \frac{1}{4}\right){L}_{D}^{* }{\rm{\nabla }}\cdot {\boldsymbol{U}}-{{bL}}_{D}^{* }{S}_{i}{S}_{j}\displaystyle \frac{\partial {U}_{i}}{\partial {x}_{j}}\\ &&\quad \pm {{bL}}_{D}^{* }\left({\rm{\nabla }}\cdot {{\boldsymbol{V}}}_{A0}+\displaystyle \frac{1}{2}\displaystyle \frac{1}{{\rho }_{\mathrm{sw}}}{{\boldsymbol{V}}}_{A0}\cdot {\rm{\nabla }}\rho \right.\\ &&\quad \left.-{S}_{i}{S}_{j}\displaystyle \frac{\partial {V}_{{Ai}}}{\partial {x}_{j}}-\displaystyle \frac{1}{2}{S}_{i}{S}_{j}{V}_{{Ai}}\displaystyle \frac{1}{\rho }\displaystyle \frac{\partial {\rho }_{\mathrm{sw}}}{\partial {x}_{j}}\right)\\ &&\quad \mp \displaystyle \frac{1}{2}\left({L}_{* }^{\pm }-\displaystyle \frac{{L}_{D}^{* }}{2}\right)\left[\displaystyle \frac{1}{\rho }{{\boldsymbol{V}}}_{A0}\cdot {\rm{\nabla }}\rho \pm \displaystyle \frac{1}{\rho }\langle {{{z}^{\infty }}^{\pm }}^{2}{\rangle }^{1/2}\hat{n}\right.\\ &&\quad \left.\cdot {\rm{\nabla }}\rho \pm 2b\left({\rm{\nabla }}\cdot {\boldsymbol{U}}-{S}_{i}{S}_{j}\displaystyle \frac{\partial {U}_{i}}{\partial {x}_{j}}\right)\right]=0;\end{eqnarray} \tag{ 87 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {L}_{D}^{* }}{\partial t}+{\boldsymbol{U}}\cdot {\rm{\nabla }}{L}_{D}^{* }-\left(b-\displaystyle \frac{1}{2}\right){L}_{D}^{* }{\rm{\nabla }}\cdot {\boldsymbol{U}}\\ &&\quad +\left(3b-\displaystyle \frac{1}{2}\right)({L}_{* }^{-}+{L}_{* }^{+}){\rm{\nabla }}\cdot {\boldsymbol{U}}\\ &&\quad -3b({L}_{* }^{-}+{L}_{* }^{+}){S}_{i}{S}_{j}\displaystyle \frac{\partial {U}_{i}}{\partial {x}_{j}}+{{bL}}_{D}^{* }{S}_{i}{S}_{j}\displaystyle \frac{\partial {U}_{i}}{\partial {x}_{j}}\\ &&\quad -2b({L}_{* }^{+}-{L}_{* }^{-})\left({\rm{\nabla }}\cdot {{\boldsymbol{V}}}_{A0}-{S}_{i}{S}_{j}\displaystyle \frac{\partial {V}_{{Ai}}}{\partial {x}_{j}}\right)\\ &&\quad -\displaystyle \frac{1}{2}({L}_{* }^{+}-{L}_{* }^{-})\left[(2b+1)\displaystyle \frac{1}{\rho }{{\boldsymbol{V}}}_{A0}\cdot {\rm{\nabla }}\rho \right)\\ &&\quad \left.-2{{bS}}_{i}{S}_{j}{V}_{{Ai}}\displaystyle \frac{1}{\rho }\displaystyle \frac{\partial \rho }{\partial {x}_{j}}\right]\\ &&\quad +\displaystyle \frac{1}{2}\left[\left({L}_{* }^{-}-\displaystyle \frac{{L}_{D}^{* }}{2}\right)\langle {{{z}^{\infty }}^{+}}^{2}{\rangle }^{1/2}\right.\\ &&\quad \left.+\left({L}_{* }^{+}-\displaystyle \frac{{L}_{D}^{* }}{2}\right)\langle {{{z}^{\infty }}^{-}}^{2}{\rangle }^{1/2}\right]\displaystyle \frac{1}{\rho }\hat{n}\cdot {\rm{\nabla }}\rho =0;\end{eqnarray} \tag{ 88 }$$

In concluding this section, we write down the 2D core transport equations that now include the source terms that correspond to the "dissipative driving" of 2D fluctuations by the nonlinear dissipation of slab turbulence. Equations (68)–(69) are modified to read

$$\begin{eqnarray}&&\displaystyle \frac{\partial \langle {{z}^{\infty \pm }}^{2}\rangle }{\partial t}+{\boldsymbol{U}}\cdot {\rm{\nabla }}\langle {{z}^{\infty \pm }}^{2}\rangle +\displaystyle \frac{1}{2}\langle {{z}^{\infty \pm }}^{2}\rangle {\rm{\nabla }}\cdot {\boldsymbol{U}}\\ &&\qquad +\left(2a-\displaystyle \frac{1}{2}\right){E}_{D}^{\infty }{\rm{\nabla }}\cdot {\boldsymbol{U}}\\ &&\qquad -\displaystyle \frac{1}{2}(\langle {{z}^{\infty \pm }}^{2}\rangle -{E}_{D}^{\infty }){\langle {{z}^{\infty \pm }}^{2}\rangle }^{1/2}\displaystyle \frac{1}{\rho }\hat{{\boldsymbol{n}}}\cdot {\rm{\nabla }}\rho \\ &&\quad =-2\displaystyle \frac{\langle {{z}^{\infty \pm }}^{2}\rangle {\langle {{z}^{\infty \mp }}^{2}\rangle }^{1/2}}{{\lambda }_{\perp }^{\pm }}+2\displaystyle \frac{{\langle {{z}^{* \mp }}^{2}\rangle }^{1/2}\langle {z}^{* {\pm }^{2}}\rangle }{{\lambda }^{* \pm }};\end{eqnarray} \tag{ 89 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {E}_{D}^{\infty }}{\partial t}+{\boldsymbol{U}}\cdot {\rm{\nabla }}{E}_{D}^{\infty }+\displaystyle \frac{1}{2}{E}_{D}^{\infty }{\rm{\nabla }}\cdot {\boldsymbol{U}}+\left(2a-\displaystyle \frac{1}{2}\right){E}_{T}^{\infty }{\rm{\nabla }}\cdot {\boldsymbol{U}}\\ &&\qquad -\displaystyle \frac{1}{4}({E}_{D}^{\infty }-{\langle {{z}^{\infty \pm }}^{2}\rangle }^{1/2}{\langle {{z}^{\infty \mp }}^{2}\rangle }^{1/2})\\ &&\qquad \times ({\langle {{z}^{\infty +}}^{2}\rangle }^{1/2}+{\langle {{z}^{\infty -}}^{2}\rangle }^{1/2})\displaystyle \frac{1}{\rho }\hat{{\boldsymbol{n}}}\cdot {\rm{\nabla }}\rho \\ &&\quad =-{E}_{D}^{\infty }\left(\displaystyle \frac{{\langle {{z}^{\infty +}}^{2}\rangle }^{1/2}}{{\lambda }_{\perp }^{+}}+\displaystyle \frac{{\langle {{z}^{\infty -}}^{2}\rangle }^{1/2}}{{\lambda }_{\perp }^{-}}\right)\\ &&\qquad +{E}_{D}^{* }\left(\displaystyle \frac{\langle {{z}^{* +}}^{2}\rangle }{{\lambda }^{* +}}+\displaystyle \frac{\langle {{z}^{* -}}^{2}\rangle }{{\lambda }^{* -}}\right).\end{eqnarray} \tag{ 90 }$$

The correlation length Equations (72)–(73) are unchanged. The 2D and slab transport equations are now fully coupled, and a set of 12 transport equations must now be solved.

By way of example, let us apply the transport model (110), (111), (84), (86), (72), (73), (87), and (88) to the supersonic solar wind and consider two simplified spherically symmetric cases (see Adhikari et al. 2015a, 2015b).

6.1. Spherically Symmetric Model Equations

For the inner heliosphere, consider a radial background magnetic field ${{\boldsymbol{B}}}_{0}={B}_{\mathrm{SW}}\hat{{\boldsymbol{r}}}$ and $\hat{{\boldsymbol{n}}}=1/\sqrt{2}(0,1,1)$, ${\boldsymbol{S}}=(1,0,0)$, where ${B}_{\mathrm{SW}}={B}_{0}{({r}_{0}/r)}^{2}$ is given by the radial component of the Parker spiral magnetic field. Assume further a constant radial solar wind velocity ${\boldsymbol{U}}=U\hat{{\boldsymbol{r}}}$, choosing $a=1/2$ and $\rho =\rho (r)={\rho }_{0}{({r}_{0}/r)}^{2}$, where the subscript refers to a reference value at a reference heliocentric distance r₀. The steady-state reduced 2D model is given by the six equations,

$$\begin{eqnarray}&&U\displaystyle \frac{d}{{dr}}\langle {{z}^{\infty \pm }}^{2}\rangle +\displaystyle \frac{U}{r}{E}_{D}^{\infty }+\displaystyle \frac{U}{r}{E}_{T}^{\infty }\\ &&\quad =-2\displaystyle \frac{\langle {{z}^{\infty \pm }}^{2}\rangle {\langle {{z}^{\infty \mp }}^{2}\rangle }^{1/2}}{{\lambda }_{\perp }^{\pm }}+2\displaystyle \frac{{\langle {{z}^{* \mp }}^{2}\rangle }^{1/2}\langle {{z}^{* \pm }}^{2}\rangle }{{\lambda }^{* \pm }};\end{eqnarray} \tag{ 91 }$$

$$\begin{eqnarray}&&U\displaystyle \frac{{{dE}}_{D}^{\infty }}{{dr}}+\displaystyle \frac{U}{r}{E}_{D}^{\infty }+\displaystyle \frac{U}{r}{E}_{T}^{\infty }\\ &&\quad =-{E}_{D}^{\infty }\left(\displaystyle \frac{{\langle {{z}^{\infty -}}^{2}\rangle }^{1/2}}{{\lambda }_{\perp }^{+}}+\displaystyle \frac{{\langle {{z}^{\infty +}}^{2}\rangle }^{1/2}}{{\lambda }_{\perp }^{-}}\right)\\ &&\qquad +{E}_{D}^{* }\left(\displaystyle \frac{\langle {{z}^{* +}}^{2}\rangle }{{\lambda }^{* +}}+\displaystyle \frac{\langle {{z}^{* -}}^{2}\rangle }{{\lambda }^{* -}}\right);\end{eqnarray} \tag{ 92 }$$

$$\begin{eqnarray}&&U\displaystyle \frac{{{dL}}_{\infty }^{\pm }}{{dr}}+\displaystyle \frac{U}{r}{L}_{\infty }^{\pm }+\displaystyle \frac{U}{2r}{L}_{D}^{\infty }=0;\end{eqnarray} \tag{ 93 }$$

$$\begin{eqnarray}&&U\displaystyle \frac{{{dL}}_{D}^{\infty }}{{dr}}+\displaystyle \frac{U}{r}{L}_{D}^{\infty }+\displaystyle \frac{U}{r}({L}_{\infty }^{+}+{L}_{\infty }^{-})=0.\end{eqnarray} \tag{ 94 }$$

The 2D core equations are coupled to the slab turbulence component through the six steady-state spherically symmetric 1D transport equations

$$\begin{eqnarray}&&\left(U\mp {V}_{A0}\displaystyle \frac{{r}_{0}}{r}\right)\displaystyle \frac{d\langle {{{z}^{* }}^{\pm }}^{2}\rangle }{{dr}}-(2b-1)\displaystyle \frac{U}{r}\langle {{{z}^{* }}^{\pm }}^{2}\rangle \\ &&\qquad +2\left(3b-\displaystyle \frac{1}{2}\right)\displaystyle \frac{U}{r}{E}_{D}^{* }\\ &&\qquad \pm (4b-1)\displaystyle \frac{{V}_{A0}}{{r}_{0}}{\left(\displaystyle \frac{{r}_{0}}{r}\right)}^{2}{E}_{D}^{* }\pm \displaystyle \frac{{V}_{A0}}{{r}_{0}}{\left(\displaystyle \frac{{r}_{0}}{r}\right)}^{2}\langle {{{z}^{* }}^{\pm }}^{2}\rangle \\ &&\quad =-2\displaystyle \frac{\langle {{{z}^{* }}^{\pm }}^{2}\rangle \langle {{{z}^{\infty }}^{\mp }}^{2}{\rangle }^{1/2}}{{\lambda }_{* }^{\pm }}-2\displaystyle \frac{\langle {{{z}^{* }}^{\pm }}^{2}\rangle \langle {{{z}^{* }}^{\mp }}^{2}{\rangle }^{1/2}}{{\lambda }_{* }^{\pm }};\end{eqnarray} \tag{ 95 }$$

$$\begin{eqnarray}&&U\displaystyle \frac{{{dE}}_{D}^{* }}{{dr}}-(2b-1)\displaystyle \frac{U}{r}{E}_{D}^{* }+2\left(3b-\displaystyle \frac{1}{2}\right)\displaystyle \frac{U}{r}{E}_{T}^{* }\\ &&\qquad -(4b-1){E}_{C}^{* }\displaystyle \frac{{V}_{A0}}{{r}_{0}}{\left(\displaystyle \frac{{r}_{0}}{r}\right)}^{2}\\ &&\quad =-{E}_{D}^{* }\left(\displaystyle \frac{\langle {{{z}^{\infty }}^{-}}^{2}{\rangle }^{1/2}}{{\lambda }_{* }^{+}}+\displaystyle \frac{\langle {{{z}^{\infty }}^{+}}^{2}{\rangle }^{1/2}}{{\lambda }_{* }^{-}}\right)\\ &&\qquad -{E}_{D}^{* }\left(\displaystyle \frac{\langle {{{z}^{* }}^{-}}^{2}{\rangle }^{1/2}}{{\lambda }_{* }^{+}}+\displaystyle \frac{\langle {{{z}^{* }}^{+}}^{2}{\rangle }^{1/2}}{{\lambda }_{* }^{-}}\right);\quad \end{eqnarray} \tag{ 96 }$$

$$\begin{eqnarray}&&\left(U\mp {V}_{A0}\displaystyle \frac{{r}_{0}}{r}\right)\displaystyle \frac{{{dL}}_{* }^{\pm }}{{dr}}-(2b-1)\displaystyle \frac{U}{r}{L}_{* }^{\pm }\\ &&\quad +2\left(\displaystyle \frac{3b}{2}-\displaystyle \frac{1}{4}\right)\displaystyle \frac{U}{r}{L}_{D}^{* }\pm \displaystyle \frac{{V}_{A0}}{{r}_{0}}{\left(\displaystyle \frac{{r}_{0}}{r}\right)}^{2}{L}_{* }^{\pm }\\ &&\quad \pm \left(2b-\displaystyle \frac{1}{2}\right)\displaystyle \frac{{V}_{A0}}{{r}_{0}}{\left(\displaystyle \frac{{r}_{0}}{r}\right)}^{2}{L}_{D}^{* }=0;\quad \end{eqnarray} \tag{ 97 }$$

$$\begin{eqnarray}&&U\displaystyle \frac{{{dL}}_{D}^{* }}{{dr}}-(2b-1)\displaystyle \frac{U}{r}{L}_{D}^{* }+2\left(3b-\displaystyle \frac{1}{2}\right)({L}_{* }^{+}+{L}_{* }^{-})\displaystyle \frac{U}{r}\\ &&\quad -(4b-1)({L}_{* }^{+}-{L}_{* }^{-})\displaystyle \frac{{V}_{A0}}{{r}_{0}}{\left(\displaystyle \frac{{r}_{0}}{r}\right)}^{2}=0.\quad \end{eqnarray} \tag{ 98 }$$

This system of equations is solved as follows, from 0.3 to 5 au for b = 0.26.

Consider now a model for which $\hat{{\boldsymbol{r}}}\,\perp \,{{\boldsymbol{B}}}_{\mathrm{SW}}$ and $\hat{{\boldsymbol{n}}}=(1,0,0)$. This is certainly appropriate to the more distant heliosphere where the Parker spiral magnetic field is sufficiently wound up so that the azimuthal magnetic field component dominates. The model 2D core equations are now given by

$$\begin{eqnarray}&&U\displaystyle \frac{d}{{dr}}\langle {{z}^{\infty \pm }}^{2}\rangle +\displaystyle \frac{U}{r}\langle {{z}^{\infty \pm }}^{2}\rangle +\displaystyle \frac{U}{r}{E}_{D}^{\infty }\\ &&\qquad +\displaystyle \frac{1}{r}{\langle {{z}^{\infty \pm }}^{2}\rangle }^{1/2}(\langle {{z}^{\infty \pm }}^{2}\rangle -{E}_{D}^{\infty })\\ &&\quad =-2\displaystyle \frac{\langle {{z}^{\infty \pm }}^{2}\rangle {\langle {{z}^{\infty \mp }}^{2}\rangle }^{1/2}}{{\lambda }_{\perp }^{\pm }}+2\displaystyle \frac{{\langle {{z}^{* \mp }}^{2}\rangle }^{1/2}\langle {{z}^{* \pm }}^{2}\rangle }{{\lambda }^{* \pm }};\end{eqnarray} \tag{ 99 }$$

$$\begin{eqnarray}&&U\displaystyle \frac{{{dE}}_{D}^{\infty }}{{dr}}+\displaystyle \frac{U}{r}{E}_{D}^{\infty }+\displaystyle \frac{U}{r}{E}_{T}^{\infty }+\displaystyle \frac{1}{2r}({E}_{D}^{\infty }-{\langle {{z}^{\infty +}}^{2}\rangle }^{1/2}{\langle {{z}^{\infty -}}^{2}\rangle }^{1/2})\\ &&\qquad \times ({\langle {{z}^{\infty +}}^{2}\rangle }^{1/2}+{\langle {{z}^{\infty -}}^{2}\rangle }^{1/2})\\ &&\quad =-{E}_{D}^{\infty }\left(\displaystyle \frac{{\langle {{z}^{\infty -}}^{2}\rangle }^{1/2}}{{\lambda }_{\perp }^{+}}+\displaystyle \frac{{\langle {{z}^{\infty +}}^{2}\rangle }^{1/2}}{{\lambda }_{\perp }^{-}}\right)\\ &&\qquad +{E}_{D}^{* }\left(\displaystyle \frac{\langle {{z}^{* +}}^{2}\rangle }{{\lambda }^{* +}}+\displaystyle \frac{\langle {{z}^{* -}}^{2}\rangle }{{\lambda }^{* -}}\right);\end{eqnarray} \tag{ 100 }$$

$$\begin{eqnarray}&&U\displaystyle \frac{{{dL}}_{\infty }^{\pm }}{{dr}}+\displaystyle \frac{U}{r}{L}_{\infty }^{\pm }+\displaystyle \frac{U}{2r}{L}_{D}^{\infty }\\ &&\quad -\displaystyle \frac{1}{2r}{\langle {{z}^{\infty \pm }}^{2}\rangle }^{1/2}({L}_{D}^{\infty }-2{L}_{\infty }^{\pm })=0;\end{eqnarray} \tag{ 101 }$$

$$\begin{eqnarray}&&U\displaystyle \frac{{{dL}}_{D}^{\infty }}{{dr}}+\displaystyle \frac{U}{r}{L}_{D}^{\infty }+\displaystyle \frac{U}{r}({L}_{\infty }^{+}+{L}_{\infty }^{-})+\displaystyle \frac{1}{2r}({L}_{D}^{\infty }({\langle {{z}^{\infty +}}^{2}\rangle }^{1/2}\\ &&\quad +{\langle {{z}^{\infty -}}^{2}\rangle }^{1/2})-2{L}_{\infty }^{+}{\langle {{z}^{\infty -}}^{2}\rangle }^{1/2}-2{L}_{\infty }^{-}{\langle {{z}^{\infty +}}^{2}\rangle }^{1/2})=0.\end{eqnarray} \tag{ 102 }$$

The steady-state slab turbulence component spherically symmetric equations for a perpendicular magnetic field geometry are given by

$$\begin{eqnarray}&&\left(U\mp {V}_{A0}\displaystyle \frac{{r}_{0}}{r}\right)\displaystyle \frac{d\langle {{{z}^{* }}^{\pm }}^{2}\rangle }{{dr}}-(2b-1)\displaystyle \frac{U}{r}\langle {{{z}^{* }}^{\pm }}^{2}\rangle \\ &&\qquad +2\left(3b-\displaystyle \frac{1}{2}\right)\displaystyle \frac{U}{r}{E}_{D}^{* }\\ &&\qquad \pm (4b-1)\displaystyle \frac{{V}_{A0}}{{r}_{0}}{\left(\displaystyle \frac{{r}_{0}}{r}\right)}^{2}{E}_{D}^{* }\pm \displaystyle \frac{{V}_{A0}}{{r}_{0}}{\left(\displaystyle \frac{{r}_{0}}{r}\right)}^{2}\langle {{{z}^{* }}^{\pm }}^{2}\rangle \\ &&\qquad +\displaystyle \frac{1}{r}(\langle {{{z}^{* }}^{\pm }}^{2}\rangle -{E}_{D}^{* })\langle {{{z}^{\infty }}^{\pm }}^{2}{\rangle }^{1/2}\\ &&\quad =-2\displaystyle \frac{\langle {{{z}^{* }}^{\pm }}^{2}\rangle \langle {{{z}^{\infty }}^{\mp }}^{2}{\rangle }^{1/2}}{{\lambda }_{* }^{\pm }}-2\displaystyle \frac{\langle {{{z}^{* }}^{\pm }}^{2}\rangle \langle {{{z}^{* }}^{\mp }}^{2}{\rangle }^{1/2}}{{\lambda }_{* }^{\pm }};\end{eqnarray} \tag{ 103 }$$

$$\begin{eqnarray}&&U\displaystyle \frac{{{dE}}_{D}^{* }}{{dr}}-(2b-1)\displaystyle \frac{U}{r}{E}_{D}^{* }+2\left(3b-\displaystyle \frac{1}{2}\right)\displaystyle \frac{U}{r}{E}_{T}^{* }\\ &&\qquad -(4b-1){E}_{C}^{* }\displaystyle \frac{{V}_{A0}}{{r}_{0}}{\left(\displaystyle \frac{{r}_{0}}{r}\right)}^{2}\\ &&\qquad +\displaystyle \frac{1}{2r}[{E}_{D}^{* }(\langle {{{z}^{\infty }}^{+}}^{2}{\rangle }^{1/2}+\langle {{{z}^{\infty }}^{-}}^{2}{\rangle }^{1/2})\\ &&\qquad -\langle {{{z}^{* }}^{-}}^{2}\rangle \langle {{{z}^{\infty }}^{+}}^{2}{\rangle }^{1/2}-\langle {{{z}^{* }}^{+}}^{2}\rangle \langle {{{z}^{\infty }}^{-}}^{2}{\rangle }^{1/2}]\\ &&\quad =-{E}_{D}^{* }\left(\displaystyle \frac{\langle {{{z}^{\infty }}^{-}}^{2}{\rangle }^{1/2}}{{\lambda }_{* }^{+}}+\displaystyle \frac{\langle {{{z}^{\infty }}^{+}}^{2}{\rangle }^{1/2}}{{\lambda }_{* }^{-}}\right)\\ &&\qquad -{E}_{D}^{* }\left(\displaystyle \frac{\langle {{{z}^{* }}^{-}}^{2}{\rangle }^{1/2}}{{\lambda }_{* }^{+}}+\displaystyle \frac{\langle {{{z}^{* }}^{+}}^{2}{\rangle }^{1/2}}{{\lambda }_{* }^{-}}\right);\end{eqnarray} \tag{ 104 }$$

$$\begin{eqnarray}&&\left(U\mp {V}_{A0}\displaystyle \frac{{r}_{0}}{r}\right)\displaystyle \frac{{{dL}}_{* }^{\pm }}{{dr}}-(2b-1)\displaystyle \frac{U}{r}{L}_{* }^{\pm }+2\left(\displaystyle \frac{3}{2}b-\displaystyle \frac{1}{4}\right)\displaystyle \frac{U}{r}{L}_{D}^{* }\\ &&\quad \pm \left(2b-\displaystyle \frac{1}{2}\right)\displaystyle \frac{{V}_{A0}}{{r}_{0}}{\left(\displaystyle \frac{{r}_{0}}{r}\right)}^{2}{L}_{D}^{* }\pm \displaystyle \frac{{V}_{A0}}{{r}_{0}}{\left(\displaystyle \frac{{r}_{0}}{r}\right)}^{2}{L}_{* }^{\pm }\\ &&\quad +\displaystyle \frac{1}{r}\left({L}_{* }^{\pm }-\displaystyle \frac{{L}_{D}^{* }}{2}\right)\langle {{{z}^{\infty }}^{\pm }}^{2}{\rangle }^{1/2}=0;\end{eqnarray} \tag{ 105 }$$

$$\begin{eqnarray}&&U\displaystyle \frac{{{dL}}_{D}^{* }}{{dr}}-2\left(b-\displaystyle \frac{1}{2}\right)\displaystyle \frac{U}{r}{L}_{D}^{* }+2\left(3b-\displaystyle \frac{1}{2}\right)({L}_{* }^{+}+{L}_{* }^{-})\displaystyle \frac{U}{r}\\ &&\quad -(4b-1)({L}_{* }^{+}-{L}_{* }^{-})\displaystyle \frac{{V}_{A0}}{{r}_{0}}{\left(\displaystyle \frac{{r}_{0}}{r}\right)}^{2}\\ &&\quad -\displaystyle \frac{1}{r}\left[\left({{L}_{* }}^{-}-\displaystyle \frac{{L}_{D}^{* }}{2}\right)\langle {{{z}^{\infty }}^{+}}^{2}{\rangle }^{1/2}\right.\\ &&\quad \left.+\left({{L}_{* }}^{+}-\displaystyle \frac{{L}_{D}^{* }}{2}\right)\langle {{{z}^{\infty }}^{-}}^{2}{\rangle }^{1/2}\right]=0.\end{eqnarray} \tag{ 106 }$$

The system of Equations (99)–(106) is solved below for two cases, one from 0.3 to 5 au and the other for the outer heliosphere from 1 to 75 au, with b = 0.26.

The previously provided turbulence transport models are augmented by the density variance transport Equation (78). If we again assume spherical symmetry and suppose that the flow velocity corresponds to that of the solar wind with ${\boldsymbol{U}}=U\hat{{\boldsymbol{r}}}$, we can consider the two cases as follows:

$$\begin{eqnarray}&&\hat{{\boldsymbol{r}}}\parallel {{\boldsymbol{B}}}_{0}\,:\quad U\displaystyle \frac{d}{{dr}}\langle {{\rho }^{\infty }}^{2}\rangle +\displaystyle \frac{4U}{r}\langle {{\rho }^{\infty }}^{2}\rangle =-\displaystyle \frac{{\langle {{u}^{\infty }}^{2}\rangle }^{1/2}\langle {{\rho }^{\infty }}^{2}\rangle }{{{\ell }}_{u}};\end{eqnarray} \tag{ 107 }$$

$$\begin{eqnarray}&&\hat{{\boldsymbol{r}}}\perp {{\boldsymbol{B}}}_{0}\,:\quad U\displaystyle \frac{d}{{dr}}\langle {{\rho }^{\infty }}^{2}\rangle +\displaystyle \frac{4U}{r}\langle {{\rho }^{\infty }}^{2}\rangle +\displaystyle \frac{4}{r}{\langle {{u}^{\infty }}^{2}\rangle }^{1/2}\langle {{\rho }^{\infty }}^{2}\rangle \\ &&\quad =-\displaystyle \frac{{\langle {{u}^{\infty }}^{2}\rangle }^{1/2}\langle {{\rho }^{\infty }}^{2}\rangle }{{{\ell }}_{u}}.\end{eqnarray} \tag{ 108 }$$

We can evaluate $\langle {{u}^{\infty }}^{2}\rangle$ from the turbulence transport models (91)–(94) or (99)–(102) for (107) or (108), respectively, since

$$\begin{eqnarray*}&&\langle {{u}^{\infty }}^{2}\rangle =\displaystyle \frac{\langle {{z}^{\infty +}}^{2}\rangle +\langle {{z}^{\infty -}}^{2}\rangle +2{E}_{D}^{\infty }}{4}.\end{eqnarray*}$$

An important point is immediately evident, as pointed out in Zank et al. (2012b), which is that when mixing is neglected (i.e., when the RHS of both (107) and (108) and the last term on the LHS of (108) are not present), then the density fluctuation variance decays as $\langle {{\rho }^{\infty }}^{2}\rangle \propto {r}^{-4}$. This result is simply a consequence of the adiabatically expanding solar wind. However, in the absence of driving, the decay of the density fluctuation variance must proceed more rapidly than ${r}^{-4}$ within 5 au, because of the added effect of dissipation due to mixing.

One of the critical results of Zank et al. (1996) was the explicit demonstration that a driven turbulence transport model that included the dissipation of fluctuations can explain the WKB-like decay of magnetic field fluctuations. Simultaneously, the corresponding dissipative heating of the background solar wind accounts for the observed non-adiabatic radial temperature profile of the solar wind (Matthaeus et al. 1999; Smith et al. 2001). However, the precise effect of the form of the source terms on possible solutions to the driven turbulence models has not been investigated extensively. Adhikari et al. (2014) extended the Zank et al. (1996) source terms by incorporating time dependence, and Adhikari et al. (2015a, 2015b) extended the source terms to the Zank et al. (2012a) model by including source terms for the energy in the forward and backward modes and the residual energy. A drawback of the original source terms is that, for analytic convenience, they were modeled as amplification terms rather than pure source terms (unlike the pickup ion [PUI] driven turbulence source term). Although still parameterized and based on dimensionality, we consider source terms for the inner heliosphere to have the form

$$\begin{eqnarray}&&{S}_{\langle {z}^{\infty 2}\rangle }={C}_{\mathrm{sh}}^{\pm }\displaystyle \frac{{r}_{0}^{c-1}| {\rm{\Delta }}U| {V}_{A0}^{2}}{{r}^{c}};\qquad {S}_{{E}_{D}^{\infty }}={C}_{\mathrm{sh}}^{{E}_{D}}\displaystyle \frac{{r}_{0}^{c-1}| {\rm{\Delta }}U| {V}_{A0}^{2}}{{r}^{c}},\end{eqnarray} \tag{ 109 }$$

where ${C}_{\mathrm{sh}}^{+}\ne {C}_{\mathrm{sh}}^{-}\geqslant 0$; ${C}_{\mathrm{sh}}^{{E}_{D}}$ can be 0, positive, or negative; ${\rm{\Delta }}U$ is the change in flow speed between a fast and slow stream; V_A0 is a reference background Alfvén speed; and c = 2. We have associated the source terms somewhat loosely with instabilities driven by the shear of fast streams interacting with slow streams in the inner heliosphere, and we assume that the shear source of turbulence falls off with increasing heliocentric distance as ${r}^{-c}$. The shear source term is included in the core 2D transport equations only, and not the slab or NI transport equations (although in principle it could be included if we suppose that the dominant source introduces primarily Alfvënic fluctuations).

In the more distant heliosphere, shear driving is no longer as important as the driving of Alfvénic fluctuations by the creation of interstellar PUIs (Lee & Ip 1987; Williams & Zank 1994; Zank 1999). We use the PUI-driven turbulence expressions given in Adhikari et al. (2015a) Equations (19)–(21), assuming that the residual energy source term is zero (Alfvénic fluctuations). Since the fluctuations are Alfvënic, the source terms enter only the NI or slab transport equations.

In computing solutions to the radial and perpendicular models for the density variance $\langle {\rho }^{\infty 2}\rangle$ (107) and (108), we assume that both stream shear in the inner heliosphere and pickup ion creation in the outer heliosphere will generate density fluctuations according to

$$\begin{eqnarray*}&&{S}_{\mathrm{shear}}={\eta }_{1}{C}_{\mathrm{sh}}^{{\rho }^{\infty }}\displaystyle \frac{{r}_{0}^{d-1}| {\rm{\Delta }}U| }{{r}^{d}},\quad {S}_{\mathrm{PUI}}={\eta }_{2}{C}_{\mathrm{PUI}}^{{\rho }^{\infty }}\exp \left(\displaystyle \frac{-L}{r}\right),\end{eqnarray*}$$

similar to the source terms introduced by Zank et al. (2012b) with ${C}_{\mathrm{sh}}^{{\rho }^{\infty }}={\langle {{\rho }^{\infty }}^{2}\rangle }_{0}$, ${\eta }_{1}={10}^{-3}-{10}^{-2}$, and ${\eta }_{2}={10}^{-4}-{10}^{-3}$, ${C}_{\mathrm{PUI}}={\langle {{\rho }^{\infty }}^{2}\rangle }_{0}({U}_{0}/{V}_{A0}){n}_{H}^{\infty }/({\tau }_{\mathrm{ion}}^{0}{n}_{\mathrm{SW}}^{0})$ at a reference distance r₀. We take d = 3. Here, ${n}_{H}^{\infty }$ is the neutral interstellar hydrogen number density crossing the heliospheric termination shock, ${\tau }_{\mathrm{ion}}^{0}$ is the characteristic ionization time, and n_SW⁰ is a reference solar wind number density. An important parameter is the characteristic ionization cavity scale L, within which distance the neutral H population decreases significantly (see Zank 1999).

6.2. Inner Heliospheric Solutions and the Effect of Shear-Driven Source Terms

Solutions to the inner heliospheric 1D model for radial (91)–(98) and perpendicular (99)–(106) magnetic field configurations are illustrated in Figures 2–7. The boundary values that we assume at 0.29 au for this set of solutions are given in Table 1.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1.  Assumed Boundary Values at 0.29 au

2D Core Model Data		Slab Model Data
$\langle {{{z}^{\infty }}^{+}}^{2}\rangle$	13515 km2s−2	$\langle {{{z}^{* }}^{+}}^{2}\rangle$	3378.75 km2s−2
$\langle {{{z}^{\infty }}^{+}}^{2}\rangle$	753 km2s−2	$\langle {{{z}^{* }}^{-}}^{2}\rangle$	188.25 km2s−2
${E}_{D}^{\infty }$	−57.07 km2s−2	${E}_{D}^{* }$	−14.27 km2s−2
${L}_{\infty }^{+}$	1.57 × 109km3s−2	${L}_{* }^{+}$	7.84 × 108km3s−2
${L}_{\infty }^{-}$	1.6 × 108km3s−2	${L}_{* }^{+}$	8.02 × 107km3s−2
${L}_{D}^{\infty }$	−1.73 × 108km3s−2	${L}_{D}^{* }$	−8.67 × 107km3s−2
$\langle {\rho }^{\infty 2}\rangle$	113.11 cm${}^{-6}$

Download table as:  ASCIITypeset image

In general, the boundary conditions for the energy in forward and backward propagating modes (i.e., $\langle {{{z}^{\infty }}^{\pm }}^{2}\rangle$), the residual energy ${E}_{D}^{\infty }$, and the corresponding correlation functions (i.e., ${L}_{\infty }^{\pm }$ and ${L}_{D}^{\infty }$) for the 2D core model shown in Tables 1 and 2 (for the 1–75 au case) are taken from Adhikari et al. (2015a). For the slab model, the boundary conditions for the energy in forward and backward modes (i.e., $\langle {{{z}^{* }}^{\pm }}^{2}\rangle$) and the residual energy ${E}_{D}^{* }$ were obtained by assuming an 80:20 ratio for the 2D and slab energies, as discussed previously. The boundary conditions for the slab correlation functions (i.e., ${L}_{* }^{\pm }$ and ${L}_{D}^{* }$) in Tables 1 and 2 were chosen to satisfy observed values of the correlation lengths. The initial condition for the density fluctuations $\langle {\rho }^{\infty 2}\rangle$ at 1 au, Table 2, is obtained from WIND spacecraft data sets, and at 0.29 au, Table 1, $\langle {\rho }^{\infty 2}\rangle$ is obtained by assuming $\langle {\rho }^{\infty 2}\rangle \sim {r}^{-4}$ and interpolating inwards from 1 au, where r is a heliocentric distance.

Table 2.  Boundary Values at 1 au

2D Core Model Data		Slab Model Data
$\langle {{{z}^{\infty }}^{+}}^{2}\rangle$	968.33 km2s−2	$\langle {{{z}^{* }}^{+}}^{2}\rangle$	242.10 km2s−2
$\langle {{{z}^{\infty }}^{+}}^{2}\rangle$	183.88 km2s−2	$\langle {{{z}^{* }}^{-}}^{2}\rangle$	45.97 km2s−2
${E}_{D}^{\infty }$	−419.66 km2s−2	${E}_{D}^{* }$	−104.92 km2s−2
${L}_{\infty }^{+}$	7.13 × 108km3s−2	${L}_{* }^{+}$	7.13× 107km3s−2
${L}_{\infty }^{-}$	3.04 × 108km3s−2	${L}_{* }^{-}$	6.08 × 107km3s−2
${L}_{D}^{\infty }$	−8.86 × 108km3s−2	${L}_{D}^{* }$	−8.86 × 107km3s−2
$\langle {\rho }^{\infty 2}\rangle$	0.8 cm${}^{-6}$

Download table as:  ASCIITypeset image

To illustrate the nature of the solutions to the transport equations, Figures 2–4 compare the transport of purely decaying turbulence to that of shear-driven turbulence models for a radial magnetic field configuration (i.e., (91)–(98)) from 0.29–5 au. For stream shear-driven turbulence, we assume ${C}_{\mathrm{sh}}^{+}=1.0={C}_{\mathrm{sh}}^{-}$ and consider two possibilities: ${C}_{\mathrm{sh}}^{{E}_{D}}=\pm 2$, where the positive case describes a source that generates primarily velocity fluctuations and the negative primarily magnetic field fluctuations. Although we do not try to explicitly match the observational data derived by Adhikari et al. (2015a), we plot the derived observational values for some of the turbulence quantities (denoted by the "+" symbols). Our focus here is on understanding the nature of the solutions for reasonable boundary conditions at 0.29 au. Figures 5–7 correspond to Figures 2–4, except that a perpendicular magnetic field configuration is assumed—that is, Equations (99)–(106) are solved from 0.29 to 5 au. For each case, we include the appropriate density variance transport Equations (107) or (108). The PUI source terms are therefore neglected. We find very little difference between the radial and perpendicular models between 0.29 and 5 au. The solid lines in the figures show various turbulence quantities for the core 2D ("$\infty$") variables, and the dashed lines depict the corresponding quantities for the NI or slab ("∗") variables.

The energy in forward and backward propagating Elsässer fluctuations falls off with increasing heliocentric distance. There is little difference in the rate of decay of $\langle {z}^{+2}\rangle$ (referring to the 2D $\langle {z}^{\infty +2}\rangle$ and NI $\langle {z}^{* +2}\rangle$ variables and the total energy density $\langle {z}^{+2}\rangle \equiv \langle {z}^{\infty +2}\rangle +\langle {z}^{* +2}\rangle$) between the two driven models, whereas the undriven model decays more rapidly. The energy in core backward modes $\langle {z}^{\infty -2}\rangle$ experiences a slight increase within 1 au, due to the generation of these modes by forward modes as they propagate through the expanding solar wind. The energy in $\langle {z}^{\infty -2}\rangle$ for the driven solutions is about five times larger than that of the decaying case. The absence of direct source terms for the slab turbulence ensures that the decaying and driven solutions are very similar for the energy in forward and backward propagating NI Elsässer fluctuations.

The normalized residual energy ${\sigma }_{D}^{\infty }\equiv {E}_{D}^{\infty }/{E}_{T}^{\infty }$ tends to −1 within 5 au for the 2D core decaying solution and the driven ${C}_{\mathrm{sh}}^{{E}_{D}}\lt 0$ solution, indicating the increasing dominance of magnetic energy over kinetic energy for both the radial and perpendicular magnetic field configurations. A similar decrease is exhibited by the minority slab component. By contrast, the driven solution with ${C}_{\mathrm{sh}}^{{E}_{D}}\gt 0$ (i.e., a source term that introduces primarily velocity fluctuations) has ${\sigma }_{D}^{\infty }\to +1$. The slab component, because we assumed no driving, follows the decay of the undriven case. Despite observations that generally suggest otherwise (e.g., Roberts et al. 1987a, 1987b), it is sometimes thought that equipartition between kinetic and magnetic energy (the "Alfvén effect") should occur. The behavior of ${\sigma }_{D}^{\infty ,* }$ indicates that equipartition between kinetic and magnetic energy does not occur, and this point is reinforced by the plot of the Alfvén ratio ${r}_{A}^{\infty ,* }$ (Equation (74)) in Figures 2(g) and 5(g). The Alfvén ratio ${r}_{A}^{\infty }$ tends to 0 with increasing heliocentric distance for the decaying solution and for solutions with ${C}_{\mathrm{sh}}^{{E}_{D}}\lt 0$, or becomes larger than 1 if ${C}_{\mathrm{sh}}^{{E}_{D}}\gt 0$. Figures 2 and 5 illustrate that magnetic energy is dominant, even within 1 au, and certainly beyond for both the decaying and ${C}_{\mathrm{sh}}^{{E}_{D}}\lt 0$ cases. If stream shear drives predominantly velocity fluctuations, then the kinetic energy dominates within 5 au. This is illustrated further in the plots for $\langle {u}^{\infty 2}\rangle$ and $\langle {B}^{\infty 2}\rangle$ in Figures 2 and 5, which show that one or the other dominates, depending on whether ${C}_{\mathrm{sh}}^{{E}_{D}}\mathop{\gt }\limits^{\lt }0$. Little difference exists in the minority slab fluctuations.

The plots of the normalized cross helicity ${\sigma }_{c}\equiv {E}_{C}/{E}_{T}$ show that minority slab and decaying majority component fluctuations propagate primarily outward, whereas ${\sigma }_{c}^{\infty }$ for the majority driven 2D component tends to zero with increasing heliocentric distance.

The decrease in the total energy ${E}_{T}^{\infty }$ of the driven core component fluctuations is virtually the same for both stream-shear-driven cases, and of course decreases more slowly than the decay case with increasing heliocentric distance.

Overall, the theoretical curves for the ${C}_{\mathrm{sh}}^{{E}_{D}}\lt 0$ example and the corresponding observed values are in reasonable agreement between 0.29 and 5 au.

Figures 3 and 6 show ${L}_{\infty ,* }^{\pm }$, ${L}_{D}^{\infty ,* }$, and the corresponding correlation lengths. For decaying and driven turbulence, ${\lambda }_{\perp ,* }^{+}$ increases with increasing heliocentric distance. The correlation length ${\lambda }_{\perp }^{-}$ exhibits a small peak close to 1 au (when backward modes are generated), after which it increases with increasing heliocentric distance. Similarly ${\lambda }_{D}^{\infty }$, although not a well-defined quantity, decreases initially before increasing with increasing heliocentric distance for the decaying and ${C}_{\mathrm{sh}}^{{E}_{D}}\lt 0$ cases (i.e., for the magnetically dominated models). The theoretical correlation length curves for ${C}_{\mathrm{sh}}^{{E}_{D}}\lt 0$ are also in reasonable accord with their observed values.

Finally, Figures 4 and 7 show the velocity, magnetic field, and mixed velocity-magnetic field fluctuation correlation lengths ℓ_u, ℓ_b, and ℓ_ub. All three correlation lengths increase monotonically with increasing heliocentric distance for the decaying solutions. Similarly, the magnetic correlation length ℓ_b increases monotonically with increasing heliocentric radial distance for all three models. By contrast, ℓ_u for the driven model ${C}_{\mathrm{sh}}^{{E}_{D}}\gt 0$ decreases with increasing distance, unlike the other cases. Finally, ℓ_ub evolves almost identically for the two driven cases within 5 au.

The behavior of ℓ_u and $\langle {u}^{\infty 2}\rangle$ play an important role in the radial evolution of $\langle {\rho }^{\infty 2}\rangle$. Figures 4(d) and 7(d) show the evolution of the density fluctuation variance $\langle {\rho }^{\infty 2}\rangle$ with heliocentric distance r (the solid curves). The dashed line is a reference line showing an ${r}^{-4}$ slope that corresponds to purely adiabatic expansion of the density fluctuation variance. Within ∼0.4 au, all models exhibit a faster than ${r}^{-4}$ decrease in $\langle {\rho }^{\infty 2}\rangle$, since the core 2D or majority velocity fluctuations are still sufficiently strong to drive the rapid dissipative mixing of density fluctuations within this distance. After the initial rapid decay, $\langle {\rho }^{\infty 2}\rangle$ continues to dissipate beyond ∼0.5 au at a rate that is only slightly different than ${r}^{-4}$, reflecting the gradual weakening of $\langle {u}^{\infty 2}\rangle$ with increasing heliocentric distance and the increasing dominance of the magnetic field fluctuations. The undriven case tends to an ${r}^{-4}$ decrease with increasing r as ℓ_u increases. The driven ${C}_{\mathrm{sh}}^{{E}_{D}}\lt 0$ model slowly flattens because of the continued slow addition of $\langle {\rho }^{\infty 2}\rangle$, despite the radial decrease in ℓ_u. This is because the change in ℓ_u is governed by the transport equations for $\langle {z}^{\infty \pm 2}\rangle$ and ${E}_{D}^{\infty }$ and not by any self-similar behavior between $\langle {\rho }^{\infty 2}\rangle$ and ℓ_u—that is, unlike the fully developed turbulence models embracing a Kolmogorov- or IK-phenomenology, the correlation length ℓ_u does not adjust in response to the driving of density fluctuations. For the ${C}_{\mathrm{sh}}^{{E}_{D}}\gt 0$ model, $\langle {\rho }^{\infty 2}\rangle$ tends to a ${r}^{-4}$ adiabatic expansion decay law after about 0.6 au, because ℓ_u grows rapidly with distance.

Finally, we note that there is very little difference between the radial and perpendicular models, suggesting that the results within the 1D spherically symmetric formulation are rather robust.

6.3. Outer Heliospheric Solutions and the Effect of Shear-Driven and PUI Source Terms

We conclude by presenting solutions to a system of 1D transport equations from 1 to 75 au under the assumption that the interplanetary magnetic field is perpendicular to the supersonic radial flow. The same shear source terms are assumed as noted previously, and we incorporate driving of slab turbulence by the creation of pickup ions (PUIs) in the outer heliosphere. The boundary values at 1 au that we use for the solutions are given in Table 2. Consequently, unlike the inner heliospheric model of Section 6.2., the transport equations for the NI variables now acquire a source term; this is because, as discussed previously, the PUI ring-beam instability (Lee & Ip 1987; Williams & Zank 1994) generates Alfvén waves and not 2D fluctuations (see Equations (19)–(21) of Adhikari et al. (2015a) for the specific source terms; we use the same parameters).

Figure 8 follows the same format as Figure 2, and we include the observational data derived by Adhikari et al. (2015a, 2015b) to illustrate the correspondence of the data and the theoretical models. We do not present a detailed parameter analysis to best fit the theory and observations, since this requires further data analysis; however, this will be addressed by L. Adhikari et al. (2016, in preparation). The same values for ${C}_{\mathrm{sh}}^{\pm }=1$ and ${C}_{\mathrm{sh}}^{{E}_{D}}=\pm 2$ have been adopted, and the same three cases (one undriven and two driven models) are plotted in Figures 8–10.

Zoom In Zoom Out Reset image size

The heliospheric evolution of the energy densities $\langle {z}^{\infty \pm 2}\rangle$ and ${E}_{T}^{\infty }$ are very similar for both driven cases. A noticeable difference from the inner heliospheric cases is in the evolution of slab turbulence with increasing heliocentric distance. PUI-driven Alfvén waves, because of the perpendicular interplanetary magnetic configuration, are introduced equally in both directions (${\sigma }_{c}^{* }=0$ for the driven Alfvén waves), contributing equally to $\langle {z}^{* \pm 2}\rangle$. Consequently, $\langle {z}^{* \pm 2}\rangle$ begins to increase with heliocentric distance after ∼3–4 au, tending to a plateau beyond about 40–50 au. Since the slab component contributes ∼20% of the total energy, a modest increase with heliocentric distance in $\langle {z}^{\pm 2}{\rangle }_{\mathrm{tot}}$ occurs beyond about 30 au. The theoretical curves compare reasonably with the observed values.

The residual energy for the majority 2D and minority slab components behaves quite differently in the outer heliosphere than within ∼5 au. Within ∼8–10 au (roughly the size of the ionization cavity), ${\sigma }_{D}^{\infty }$ for the ${C}_{\mathrm{sh}}^{{E}_{D}}\lt 0$ case (i.e., stream shear drives primarily 2D magnetic fluctuations) is very close to −1. That is, turbulence within the ionization cavity is primarily 2D and magnetic. This is true too of the slab component within ∼2–3 au. However, as the PUI-driven slab turbulence source terms begin to strengthen, ${\sigma }_{D}^{* }$ increases to 0 by about 8–10 au. With the nonlinear transfer of slab fluctuations to 2D turbulence, ${\sigma }_{D}^{\infty }$ also increases, from −1 to a little less than 0 by 75 au. Correspondingly, the theoretical Alfvén ratios ${r}_{A}^{* ,\infty }$ exhibit similar decreases followed by an increase with increasing heliocentric distance. There is a large spread in the observed values of ${\sigma }_{D}$, but the undriven and driven ${C}_{\mathrm{sh}}^{{E}_{D}}\gt 0$ models are certainly inconsistent with the data. Observed values of r_A are more tightly constrained and are approximately constant in the outer heliosphere. The Alfvén ratios ${r}_{A}^{* ,\infty }$ are also more consistent with the driven ${C}_{\mathrm{sh}}^{{E}_{D}}\lt 0$ model than the other models. Thus the "Alfvén effect" appears to hold in the outer heliosphere but only because of the presence of Alfvén waves introduced by interstellar PUI-driving. The "Alfvén effect" is not a consequence of evolving turbulence, as demonstrated explicitly by the undriven solutions plotted in Figure 8.

The evolution of the cross helicities ${\sigma }_{c}^{\infty ,* }$ to 0 with increasing heliocentric distance is evident, both from the models and the data. Both driven models decay faster than the undriven model, but the initial rapid decrease is arrested by PUI-driving of turbulence in the outer heliosphere.

Although the ${C}_{\mathrm{sh}}^{{E}_{D}}\gt 0$ solution provides a better fit to the observed $\langle {u}^{\infty 2}\rangle$, the PUI-driven ${C}_{\mathrm{sh}}^{{E}_{D}}\lt 0$ solution increases after ∼3 au to be more consistent with observations. However, the theoretical evolution of $\langle {B}^{\infty ,* 2}\rangle$ for the ${C}_{\mathrm{sh}}^{{E}_{D}}\lt 0$ case follows the observed values very closely, as it decreases with increasing heliocentric distance. These two results suggest that a detailed fitting of the theory to observations may be quite sensitive to the exact form and heliocentric dependence of the shear source terms.

Figure 9 shows plots of ${L}_{\infty ,* }^{\pm }$ and ${L}_{D}^{\infty ,* }$ and the corresponding correlation lengths when properly defined. For both driven cases, ${\lambda }_{\perp }^{\pm }$ evolve almost identically from 1 to 75 au. Furthermore, beyond a few au, ${\lambda }_{\perp }^{+}\sim {\lambda }_{\perp }^{-}$, suggesting that this is a reasonable assumption to simplify the model equations. The same is true of ${\lambda }_{* }^{\pm }$, which increases initially within ∼3–4 au, after which turbulence driven by PUIs leads to a decrease and eventual plateau in the evolution of ${\lambda }_{* }^{\pm }$. The evolution of the related correlation lengths ℓ_u, ℓ_b, and ℓ_ub is illustrated in Figure 10, and their behavior is understood easily on the basis of PUI-driven turbulence in the outer heliosphere.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Finally, the evolution of the fluctuating density variance $\langle {\rho }^{\infty 2}\rangle$ is predicted to be quite different between the two driven cases. The ${C}_{\mathrm{sh}}^{{E}_{D}}\lt 0$ solution follows the adiabatic solution ${r}^{-4}$ very closely within the ionization cavity (as does the undriven case). In this regard, the density turbulence has "frozen," since kinetic turbulence levels are too low to actively mix and dissipate the density fluctuations. With the reintroduction of velocity fluctuations due to turbulence driven by PUIs, the density fluctuations can again experience mixing and therefore dissipate beyond ∼10 au. The decay of $\langle {\rho }^{\infty 2}\rangle$ therefore departs from the adiabatic ${r}^{-4}$ power law, decreasing more rapidly despite the addition of NI density fluctuations in the outer heliosphere as a consequence of the creation of PUIs. Relating the observed evolution of $\langle {\rho }^{\infty 2}\rangle$ to the theoretical models (Adhikari et al. 2016, in preparation) will be an important test to identify the nature of the turbulence driven by stream shear. (For instance, is it primarily kinetic or magnetic energy turbulence that is excited?)

A detailed development of the theory of nearly incompressible MHD turbulence for homogeneous and inhomogeneous flows has been presented. We re-expressed the original dimensionless form of the NI MHD equations in dimensional units for the plasma beta $\beta \sim O(1)$ or $\ll 1$ regimes, showing that NI MHD in these regimes comprises a superposition of majority 2D and minority slab fluctuations. Both the homogeneous and inhomogeneous NI MHD systems of equations were expressed in terms of either the core 2D (majority) or slab (minority) Elsässer variables, allowing us to investigate spectral and wave characteristics of the basic homogeneous system and to develop transport formalisms for the evolution of various turbulence quantities in an expanding flow. Our basic conclusions can be summarized in terms of three categories.

• A.  
  Homogeneous NI MHD conclusions:

  • 1.  
    For $\beta \sim 1$ or $\ll 1$, the majority fluctuating component is fully 2D in planes perpendicular to the large-scale mean magnetic field, being advected by the large-scale background flow. The Alfvén velocity is absent, ensuring that the majority component is composed entirely of zero-frequency properly 2D fluctuations. Since nonlinear interactions only occur between 2D modes, the spectrum for the inertial range of the energy density for Elsässer modes is Kolmogorov in the perpendicular wave number ${k}_{\perp }=| {{\boldsymbol{k}}}_{\perp }|$ (i.e., $\propto {k}_{\perp }^{-5/3}$; Equations (25) and footnote).
  • 2.  
    2D magnetic islands or vortex structures are explicitly predicted to be a nonlinear component of the dominant core 2D incompressible MHD fluctuations. Observations of magnetic islands on a variety of scales in the supersonic solar wind have been reported, providing some support for this expectation.
  • 3.  
    For the homogeneous NI component, we show that the component of vorticity parallel to the mean magnetic field satisfies the 1D Alfvén wave equation, where the spatial derivative is defined by the mean magnetic field direction. The remaining NI modes are magnetosonic. By considering only the incompressible NI corrections, we find that all three components of the vorticity satisfy 1D Alfvén wave equations, thus explicitly demonstrating that the incompressible NI corrections correspond to a minority slab component. NI MHD in the $\beta \sim 1$ or $\ll 1$ limits is therefore a superposition of 2D plus slab fluctuations.
  • 4.  
    The NI Elsässer formulation shows that the slab components are mixed passively by the coupling of the higher-order NI corrections to the leading-order 2D incompressible fluctuations, this being the dominant nonlinear interaction. The nonlinear interaction of counter-propagating Alfvén wave packets, which we have included in the NI transport equation, enters at the next order only, but is retained in order to isolate the "richest" transport equation. Analogous with the results of Shebalin et al. (1983), the nonlinear coupling of counter-propagating Alfvén waves in the NI formalism proceeds via the generation of zero-frequency 2D modes. This nonlinear interaction is a source of 2D fluctuations in the leading-order 2D incompressible description, serving to couple the NI description back to the majority component.
  • 5.  
    In considering the spectral characteristics of the slab component, we identified the inverse triple correlation time as the sum of the inverse passive scalar timescale associated with the leading-order 2D Elsässer variables and the inverse Alfvén timescale. Rather remarkably, the "critical balance" parameter falls out of the triple correlation time as the term that determines whether the slab energy flux in wave number space is dominated by either passive scalar convection by leading-order 2D Elsässer fluctuations or Alfvén wave sweeping. This is a quite different physical interpretation of the critical balance parameter than that of Goldreich & Sridhar (1995), and does not emerge as a conjecture.
  • 6.  
    Use of the triple correlation timescale and the leading-order 2D spectrum allows us to derive the anisotropic slab energy spectrum ${E}^{* }({k}_{\perp },{k}_{\parallel }){k}_{\perp }^{2}$ (Equation (53)), in which the critical balance parameter enters by demarcating the wave number regime that has a Kolmogorov spectrum from the region with an IK spectrum. The complete wave number spectrum can be expressed as a superposition of the majority 2D and minority slab spectra (Equation (54)). The composite Elsässer energy spectrum is anisotropic and does not possess a simple power law in either k_⊥ or k_∥, sometimes exhibiting a somewhat complex "broken power law" structure.
  • 7.  
    Although NI density fluctuations enter at different orders in the turbulent Mach number for the homogeneous and inhomogeneous formulations ($O({M}_{s0}^{2})$ and $O({M}_{s0})$, respectively), NI MHD shows that density fluctuations behave as passive scalars in response to advection by the majority 2D incompressible velocity fluctuations. This is in contrast to the conjecture by Lithwick & Goldreich (2001) that ascribes the mixing of density fluctuations to shear Alfvén waves. The predicted isotropic and anisotropic wave number spectra for the fluctuating density variance are $\propto {k}_{\perp }^{-5/3}$ (Equation (59)) or $\propto {k}_{\perp }^{-2/3}{k}_{\parallel }^{-1}$ (Equation (61)), respectively. If at some point the majority 2D component of the magnetic energy dominates the kinetic energy (i.e., ${\sigma }_{D}^{\infty }\simeq -1$), the density turbulence "freezes" into a non-evolving statistical state.
• B.  
  Inhomogeneous NI MHD conclusions:

  • 1.  
    A transport Equation (63) for the majority 2D forward and backward Elsässer variables was derived for flows with a large-scale inhomogeneous background in velocity, magnetic field, pressure, and density. In some respects, the new transport equation is similar to a previously derived transport equation (Marsch & Tu 1989, 1990a, 1990b; Zhou & Matthaeus 1990b), which is, however, correct only for a $\beta \gg 1$ plasma. The transport equation appropriate to $\beta \sim 1$ or $\ll 1$ (63) differs from the previous form in that the Alfvén speed and Alfvén propagation effects are absent (as is therefore the "Alfvén critical point," ${\boldsymbol{U}}=\pm {{\boldsymbol{V}}}_{A0}$) and the fluctuations are properly 2D.
  • 2.  
    The taking of suitable "moments" of the majority 2D inhomogeneous transport equation for the Elsässer variables yields a complete transport formulation for leading-order turbulence variables such as the energy density in forward and backward propagating Elsässer variables, the total energy density, the residual energy, the cross-helicity and Alfvén ratio, the kinetic and magnetic energy densities, and the correlation lengths for the Elsässer energies, residual energy, and kinetic and magnetic energy densities.
  • 3.  
    By showing that the majority 2D incompressible fluctuations drive only low-frequency compressible fluctuations (a generalization of the Lighthill mechanism for the generation of sound to MHD) and not incompressible higher-order fluctuations, we derived an inhomogeneous transport equation describing the evolution of the slab (incompressible) Elsässer variables ${{\boldsymbol{z}}}^{* \pm }$. As for the homogeneous case, the ${{\boldsymbol{z}}}^{* \pm }$ are mixed passively by the majority lower-order 2D Elsässer fluctuations ${{\boldsymbol{z}}}^{\infty \pm }$. Unlike the inhomogeneous transport equation for 2D fluctuations ${{\boldsymbol{z}}}^{\infty \pm }$, the slab transport equation contains Alfvén velocity terms ${{\boldsymbol{V}}}_{A0}$, including the "Alfvén critical point," ${\boldsymbol{U}}=\pm {{\boldsymbol{V}}}_{A0}$.
  • 4.  
    In a similar vein, we took suitable moments of the slab transport equation to derive a complete description of the evolution of slab turbulence variables such as the energy density in forward and backward propagating Elsässer variables, the total energy density, the residual energy, the cross-helicity and Alfvén ratio, the kinetic and magnetic energy densities, and the slab correlation lengths for the Elsässer energies, residual energy, and kinetic and magnetic energy densities.
  • 5.  
    The final system of turbulence transport equations for the majority 2D fluctuations (110)–(111) includes a source term that describes the dissipative generation of 2D fluctuations by the nonlinear dissipation of slab turbulence induced by counter-propagating Alfvén wave packets. Because of the inclusion of the higher-order nonlinear terms in the slab transport equation, the 2D and slab transport equations are fully coupled in both directions (i.e., the 2D fluctuations mix the slab fluctuations passively and the nonlinearity of the slab fluctuation couplings generates 2D fluctuations), unlike the original NI MHD formulation.
  • 6.  
    A transport equation describing the evolution of the variance of the $O({M}_{s})$ density fluctuations was derived. The density fluctuations are dissipated by mixing due to the leading-order 2D velocity fluctuations, and experience adiabatic expansion in a divergent flow.
• C.  
  Spherically symmetric 1D solution conclusions:The full coupled system of transport equations describing the evolving majority 2D and minority slab turbulence variables was reduced to a simpler spherically symmetric 1D set of coupled equations under the assumption that the large-scale interplanetary magnetic field is oriented either parallel or perpendicularly to the large-scale background radially expanding flow. The flow was assumed to be supersonic and of constant velocity with parameters, representative of the supersonic solar wind. A corresponding pair of transport equations for the density variance in a spherically expanding radial flow was similarly derived.

  • 1.  
    To accommodate the more complete set of transport equations, a new set of turbulence source terms due to stream shear in the supersonic solar wind was derived. These differ from those used previously, which treated the source terms as amplification terms. The new expressions allow for sources that drive unequal intensities in forward and backward directions and can generate predominantly either kinetic or magnetic energy. For the present, we assumed that stream shear drove only 2D turbulence, but this assumption is easily relaxed. Turbulence driven by pickup ions in the outer heliosphere is assumed to be in the form of Alfvén waves, implying that, unlike previous treatments, a source term for turbulence driven by PUIs in the slab turbulence transport equations is necessary. Source terms associated with both stream shear and PUI creation were included in the density variance transport equation.
  • 2.  
    The spherically symmetric 1D models were solved for the inner (0.29–5 au) and outer heliosphere (1–75 au). For the inner heliosphere, little difference was found between the solutions for the radial and perpendicular interplanetary magnetic field configurations, suggesting that the model results are rather robust. Despite the inclusion of turbulence source terms, the majority energy densities of the forward and backward Elsässer variables, the total energy, the kinetic energy, and magnetic energy all decay more slowly than the undriven solutions. However, in the outer heliosphere, turbulence generated by PUIs causes the various slab energies to increase in the outer heliosphere, with increasing heliocentric distance.
  • 3.  
    Within ∼5 au, the normalized residual energy ${\sigma }_{D}$ is sensitive to whether turbulence energy generated by stream shear is primarily magnetic (${C}_{\mathrm{sh}}^{{E}_{D}}\lt 0$) or kinetic (${C}_{\mathrm{sh}}^{{E}_{D}}\gt 0$). For the former case, the majority normalized residual energy tends to −1 within 5 au, implying that the turbulence is predominantly magnetic rather than kinetic. For the latter case, ${\sigma }_{D}^{\infty }\to +1$, and the majority 2D turbulence component is predominantly kinetic rather than magnetic. In the absence of sources of turbulence (i.e., purely decaying turbulence), ${\sigma }_{D}^{\infty }\to -1$. Thus for neither decaying nor driven turbulence can equipartition between magnetic or kinetic energy occur in the solar wind within ∼5 au (i.e., within the ionization cavity).
  • 4.  
    The normalized cross helicity of the majority 2D fluctuations tends to 0 within 5 au.
  • 5.  
    In our limited comparison to data, the solutions for the driven ${C}_{\mathrm{sh}}^{{E}_{D}}\lt 0$ case most closely conform to observations between 0.29 and 5 au.
  • 6.  
    In the absence of density fluctuation sources and dissipation by mixing, the variance in the density fluctuations decays adiabatically as ${r}^{-4}$, where r denotes heliocentric radial distance. Within 5 au, the variance in the density fluctuations decays more rapidly than ${r}^{-4}$, initially because the intensity of the kinetic energy $\langle {u}^{\infty 2}\rangle$ in 2D fluctuations ensures dissipation of the density fluctuation variance by mixing. However, if the magnetic energy comes to dominate, the rate of dissipation of $\langle {\rho }^{\infty 2}\rangle$ slows, the variance "freezes" statistically, and $\langle {\rho }^{\infty 2}\rangle$ expands adiabatically as ${r}^{-4}$.
  • 7.  
    Solutions analogous to those from 0.29 to 5 au were presented for the outer heliosphere, from 1 to 75 au. These solutions now include the driving of slab turbulence generated by the creation of PUIs. The decay of the various energies is in reasonable accord with the observations of Adhikari et al. (2015a). An important difference from the inner heliosphere solutions involves the behavior of the normalized residual energy. The inclusion of PUI-driven turbulence leads to the approximate equipartition of magnetic and kinetic energy in the distant heliosphere for both the majority 2D and minority slab components. This is despite the source term being Alfvénic and therefore entering the slab turbulence transport equations only. However, the subsequent nonlinear generation of zero-frequency 2D modes by driven slab turbulence drives the majority fluctuations toward approximate equipartition beyond ∼30 au. It is important to emphasize that the "Alfvén effect" in the outer heliosphere is a direct consequence of turbulence generation by PUIs and not the natural evolution of turbulence in an expanding supersonic flow.
  • 8.  
    Like the inner heliosphere, the normalized cross helicity tends to zero with increasing heliocentric distance.
  • 9.  
    The fluctuating density variance is found to evolve very differently according as ${C}_{\mathrm{sh}}^{{E}_{D}}\lt 0$ or $\gt 0$. The former leads to $\langle {\rho }^{\infty 2}\rangle$, evolving as ${r}^{-4}$ until ∼8 au, at which point the kinetic energy in the 2D fluctuations generated by PUI-driven turbulence is large enough to "thaw" the density turbulence and allow for further mixing and dissipation. By contrast, the kinetic energy dominated solution leads to rapid and strong dissipation of $\langle {\rho }^{\infty 2}\rangle$ within 8 au, which is eventually arrested by the slow addition of new density fluctuations generated by turbulence driven by PUI creation.

We have presented an extensive study of NI MHD turbulence for both homogeneous and inhomogeneous flows. The theory provides a comprehensive framework within which to investigate a broad range of solar wind observations in a manner not hitherto possible with the decomposition of fluctuations into a consistent majority 2D and minority slab description. The models presented here also provide a framework to investigate turbulence in sub-Alfvénic $\beta \ll 1$ and ∼1 flows such as the solar corona.

We acknowledge the partial support of NASA grants NNX08AJ33G, Subaward 37102-2, NNX14AC08G, NNX14AJ53G, A99132BT, RR185-447/4944336, and NNX12AB30G. G.P.Z. thanks R. Bruno for his kind hospitality while visiting the INAF-IAPS.

The dominant 2D incompressible MHD Equations (2)–(4) admit a class of exact nonlinear solutions that are the analogue of hydrodynamic vortices. We present a more specialized derivation of these vortical solutions, essentially following the analysis of Kadomtsev & Pogutse (1973), Verkhoglyadova et al. (2003), and Alexandrova (2008).

On introducing the vector potential ${{\boldsymbol{A}}}^{\infty }$ by ${{\boldsymbol{B}}}^{\infty }\equiv {\rm{\nabla }}\,\times {{\boldsymbol{A}}}^{\infty }={\rm{\nabla }}\times {A}^{\infty }\hat{z}$, the induction Equation (4) shows that ${{\boldsymbol{A}}}^{\infty }$ is a passive scalar in the 2D flow, satisfying

$$\begin{eqnarray}&&\displaystyle \frac{\partial {A}^{\infty }}{\partial t}+{{\boldsymbol{u}}}^{\infty }\cdot {\rm{\nabla }}{A}^{\infty }=0.\end{eqnarray} \tag{ 110 }$$

The flow vorticity is expressed as ${\xi }^{\infty }\equiv {\rm{\nabla }}\times {{\boldsymbol{u}}}^{\infty }$, and a flux function ${{\rm{\Psi }}}^{\infty }$ can be introduced through ${{\boldsymbol{u}}}^{\infty }\equiv \hat{z}\times {\rm{\nabla }}{{\rm{\Psi }}}^{\infty }$. It then follows that

$$\begin{eqnarray}&&\displaystyle \frac{\partial {A}^{\infty }}{\partial t}+\{{{\rm{\Psi }}}^{\infty },{A}^{\infty }\}=0,\end{eqnarray} \tag{ 111 }$$

where $\{a,b\}\equiv (\partial a/\partial x)(\partial b/\partial y)-(\partial a/\partial y)(\partial b/\partial x)$ denotes a Poisson bracket. Since the current ${{\boldsymbol{J}}}^{\infty }$ is defined by ${\mu }_{0}{{\boldsymbol{J}}}^{\infty }={\rm{\nabla }}\times {{\boldsymbol{B}}}^{\infty }$, for 2D flows we find that ${\xi }^{\infty }\parallel {{\boldsymbol{J}}}^{\infty }$ and

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\xi }^{\infty }}{\partial t}+{{\boldsymbol{u}}}^{\infty }\cdot {{\rm{\nabla }}}_{\perp }{\xi }^{\infty }=\displaystyle \frac{1}{\rho }{{\boldsymbol{B}}}^{\infty }\cdot {\rm{\nabla }}{{\boldsymbol{J}}}^{\infty }.\end{eqnarray} \tag{ 112 }$$

Since ${\xi }_{z}^{\infty }={{\rm{\nabla }}}_{\perp }^{2}{{\rm{\Psi }}}^{\infty }$ and ${\mu }_{0}{J}_{z}^{\infty }=-{{\rm{\nabla }}}_{\perp }^{2}{A}^{\infty }$, we can rewrite (112) as

$$\begin{eqnarray}\displaystyle \frac{\partial {\xi }_{z}^{\infty }}{\partial t}+\{{{\rm{\Psi }}}^{\infty },{\xi }_{z}^{\infty }\} & = & -\,\displaystyle \frac{1}{\rho }\{{A}^{\infty },{J}_{z}^{\infty }\};\\ \displaystyle \frac{\partial }{\partial t}({{\rm{\nabla }}}_{\perp }^{2}{{\rm{\Psi }}}^{\infty })+\{{{\rm{\Psi }}}^{\infty },{{\rm{\nabla }}}_{\perp }^{2}{{\rm{\Psi }}}^{\infty }\} & = & \displaystyle \frac{1}{{\mu }_{0}\rho }\{{A}^{\infty },{{\rm{\nabla }}}_{\perp }^{2}\}.\end{eqnarray} \tag{ 113 }$$

For non-propagating structures, Equations (111) and (113) become

$$\begin{eqnarray}\{{{\rm{\Psi }}}^{\infty },{A}^{\infty }\} & = & 0;\\ \{{{\rm{\Psi }}}^{\infty },{{\rm{\nabla }}}_{\perp }^{2}{{\rm{\Psi }}}^{\infty }\} & = & \displaystyle \frac{1}{{\mu }_{0}\rho }\{{A}^{\infty },{{\rm{\nabla }}}_{\perp }^{2}\},\end{eqnarray} \tag{ 114 }$$

the first condition of which ensures that ${A}^{\infty }=f({{\rm{\Psi }}}^{\infty })$ for an arbitrary function f. The second condition of (114) implies that

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{\perp }^{2}{{\rm{\Psi }}}^{\infty }=-\displaystyle \frac{1}{\rho }{f}^{\prime }{J}_{z}^{\infty }+{f}_{1}({{\rm{\Psi }}}^{\infty }),\end{eqnarray} \tag{ 115 }$$

where f₁ is another arbitrary function. We consider the simplest spherically symmetric solution of (115) by choosing ${A}^{\infty }={{\rm{\Psi }}}^{\infty }$ (which implies ${f}^{\prime }=1$) and suppose ${J}_{z}^{\infty }$ is linear in ${A}^{\infty }$—say ${\mu }_{0}{J}_{z}^{\infty }={k}^{2}{A}^{\infty }$ for some constant k. On rewriting (115) in polar coordinates, subject to the previous assumptions, (115) reduces to Bessel's equation of order 0,

$$\begin{eqnarray}&&\displaystyle \frac{{d}^{2}{A}^{\infty }}{{{dx}}^{2}}+\displaystyle \frac{1}{x}\displaystyle \frac{{{dA}}^{\infty }}{{dx}}+{A}^{\infty }=0,\end{eqnarray} \tag{ 116 }$$

after setting $x=\alpha r$, where $\alpha =k/\sqrt{{\mu }_{0}\rho }$. Solutions are therefore ${A}^{\infty }(r)={{CJ}}_{0}(k/\sqrt{{\mu }_{0}\rho }r)$, where J₀ is Bessel's function of order 0. The parameter k is chosen to ensure that $A({r}_{0})=0$. This then yields the vortical solution (26).