A GENERALIZED TWO-COMPONENT MODEL OF SOLAR WIND TURBULENCE AND AB INITIO DIFFUSION MEAN-FREE PATHS AND DRIFT LENGTHSCALES OF COSMIC RAYS

We begin by introducing our notation for the large-scale and small-scale fields. The total solar wind velocity is written ${\boldsymbol{U}}({\boldsymbol{r}})+{\boldsymbol{v}}({\boldsymbol{r}},{\boldsymbol{x}})$, the sum of a large-scale piece dependent upon the heliocentric position vector ${\boldsymbol{r}}$, and a small-scale contribution that also depends upon local small-scale coordinates ${\boldsymbol{x}}$, relative to each ${\boldsymbol{r}}$. Similarly, the total magnetic field is ${\boldsymbol{B}}({\boldsymbol{r}})+{\boldsymbol{b}}({\boldsymbol{r}},{\boldsymbol{x}})$, with associated large-scale Alfvén speed ${{\boldsymbol{V}}}_{{\rm{A}}}={\boldsymbol{B}}/\sqrt{4\pi \rho }$, where $\rho ({\boldsymbol{r}})$ is the large-scale mass density. The small-scale dynamics is treated as incompressible (see Zank et al. 2012b for a discussion of transport of density fluctuations). As a simplifying assumption, the fluctuation amplitudes, ${\boldsymbol{v}}$ and ${\boldsymbol{b}}$, are restricted to be transverse to ${\boldsymbol{B}};$ that is, parallel variances are neglected. Solar wind observations indicate that this is often a reasonable approximation (e.g., Belcher & Davis 1971; Klein et al. 1991). In general, the above quantities are also time-dependent.

The large-scale wind velocity ${\boldsymbol{U}}$ is with respect to an inertial frame; in the frame corotating with the Sun the large-scale velocity is ${\boldsymbol{V}}={\boldsymbol{U}}-{\boldsymbol{\Omega }}\times {\boldsymbol{r}}$, where ${\boldsymbol{\Omega }}$ is the solar angular rotation rate. Our numerical computations are often performed in this corotating frame. In obtaining the transport equations in this frame, we make use of the relation ${\rm{\nabla }}\cdot {\boldsymbol{V}}={\rm{\nabla }}\cdot {\boldsymbol{U}}$, which holds because ${\rm{\nabla }}\cdot ({\boldsymbol{\Omega }}\times {\boldsymbol{r}})=0$.

The two-component aspect of the model involves separating the fluctuations into two precisely defined incompressible elements: quasi-2D turbulence and a complementary wave-like component (Oughton et al. 2006, 2011). Specifically, employing Elasser variables, ${{\boldsymbol{z}}}^{\pm }={\boldsymbol{v}}\pm {\boldsymbol{b}}/\sqrt{4\pi \rho ({\boldsymbol{r}})}$, we express the fluctuations as

$$\begin{eqnarray}&&{{\boldsymbol{z}}}^{\pm }({\boldsymbol{r}},{\boldsymbol{x}})={{\boldsymbol{q}}}^{\pm }+{{\boldsymbol{w}}}^{\pm },\end{eqnarray} \tag{ 1 }$$

where ${{\boldsymbol{q}}}^{\pm }$ and ${{\boldsymbol{w}}}^{\pm }$ are the quasi-2D and wave-like components, respectively; both quantities are functions of the (large-scale) heliocentric radius ${\boldsymbol{r}}$ and the small-scale displacements ${\boldsymbol{x}}$ from each ${\boldsymbol{r}}$.

Table 1 summarizes the definitions of the major energy-related fluctuation quantities, which appear in the transport model. For the quasi-2D component, ${\sigma }_{c,z}$ is the normalized cross helicity, and ${\sigma }_{D}^{z}$ the normalized energy difference, equal to the (normalized) kinetic energy less the magnetic energy all divided by the sum of these. In general, the analogous quantity for the wave-like component is indicated by a subscript or superscript w.

Table 1.  Definitions of Some Important Physical Variables for the Quasi-2D and Wave-like Components

Quasi-2D		Wave-like
Fluctuations	Quantity	Fluctuations
${Z}_{\pm }^{2}=\langle {{\boldsymbol{q}}}_{\pm }\cdot {{\boldsymbol{q}}}_{\pm }\rangle$	Elasser "energies"	${W}_{\pm }^{2}=\langle {{\boldsymbol{w}}}_{\pm }\cdot {{\boldsymbol{w}}}_{\pm }\rangle$
$2{Z}^{2}={Z}_{+}^{2}+{Z}_{-}^{2}$	total "energies"	$2{W}^{2}={W}_{+}^{2}+{W}_{-}^{2}$
$2{H}_{c}^{z}={Z}_{+}^{2}-{Z}_{-}^{2}$	cross helicities	$2{H}_{c}^{w}={W}_{+}^{2}-{W}_{-}^{2}$
${\sigma }_{c,z}=\displaystyle \frac{{Z}_{+}^{2}-{Z}_{-}^{2}}{{Z}_{+}^{2}+{Z}_{-}^{2}}$	normalized cross helicities	${\sigma }_{c,w}=\displaystyle \frac{{W}_{+}^{2}-{W}_{-}^{2}}{{W}_{+}^{2}+{W}_{-}^{2}}$
${\sigma }_{D}^{z}=\displaystyle \frac{\langle {{\boldsymbol{q}}}_{+}\cdot {{\boldsymbol{q}}}_{-}\rangle }{{Z}^{2}}$	normalized energy differences	${\sigma }_{D}^{w}=\displaystyle \frac{\langle {{\boldsymbol{w}}}_{+}\cdot {{\boldsymbol{w}}}_{-}\rangle }{{W}^{2}}$

Note. Angle brackets $\langle \cdots \rangle$ indicate averaging over the small-scale coordinate ${\boldsymbol{x}}$ (at each large-scale coordinate ${\boldsymbol{r}}$). Note that ${H}_{c}^{z}$ and ${H}_{c}^{w}$ differ by a factor of two from the definitions used in Oughton et al. (2011).

Download table as:  ASCIITypeset image

Along with the energies (per mass) of the fluctuations, ${Z}_{\pm }^{2}$ and ${W}_{\pm }^{2}$, it is also necessary to consider their characteristic lengthscales, typically defined using correlation lengths. In general, these are distinct for each type of field; for example, ${{\ell }}_{+}$ for ${Z}_{+}^{2}$ and ${{\ell }}_{-}$ fo ${Z}_{-}^{2}$. Here we make the simplifying assumption that these scales are equal and denote the characteristic lengthscale of Z² as ${\ell }$ and that of W² as λ. In addition, the typical parallel scale of the wave-like component, ${\lambda }_{\parallel }$, is needed, particularly in connection with driving by pickup ions. (For one-component transport models that consider the ± lengthscales separately see Zank et al. 2012a and Adhikari et al. 2015.)

Finally, in this section, we address the suitability of using incompressible MHD to model solar wind fluctuations. Naturally, the actual solar wind fluctuations will often display some compressive activity. Here, however, from the outset we approximate them as being incompressible and thus neglect small-scale compressive behavior. On the observational side, density fluctuations are often found to be ∼10% of the mean value (e.g., Roberts et al. 1987; Matthaeus et al. 1991), providing motivation for neglecting compressive activity at this level. On the theory side, the nearly incompressible (NI) approach for systems with small Mach numbers (Zank & Matthaeus 1992, 1993), leads to a leading-order description that is either incompressible 3D MHD (large plasma beta) or incompressible 2D MHD (beta small or order unity). The next order corrections are termed "NI" and support MHD waves. In particular, when beta is of the order of unity, as is typical for the solar wind, the NI solutions include Alfvén waves with timescales shorter than those associated with the leading-order incompressible behavior. Thus, modeling the system as we do herein, i.e., using incompressible quasi-2D and incompressible wave-like components, is consistent with the NI results.

The transport and driving terms—for the energy, cross helicity, and characteristic lengthscales of the fluctuations—have been derived and discussed in various works (e.g., Matthaeus et al. 1994; Usmanov et al. 2011; Zank et al. 2012a). Here we largely follow the approach of Matthaeus et al. (1994) and Usmanov et al. (2011), extended to incorporate the homogeneous two-component phenomenology presented in Oughton et al. (2011) and also retaining terms of the order of ${V}_{{\rm{A}}}/U$ (Adhikari et al. 2015; Wiengarten et al. 2015).

This leads to the following equations for the fluctuation energies, in the frame corotating with the Sun,

$$\begin{eqnarray}\displaystyle \frac{\partial {Z}^{2}}{\partial t} & = & -\,{\rm{\nabla }}\cdot ({\boldsymbol{V}}{Z}^{2}+{{\boldsymbol{V}}}_{{\rm{A}}}{H}_{c}^{z})+2{{\boldsymbol{V}}}_{{\rm{A}}}\cdot {\rm{\nabla }}{H}_{c}^{z}+\displaystyle \frac{1}{2}({\rm{\nabla }}\cdot {\boldsymbol{U}}){Z}^{2}\\ & & -{\sigma }_{D}^{z}{Z}^{2}\left[\displaystyle \frac{{\rm{\nabla }}\cdot {\boldsymbol{U}}}{2}-\hat{{\boldsymbol{B}}}\cdot (\hat{{\boldsymbol{B}}}\cdot {\rm{\nabla }}){\boldsymbol{U}}\right]\\ & & -{\alpha }_{z}\left[\displaystyle \frac{{Z}^{3}}{{\ell }}{f}_{{zz}}^{+}+\displaystyle \frac{2{{WZ}}^{2}}{{\ell }}\displaystyle \frac{{f}_{{zw}}^{+}}{1+Z/W}\right]\ +\ {\alpha }_{z}{X}^{+}\\ & & \ +\ \displaystyle \frac{{Z}^{2}}{r}{C}_{\mathrm{sh}}^{Z}| {\boldsymbol{U}}| ,\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}\displaystyle \frac{\partial {W}^{2}}{\partial t} & = & -\,{\rm{\nabla }}\cdot ({\boldsymbol{V}}{W}^{2}+{{\boldsymbol{V}}}_{{\rm{A}}}{H}_{c}^{w})+2{{\boldsymbol{V}}}_{{\rm{A}}}\cdot {\rm{\nabla }}{H}_{c}^{w}+\displaystyle \frac{1}{2}({\rm{\nabla }}\cdot {\boldsymbol{U}}){W}^{2}\\ & & -{\sigma }_{D}^{w}{W}^{2}\left[\displaystyle \frac{{\rm{\nabla }}\cdot {\boldsymbol{U}}}{2}-\hat{{\boldsymbol{B}}}\cdot (\hat{{\boldsymbol{B}}}\cdot {\rm{\nabla }}){\boldsymbol{U}}\right]\\ & & \ -\ {\alpha }_{w}\left[\displaystyle \frac{2{W}^{2}Z}{\lambda }\displaystyle \frac{{f}_{{wz}}^{+}}{1+\lambda /{\ell }}+\displaystyle \frac{2{W}^{4}{\lambda }_{\parallel }}{{\lambda }^{2}{V}_{{\rm{A}}}}(1-{\sigma }_{c,w}^{2})\right]-{\alpha }_{z}{X}^{+}\\ & & \ +\ \displaystyle \frac{{W}^{2}}{r}{C}_{\mathrm{sh}}^{W}| {\boldsymbol{U}}| \ +\ {\dot{E}}_{\mathrm{PI}},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}\displaystyle \frac{\partial {H}_{c}^{z}}{\partial t} & = & -\,{\rm{\nabla }}\cdot ({\boldsymbol{V}}{H}_{c}^{z}+{{\boldsymbol{V}}}_{{\rm{A}}}{Z}^{2})+2{{\boldsymbol{V}}}_{{\rm{A}}}\cdot {\rm{\nabla }}{Z}^{2}+\displaystyle \frac{1}{2}({\rm{\nabla }}\cdot {\boldsymbol{U}}){H}_{c}^{z}\\ & & +{\sigma }_{D}^{z}{Z}^{2}\left[{\rm{\nabla }}\cdot {{\boldsymbol{V}}}_{{\rm{A}}}+\displaystyle \frac{\hat{{\boldsymbol{B}}}\cdot (\hat{{\boldsymbol{B}}}\cdot {\rm{\nabla }}){\boldsymbol{B}}}{\sqrt{4\pi \rho }}\right]\\ & & \ -\ {\alpha }_{z}\left[\displaystyle \frac{{Z}^{3}}{{\ell }}{f}_{{zz}}^{-}+\displaystyle \frac{2{{WZ}}^{2}}{{\ell }}\displaystyle \frac{{f}_{{zw}}^{-}}{1+Z/W}\right]+{\alpha }_{z}{X}^{-},\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}\displaystyle \frac{\partial {H}_{c}^{w}}{\partial t} & = & -\,{\rm{\nabla }}\cdot ({\boldsymbol{V}}{H}_{c}^{w}+{{\boldsymbol{V}}}_{{\rm{A}}}{W}^{2})+2{{\boldsymbol{V}}}_{{\rm{A}}}\cdot {\rm{\nabla }}{W}^{2}+\displaystyle \frac{1}{2}({\rm{\nabla }}\cdot {\boldsymbol{U}}){H}_{c}^{w}\\ & & +{\sigma }_{D}^{w}{W}^{2}\left[{\rm{\nabla }}\cdot {{\boldsymbol{V}}}_{{\rm{A}}}+\displaystyle \frac{\hat{{\boldsymbol{B}}}\cdot (\hat{{\boldsymbol{B}}}\cdot {\rm{\nabla }}){\boldsymbol{B}}}{\sqrt{4\pi \rho }}\right]\\ & & \ -\ {\alpha }_{w}\left[\displaystyle \frac{2{W}^{2}Z}{\lambda }\displaystyle \frac{{f}_{{wz}}^{-}}{1+\lambda /{\ell }}\right]-{\alpha }_{z}{X}^{-},\end{eqnarray} \tag{ 5 }$$

where

$$\begin{eqnarray}&&{X}^{\pm }={Y}^{+}\pm {Y}^{-},\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}{Y}^{\pm } & = & {W}_{\pm }{Z}_{\pm }\left[\displaystyle \frac{{Z}_{\mp }}{\lambda }{{\rm{\Gamma }}}_{w}^{{z}_{\mp }{w}_{\pm }}+\displaystyle \frac{{W}_{\mp }}{\lambda }{{\rm{\Gamma }}}_{w}^{{w}_{\mp }{w}_{\pm }}\right.\\ & & -\left.\displaystyle \frac{{Z}_{\mp }}{{\ell }}{{\rm{\Gamma }}}_{w}^{{z}_{\mp }{z}_{\pm }}-\displaystyle \frac{{W}_{\mp }}{{\ell }}{{\rm{\Gamma }}}_{w}^{{w}_{\mp }{z}_{\pm }}\right],\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{c}^{{ab}}=\displaystyle \frac{1}{1+{\tau }_{\mathrm{nl}}^{{ab}}/{\tau }_{{\rm{A}}}^{c}},\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&2{f}_{{ab}}^{\pm }\,=\,(1+{\sigma }^{a})\sqrt{1-{\sigma }^{b}}\pm (1-{\sigma }^{a})\sqrt{1+{\sigma }^{b}},\end{eqnarray} \tag{ 9 }$$

with ${\sigma }^{a}\equiv {\sigma }_{c,z}$ or ${\sigma }_{c,w}$ for the component a = Z or W. Equation (9) defines various "f" functions, bounded by ±1. These act as attenuation factors for the modeled nonlinear terms when the cross helicities are non-zero, as is appropriate (see, e.g., Dobrowolny et al. 1980a).

The ${Y}^{+}$ term, which may be positive or negative, models the exchange of excitation between ${Z}_{+}^{2}$ and ${W}_{+}^{2}$, and similarly for ${Y}^{-}$. The ${{\rm{\Gamma }}}_{c}^{{ab}}$ are associated with the decay rate of the triple correlation for the term being modeled, and involve the nonlinear (${\tau }_{\mathrm{nl}}$) and Alfvén (${\tau }_{{\rm{A}}}$) timescales of the appropriate components. Further details are given in Oughton et al. (2006).

Structurally, we have written Equations (2)–(5) so that different types of physics appear on separate lines. On the first lines, we have advection, expansion, and propagation effects (essentially the WKB terms). The "mixing" terms, proportional to a ${\sigma }_{D}$ (Zhou & Matthaeus 1990), are on the second lines. The third line in each equation presents the homogeneous decay phenomenology terms. If there is any forcing, the terms modeling those effects appear as a fourth line. For example, the quasi-2D and wave-like energies are driven by large-scale velocity shear—modeled using either self-consistently computed velocity gradients (Wiengarten et al. 2015) or ad hoc terms in the manner of earlier models (e.g., Zank et al. 1996; Breech et al. 2008)—and W² is also forced by waves generated during the near isotropization of pickup ions (${\dot{E}}_{\mathrm{PI}}$).

Note that the right-most mixing terms in Equations (4) and (5) are absent from the model of Zank et al. (2012a) on setting their suggested structural similarity parameters for axisymmetric quasi-2D fluctuations, namely $a=1/2$, b = 0. We find, however, that in order to recover the model of Matthaeus et al. (1994) it is appropriate to choose $a=b=1/2$.

Since in each of the Equations (2)–(5) the final line arises from a turbulence phenomenology (Oughton et al. 2011), the terms on these lines are only determined to within O(1) multiplying constants. There are some constraints on these constants; for example, when adding the Z² and W² equations, we require that the exchange terms cancel. Here, we adopt the simplest approach of using a single constant in each equation (except for the variations required in connection with the exchange terms), denoted ${\alpha }_{z}$ and ${\alpha }_{w}$.

Transport equations for the characteristic lengthscales—${\ell }$, λ, ${\lambda }_{\parallel }$—are derived following the approach of Matthaeus et al. (1994). This is based on integrating correlation functions over the (small-scale) lag, ${\boldsymbol{\xi }}$. For example, in the case of ${\ell }$, one starts with transport equations for ${R}^{\pm }({\boldsymbol{r}},{\boldsymbol{\xi }})=\langle {{\boldsymbol{q}}}^{\pm }({\boldsymbol{r}},{\boldsymbol{x}})\cdot {{\boldsymbol{q}}}^{\pm }({\boldsymbol{r}},{\boldsymbol{x}}+{\boldsymbol{\xi }})\rangle$, defines ${L}_{\pm }={\int }_{0}^{\infty }{R}^{\pm }({\boldsymbol{r}},{\boldsymbol{\xi }})\,{\rm{d}}\xi$ and obtains their transport equations, adds these to give an equation for $L={L}_{+}+{L}_{-}=2{Z}^{2}{\ell }$, and then extracts the equation for ${\ell }$. The choice of integration direction, $\hat{{\boldsymbol{\xi }}}$, is discussed below. (See Zank et al. 2012a for a distinct approach.) With the extension to two components and retention of $O({{\boldsymbol{V}}}_{{\rm{A}}})$ terms, this leads to

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\ell }}{\partial t}=-\,{\boldsymbol{V}}\cdot {\rm{\nabla }}{\ell }+{\sigma }_{c,z}{{\boldsymbol{V}}}_{{\rm{A}}}\cdot {\rm{\nabla }}{\ell }\\ &&\,+\displaystyle \frac{{L}_{D}}{{Z}^{2}}\left[{\rm{\nabla }}\cdot \displaystyle \frac{{\boldsymbol{U}}}{2}-2{\hat{\xi }}_{i}{\hat{\xi }}_{j}\displaystyle \frac{\partial {U}_{i}}{\partial {r}_{j}}\right]\\ &&\,+{\ell }{\sigma }_{D}^{z}\left[{\rm{\nabla }}\cdot \displaystyle \frac{{\boldsymbol{U}}}{2}-\hat{{\boldsymbol{B}}}\cdot (\hat{{\boldsymbol{B}}}\cdot {\rm{\nabla }}){\boldsymbol{U}}\right]\\ &&\,+\ {\beta }_{z}\left[{f}_{{zz}}^{+}Z+{f}_{{zw}}^{+}\displaystyle \frac{2W}{1+Z/W}-\displaystyle \frac{{\ell }{X}^{+}}{{Z}^{2}}\right],\,\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}\displaystyle \frac{\partial \lambda }{\partial t} & = & -\,{\boldsymbol{V}}\cdot {\rm{\nabla }}\lambda +{\sigma }_{c,w}{{\boldsymbol{V}}}_{{\rm{A}}}\cdot {\rm{\nabla }}\lambda \\ & & +\displaystyle \frac{{\tilde{L}}_{D}}{{W}^{2}}\left[{\rm{\nabla }}\cdot \displaystyle \frac{{\boldsymbol{U}}}{2}-2{\hat{\xi }}_{i}{\hat{\xi }}_{j}\displaystyle \frac{\partial {U}_{i}}{\partial {r}_{j}}\right]\\ & & +\lambda {\sigma }_{D}^{w}\left[{\rm{\nabla }}\cdot \displaystyle \frac{{\boldsymbol{U}}}{2}-\hat{{\boldsymbol{B}}}\cdot (\hat{{\boldsymbol{B}}}\cdot {\rm{\nabla }}){\boldsymbol{U}}\right]\\ & & \ +\ {\beta }_{w}\left[\displaystyle \frac{2{f}_{{wz}}^{+}Z}{1+\lambda /{\ell }}+2(1-{\sigma }_{c,w}^{2})\displaystyle \frac{{W}^{2}{\lambda }_{\parallel }}{\lambda {V}_{{\rm{A}}}}+\displaystyle \frac{{\alpha }_{z}\lambda {X}^{+}}{{\alpha }_{w}{W}^{2}}\right],\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\lambda }_{\parallel }}{\partial t}=-\,{\boldsymbol{V}}\cdot {\rm{\nabla }}{\lambda }_{\parallel }+{\sigma }_{c,w}{{\boldsymbol{V}}}_{{\rm{A}}}\cdot {\rm{\nabla }}{\lambda }_{\parallel }\\ &&\,+0\quad (\mathrm{mixing}\ \mathrm{terms}\ \mathrm{cancel})\\ &&\,+2{\alpha }_{w}(1-{\sigma }_{c,w}^{2})\displaystyle \frac{{W}^{2}{\lambda }_{\parallel }}{{V}_{{\rm{A}}}{\lambda }^{2}}{\lambda }_{\parallel }\\ &&\,-\ ({\lambda }_{\parallel }-{\lambda }_{\mathrm{res}})\displaystyle \frac{{\dot{E}}_{\mathrm{PI}}}{{W}^{2}}.\,\end{eqnarray} \tag{ 12 }$$

Again the presentation structure has advection, expansion, and wave propagation terms on the first lines, mixing terms (when nonzero) on the second and third lines, turbulence phenomenology terms on the line after those, and any forcing on a separate line after that. In the general case, terms associated with shear driving also appear in the lengthscale equations (e.g., Matthaeus et al. 1996; Zank et al. 1996, 2012a; Breech et al. 2008; Oughton et al. 2011). Herein, however, we assume that shear driving occurs at the correlation scales and thus ${\ell }$, λ, and ${\lambda }_{\parallel }$ are unaffected by such forcings.

Because ${\ell }$ and λ are characteristic transverse lengthscales, in Equations (10) and (11), the unit vector $\hat{{\boldsymbol{\xi }}}$ must be chosen to lie in the plane perpendicular to ${\boldsymbol{B}}$, i.e., in the plane of the fluctuation amplitudes. For a ${\boldsymbol{B}}$ that lies in the R-T plane, such as the Parker spiral field, a useful choice is $\hat{{\boldsymbol{\xi }}}=\hat{{\boldsymbol{\vartheta }}}$, where ϑ is the polar angle in heliocentric spherical coordinates. (See Matthaeus et al. 1994, where $\hat{{\boldsymbol{\xi }}}$ is denoted as $\hat{{\boldsymbol{r}}}$.)

In general, one also needs equations for the energy difference lengthscales (Matthaeus et al. 1994; Zank et al. 2012a; Adhikari et al. 2015). Here we employ the closures ${L}_{D}={\ell }{\sigma }_{D}^{z}{Z}^{2}$ and ${\tilde{L}}_{D}=\lambda {\sigma }_{D}^{w}{W}^{2}$. These imply equality of the correlation lengths for the velocity and magnetic fields (${{\ell }}_{v}={{\ell }}_{b};$ ${\lambda }_{v}={\lambda }_{b}$), and induce slight simplifications of Equations (10) and (11).

In obtaining the equation for ${\lambda }_{\parallel }$, we assume that the correlation functions for the W component have the same symmetry structure as that for "slab" Alfvén waves and integrate along the mean field direction: $\hat{{\boldsymbol{\xi }}}=\hat{{\boldsymbol{B}}}$. We also make the approximation of a single parallel lengthscale, e.g., ${\lambda }_{\parallel ,D}={\lambda }_{\parallel }$. These features combine to cause cancellation of the mixing terms. The energy injection associated with (near) isotropization of pickup-ion-induced waves occurs at the gyroradius of the pickup protons, ${\lambda }_{\mathrm{res}}({\boldsymbol{r}})=2\pi U({\boldsymbol{r}})/{{\rm{\Omega }}}_{{\rm{p}}}({\boldsymbol{r}})$ with the proton gyrofrequency Ω.

To close the model, assuming that the large-scale fields like ${\boldsymbol{V}}$ and ${{\boldsymbol{V}}}_{{\rm{A}}}$ are known, we require knowledge of the normalized energy differences, ${\sigma }_{D}^{z}$, ${\sigma }_{D}^{w}$. Their transport equations are obtained in a similar fashion to the above derivations (Matthaeus et al. 1994; Zank et al. 2012a; Adhikari et al. 2015). Herein, however, we approximate ${\sigma }_{D}^{z}$ and ${\sigma }_{D}^{w}$ as constant parameters, on the basis of rough observational support (Roberts et al. 1987; Perri & Balogh 2010; Iovieno et al. 2016). This yields a closed set of equations for the fluctuations, given the large-scale fields. Transport equations for the latter are now considered.

The fluctuations in the present model consist of two different components. This leads to some modified terms in the large-scale momentum equation. The single fluctuation component form is given in Usmanov et al. (2011), see their Equation (B2), as

$$\begin{eqnarray}&&\displaystyle \frac{\partial (\rho {\boldsymbol{U}})}{\partial t}+{\rm{\nabla }}\cdot \left[\rho {\boldsymbol{V}}{\boldsymbol{U}}-\displaystyle \frac{\eta }{4\pi }{\boldsymbol{B}}{\boldsymbol{B}}+\left(P+\displaystyle \frac{{B}^{2}}{8\pi }+{p}_{\mathrm{fluct}}\right){\boldsymbol{I}}\right]\\ &&\quad =-\rho ({\boldsymbol{g}}+{\boldsymbol{\Omega }}\times {\boldsymbol{U}}),\end{eqnarray} \tag{ 13 }$$

where ${\boldsymbol{\Omega }}={\rm{\Omega }}{{\boldsymbol{e}}}_{z};$ ${\rm{\Omega }}=14\buildrel{\circ}\over{.} 71\,{\mathrm{day}}^{-1}$ (Snodgrass & Ulrich 1990), and ${\boldsymbol{g}}=({{GM}}_{\odot }/{r}^{2})\,{{\boldsymbol{e}}}_{r}$ describes the Sun's gravitational acceleration, and P is the large-scale gas pressure.

The forms of η and the pressure of the fluctuations ${p}_{\mathrm{fluct}}$ depend upon the assumed symmetries of the latter, e.g., via the modeling of the MHD Reynolds stress (Usmanov et al. 2011). For the present (transverse, axisymmetric) two-component case, they become

$$\begin{eqnarray}&&{\eta }^{2\mathrm{cpt}}=1+\displaystyle \frac{{\sigma }_{D}^{z}{Z}^{2}+{\sigma }_{D}^{w}{W}^{2}}{2{V}_{{\rm{A}}}^{2}},\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&{p}_{\mathrm{fluct}}^{2\mathrm{cpt}}=(1+{\sigma }_{D}^{z})\displaystyle \frac{\rho {Z}^{2}}{4}+(1+{\sigma }_{D}^{w})\displaystyle \frac{\rho {W}^{2}}{4},\end{eqnarray} \tag{ 15 }$$

and are used in place of η and ${p}_{\mathrm{fluct}}$ in Equation (13), which is otherwise unchanged. "Cross-component" effects like $\langle {{\boldsymbol{b}}}_{Z}{{\boldsymbol{b}}}_{W}\rangle$ with ${\boldsymbol{b}}={{\boldsymbol{b}}}_{Z}+{{\boldsymbol{b}}}_{W}$ have been neglected. Note that (15) is equivalent to the kinetic (not magnetic) pressure of the fluctuations, though this is a little misleading (physically) since the term is actually the sum of the fluctuation magnetic pressure and contributions from modeling of the MHD Reynolds stresses (Usmanov et al. 2011).

An equation for the total energy density is straightforward to obtain (e.g., Usmanov et al. 2011). However, due to a feature of the Cronos code, we work instead with the energy density associated with unforced ideal MHD,

$$\begin{eqnarray}&&e=\displaystyle \frac{\rho {U}^{2}}{2}+\displaystyle \frac{{B}^{2}}{8\pi }+\displaystyle \frac{P}{\gamma -1},\end{eqnarray} \tag{ 16 }$$

where the full energy density also includes gravitational potential energy and the turbulence energy, $\rho ({Z}^{2}+{W}^{2})/2$. These "missing" terms in e are accounted for using source terms in the energy equation (Wiengarten et al. 2015, Appendix B). Following the latter approach with a $\gamma =5/3$ adiabatic equation of state and Hollweg's heat flux ${{\boldsymbol{q}}}_{H}$ (Hollweg 1974, 1976) yields

$$\begin{eqnarray}&&{\partial }_{t}e+{\rm{\nabla }}\cdot \left[e{\boldsymbol{V}}+\left(P+\displaystyle \frac{| {\boldsymbol{B}}{| }^{2}}{8\pi }\right){\boldsymbol{U}}-({\boldsymbol{U}}\cdot {\boldsymbol{B}})\displaystyle \frac{{\boldsymbol{B}}}{4\pi }\right.\\ &&\quad -\left.{{\boldsymbol{V}}}_{{\rm{A}}}\rho \displaystyle \frac{{H}_{c}}{2}+{{\boldsymbol{q}}}_{H}\right]\\ &&\quad =-\rho {\boldsymbol{V}}\cdot {\boldsymbol{g}}-{\boldsymbol{U}}\cdot {\rm{\nabla }}{p}_{\mathrm{fluct}}^{2\mathrm{cpt}}\\ &&\,-\displaystyle \frac{{H}_{c}}{2}{{\boldsymbol{V}}}_{{\rm{A}}}\cdot {\rm{\nabla }}\rho -\rho {{\boldsymbol{V}}}_{{\rm{A}}}\cdot {\rm{\nabla }}{H}_{c}\\ &&\,+{\boldsymbol{U}}\cdot ({\boldsymbol{B}}\cdot {\rm{\nabla }})\left[({\eta }^{2\mathrm{cpt}}-1)\displaystyle \frac{{\boldsymbol{B}}}{4\pi }\right]\\ &&\,+\rho \left[{\alpha }_{z}\left(\displaystyle \frac{{f}_{{zz}}^{+}{Z}^{3}}{2{\ell }}+\displaystyle \frac{{f}_{{zw}}^{+}}{1+Z/W}\displaystyle \frac{{{WZ}}^{2}}{{\ell }}\right)\right.\\ &&\,+\left.{\alpha }_{w}\left(\displaystyle \frac{{f}_{{wz}}^{+}}{1+\lambda /{\ell }}\displaystyle \frac{{{ZW}}^{2}}{\lambda }+(1-{\sigma }_{c,w}^{2})\displaystyle \frac{{W}^{4}{\lambda }_{\parallel }}{{\lambda }^{2}{V}_{{\rm{A}}}}\right)\right],\end{eqnarray} \tag{ 17 }$$

where ${H}_{c}={H}_{c}^{z}+{H}_{c}^{w}$.

The equations describing the evolution of the large-scale density and magnetic field are unaffected by the extension to incompressible two-component fluctuations. Neglecting the turbulent electric field, one has (e.g., Usmanov et al. 2011; Wiengarten et al. 2015),

$$\begin{eqnarray}&&{\partial }_{t}\rho +{\rm{\nabla }}\cdot (\rho {\boldsymbol{V}})=0,\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&{\partial }_{t}{\boldsymbol{B}}+{\rm{\nabla }}\cdot ({\boldsymbol{V}}{\boldsymbol{B}}-{\boldsymbol{B}}{\boldsymbol{V}})={\bf{0}}.\end{eqnarray} \tag{ 19 }$$

In order to validate the implementation in Cronos, we compare results obtained in Oughton et al. (2011) with those from an appropriately restricted form of the new model's equations. Specifically, the background solar wind is prescribed to be a uniform and constant radial flow with ${\boldsymbol{U}}=440\,\mathrm{km}\,{{\rm{s}}}^{-1}\,{{\boldsymbol{e}}}_{r}$, and a proton number density profile $n={n}_{0}{({r}_{0}/r)}^{2}$ where ${n}_{0}=n({r}_{0}=0.3\,\mathrm{au})=66\,{\mathrm{cm}}^{-3}$. The large-scale magnetic field is a Parker spiral, expressed in terms of a vector potential (e.g., Wiengarten et al. 2015),

$$\begin{eqnarray}&&{\boldsymbol{A}}=-{B}_{0}{r}_{0}^{2}\sin (\vartheta )\left(\displaystyle \frac{\varphi }{r}+\displaystyle \frac{{\rm{\Omega }}}{U}\right){{\boldsymbol{e}}}_{\vartheta },\end{eqnarray} \tag{ 20 }$$

where ϑ and φ are the polar and azimuthal angles in (heliocentric) spherical polar coordinates and ${B}_{0}=43$ nT. Additionally, the turbulence transport equations are relieved of all advection and mixing terms involving the Alfvén velocity (but retain the dissipation and interchange terms), as well as the advection and mixing terms in the lengthscale equations. The energy density equation, (17), simplifies considerably and can be usefully re-expressed via $P=2{nkT}$ in Equation (16) in terms of the proton temperature (Oughton et al. 2011, Equation (14)); in the present study, however, it is the energy density equation that is solved. The equations are then formally equivalent to those of Oughton et al. (2011), where the sources of turbulence considered are stream shear (modeling the influence of, e.g., corotating interaction regions) and isotropization of pickup ion distributions. While the stream shear drives both the quasi-2D and the wave-like component (so that ${C}_{\mathrm{sh}}^{Z,W}=1$, see below), the pickup-ion driving feeds the wave-like component only and is approximated as (Zank et al. 1996; Williams et al. 1997)

$$\begin{eqnarray}&&{\dot{E}}_{\mathrm{PI}}=\displaystyle \frac{\zeta {U}^{2}{n}_{H}}{{n}_{\mathrm{sw}}{\tau }_{\mathrm{ion}}}\exp \left(-\displaystyle \frac{{L}_{\mathrm{cav}}}{r}\displaystyle \frac{{\rm{\Psi }}}{\sin ({\rm{\Psi }})}\right),\end{eqnarray} \tag{ 21 }$$

where ${n}_{H}=0.1\,{\mathrm{cm}}^{-3}$ is the interstellar neutral hydrogen density, ${\tau }_{\mathrm{ion}}=1.33\times {10}^{6}\,{\rm{s}}$ is the hydrogen ionization time at 1 au, ${L}_{\mathrm{cav}}=5.6\,\mathrm{au}$ is the characteristic scale of the ionization cavity of the Sun, and ${n}_{\mathrm{sw}}=6\,{\mathrm{cm}}^{-3}$ is the solar wind density at 1 au. The angle Ψ is that between the observation point and the upwind direction; for pickup ions entering the heliosphere along the x-axis, it corresponds to heliospheric latitude, so that above the poles the effective ionization cavity is larger by a factor of $\pi /2$, and this pushes the region where pickup-ion heating is important to larger r. The factor ζ describes the fraction of the available energy actually channeled into the fluctuations and is mainly a function of the ratio of Alfvén speed to solar wind speed according to the model of Isenberg et al. (2003) and Isenberg (2005) that is used in Section 3.2. For this validation case, we assume a constant $\zeta =0.04$. The Kármán–Taylor constants are set as ${\alpha }_{z}={\alpha }_{w}=2{\beta }_{z}=2{\beta }_{w}=0.25$ and the residual energies are assumed constant with ${\sigma }_{D}^{z}={\sigma }_{D}^{w}=-1/3$ (e.g., Roberts et al. 1987; Perri & Balogh 2010).

The computational domain extends from 0.3 to 100 au and is covered with 300 cells of increasing cell size ${\rm{\Delta }}r$ from 10 to 250 solar radii, while azimuthal symmetry is assumed and the computations are restricted to the ecliptic plane. The remaining inner boundary values at ${r}_{0}=0.3\,\mathrm{au}$ are ${Z}^{2}=1500$ km² s⁻², ${W}^{2}=150$ km² s⁻², ${\sigma }_{c,z}={\sigma }_{c,w}=0.6$, $l=\lambda =0.008\,\mathrm{au}$, ${\lambda }_{\parallel }=0.036\,\mathrm{au}$ and $T=1.6\times {10}^{5}$ K.

Figure 1 shows the resulting behavior of the turbulence quantities with radial distance. In the inner heliosphere, due to shear driving both the quasi-2D and the wave-like component's energy densities decrease less steeply and normalized cross helicities drop strongly. The latter point follows because, for example, ${\sigma }_{c,z}=({Z}_{+}^{2}-{Z}_{-}^{2})/({Z}_{+}^{2}+{Z}_{-}^{2})$, and thus adding energy equally to the ${Z}_{\pm }^{2}$ leaves the numerator unchanged but increases the denominator (Matthaeus et al. 2004; Breech et al. 2005). As shear driving diminishes with heliospheric distance, the quasi-2D component decays freely while pickup-ion driving feeds only the wave-like component, which, consequently, constitutes the dominant component in the outer heliosphere with its normalized cross helicity ${\sigma }_{c,w}$ quickly going to zero and its correlation length λ much shorter than its quasi-2D counterpart l. In the case shown, the pickup driving is strong enough to induce noticable transfer of energy from W² to Z² beyond ∼40 au. This occurs via the "exchange" term, ${X}^{+}$, in Equations (2) and (3), as discussed in Oughton et al. (2011). There is also an associated decrease of ${\ell }$ at these distances. Note that the "anti-correlated" behavior of Z² and ${\ell }$ with heliocentric distance does not hold for W and λ, which is a consequence of the pickup-ion driving. Furthermore, convergence of the parallel lengthscale toward the pickup-ion gyroradius is also evident. The decaying turbulent energy is dissipated and heats the outer heliosphere as can be seen in the temperature panel. Results obtained with Cronos (black lines) are shown alongside those obtained with the IDL code (red lines) used in Oughton et al. (2011). An implementation mistake that was present in the latter has since been corrected. The agreement validates the implementation in Cronos.

Zoom In Zoom Out Reset image size

The model presented in Section 2, and employed in the remainder of the paper, extends that by Oughton et al. (2011) of the previous section in two ways. First, the background solar wind is no longer prescribed, but computed self-consistently and in a fully three-dimensional manner alongside the turbulence transport equations. Second, the latter are augmented in several ways, namely by (1) not neglecting transport and mixing terms involving the Alfvén velocity, (2) improving the stream shear driving so that it is computed from the background wind, and (3) employing the theory from Isenberg (2005) for the efficiency of pickup-ion driving.

In consequence, the implemented model is applicable to arbitrary solar wind conditions, including sub-Alfvénic heliospheric regions such as the corona and the heliosheath. Coronal models and global heliospheric simulations are both challenging in regard to computer resources, due to the high space and time resolutions required for the former and the long propagation times needed for the latter, especially when including multi-fluid aspects and magnetic fields (e.g., Scherer et al. 2016). We leave such applications for future studies and consider here the super-Alfvénic solar wind during typical solar minimum conditions of fast polar winds and a band of slow wind occupying equatorial regions. We impose azimuthal symmetry, which allows for a considerable reduction of computational costs and thereby enables coverage of the full polar angle with one degree resolution. The radial grid is the same as in the previous section, covering the distance from 0.3 to 100 au. Figure 2 displays the applied inner boundary conditions depending on colatitude. The top row shows the background quantities (velocity, number density, magnetic field strength, and temperature), in setting which we were guided by Ulysses measurements (McComas et al. 2000). This includes a small latitudinal gradient ($\approx 1\,\mathrm{km}/{\rm{s}}/^\circ$) of solar wind speed in the fast wind regime, constant mass flux, and a Parker spiral magnetic field structure that neglects a polarity reversal and current sheet. The latter would be under-resolved in these non-AMR simulations and would affect the equatorial results more strongly as appropriate. The bottom row shows the turbulence quantities (turbulent energy density, lengthscales and cross helicities). There is considerable spread and uncertainty associated with spacecraft measurements of these quantities (see Figure 5) and boundary values were chosen to give a reasonable fit to the available data, with the 90%–10% partitioning for Z²-W² guided by observation-based studies (e.g., Bieber et al. 1996; Hamilton et al. 2008). Such studies report a range of values but typically find a dominant quasi-2D component; see Oughton et al. (2015) for a recent review.

Zoom In Zoom Out Reset image size

Turbulence driven by stream shear can be calculated self-consistently from the background wind in the present setup, as introduced in Wiengarten et al. (2015). However, the influence of corotating interaction regions, present near solar minimum, is not inherently covered in this simplified geometry with azimuthal symmetry. Moreover, we find that if additional shear is not included in the high-speed regions this results in cross helicities that increase with radial distance (see Dobrowolny et al. 1980a, 1980b), which is in contrast to Ulysses measurements (Figure 5). The source of this additional shear can be attributed to so-called microstreams (Neugebauer et al. 1995). In order to model these additional effects, we include ad hoc terms ${C}_{\mathrm{add}}^{Z,W}$ in the full driving for Z and W, so that

$$\begin{eqnarray}&&{C}_{\mathrm{sh}}^{Z,W}=\displaystyle \frac{1}{| {\boldsymbol{U}}| }({\partial }_{\vartheta }+{(\sin \vartheta )}^{-1}{\partial }_{\varphi })| {\boldsymbol{U}}| +{C}_{\mathrm{add}}^{Z,W},\end{eqnarray} \tag{ 22 }$$

with ${C}_{\mathrm{add}}^{Z,W}$ chosen such that in the band of slow wind ${C}_{\mathrm{sh}}^{Z,W}=1$, while ${C}_{\mathrm{sh}}^{Z,W}=0.25$ for the fast wind, i.e., a lower bound on shear driving is imposed at all latitudes. The transition region results in higher values and the latitudinal profile of the shear driving displayed in Figure 3 is similar to that used in Breech et al. (2008).

Zoom In Zoom Out Reset image size

The other source for driving turbulence is the excitation of waves via the near isotropization of pickup-ion distributions (${\dot{E}}_{\mathrm{PI}}$), which we use here in the same form as in Equation (21), but with the efficiency factor $\zeta ({V}_{{\rm{A}}}/U,Z/{V}_{{\rm{A}}})$ calculated using the improved formulation developed in Isenberg et al. (2003), see also Isenberg (2005).

As before, the residual energy densities are assumed constant with ${\sigma }_{D}^{z}={\sigma }_{D}^{w}=-1/3$, and the Kármán–Taylor constants are taken to be ${\alpha }_{z}={\alpha }_{w}=2{\beta }_{z}=2{\beta }_{w}=0.2$ (Breech et al. 2008; Oughton et al. 2011). The low-latitude inner boundary values at ${r}_{0}=0.3\,\mathrm{au}$ are ${Z}^{2}=900$ km² s⁻², ${W}^{2}=90$ km² s⁻², ${\sigma }_{c,z}={\sigma }_{c,w}=0.4$, $l=\lambda =0.012\,\mathrm{au}$, ${\lambda }_{\parallel }=0.03\,\mathrm{au}$, and $T=3.0\times {10}^{5}$ K, while at high latitudes these values are ${Z}^{2}=5000$ km² s⁻², ${W}^{2}=500$ km² s⁻², ${\sigma }_{c,z}={\sigma }_{c,w}=0.6$, $l=\lambda =0.018\,\mathrm{au}$, ${\lambda }_{\parallel }=0.03\,\mathrm{au}$, and $T=1.5\times {10}^{6}$ K. Simulations are performed until a steady state is reached, for which the required physical time corresponds approximately to the propagation time from the inner to the outer radial boundary, i.e., about one year. The resulting configuration of the background wind is illustrated in the top row of Figure 4, along with the turbulence quantities in the middle and bottom rows, by contour plots of two-dimensional meridional slices.

Zoom In Zoom Out Reset image size

The magnetic field exhibits the typical Parker spiral behavior of decreasing more slowly in the ecliptic (∝r⁻¹) than above the poles (∝r⁻²), resulting in a constant Alfvén speed in the former and a radially decreasing one in the latter region. The solar wind speed is approximately constant along radial spokes. The background solar wind quantities are barely affected by the inclusion of a turbulence description (Wiengarten et al. 2015), except for some additional heating, mainly occurring in the fast wind/slow wind transition region due to the strong shear there, and in the outer heliosphere due to increased pickup-ion production. The latter effect essentially only acts in the ecliptic plane, because the efficiency factor ζ tends to zero for small ${\text{}}{V}_{{\rm{A}}}/{\text{}}U$, as is the case away from the ecliptic plane. This is seen best in the panel for the wave-like turbulence component, W². Also visible are the stripes of enhanced turbulence levels in the transition region, and these are even clearer in the Z² panel.

The regions with stronger generation of turbulence are associated with cross helicities quickly going to zero in their respective component. In other regions, cross helicities that are not equal to zero are retained also at large radial distances, which is not only due to the absence of sources for turbulence, but also because of the inclusion of the additional Alfvén velocity related transport terms, as already demonstrated in Wiengarten et al. (2015) for a one-component turbulence model. Furthermore, the perpendicular lengthscales increase with radial distance as turbulence decays, while the parallel lengthscale approaches the resonant one (${\lambda }_{\mathrm{res}}$), which is inversely proportional to the magnetic field strength.

Figure 5 shows comparisons of the model results at selected colatitudes with spacecraft measurements. For the fast wind regions, we use Ulysses measurements during its first fast latitude scan (Bavassano et al. 2000a, 2000b, blue crosses) picking out latitudes higher than 35°. Although there is a mixed latitudinal and radial dependence in these data, we use it for comparison with radial dependence of the model data only and choose a colatitude of 15° (blue lines). Model output in the equatorial plane (black lines) is compared with measurements from the Voyager 2 spacecraft that have been used in previous studies (Roberts et al. 1987; Zank et al. 1996; Smith et al. 2001).

Zoom In Zoom Out Reset image size

Consider first the high-latitude results. The Ulysses measurements for the turbulent energies (assumed to reside mainly in the quasi-2D component) and temperature show little scattering and are well reproduced by the model, whereas spread in the data is large for the correlation lengths and cross helicity. However, the model results are well within the covered range. In the outer heliosphere, pickup-ion driving is evident in W² and ${\sigma }_{c,w}$ at $r\gtrsim 20\,\mathrm{au}$, but only becomes significant in terms of the total fluctuation energy for $r\gtrsim 80\,\mathrm{au}$. Since shear driving is also weak in the outer heliosphere, ${\sigma }_{c,z}$ remains significantly non-zero and there is no strong heating at these high latitudes. This is in contrast to the situation near the ecliptic.

At low latitudes, shear driving is relatively strong inside $\approx 5$ au, so the radial profiles of the turbulent energies are flatter than their high-latitude counterparts. Pickup-ion driving also becomes important closer in (around 5 au) and causes the wave-like component to become the dominant one for $r\gtrsim 10\,\mathrm{au}$. This leads to a stronger cascade of fluctuation energy and the associated dissipation yields the increasing temperature profile in the outer heliosphere. Thus, it appears that an important reason for the stronger heating near the ecliptic, compared to high latitudes, is the greater radial range where pickup-ion forcing is effective. Voyager measurements show considerable spread but there is again some agreement with the (ecliptic) model results. In particular, the model temperature is a rough lower bound to the observational data and the energy-weighted lengthscale, $L=({\ell }{Z}^{2}+\lambda {W}^{2})/({Z}^{2}+{W}^{2})$, passes close to most of the ecliptic data values. Recall that here (and in Wiengarten et al. 2015), Alfvén velocity terms are retained in the transport equations. As Wiengarten et al. (2015) note, this is associated with a shallower radial decrease of ${\sigma }_{c,z}$ and ${\sigma }_{c,w}$, compared to transport models, which neglect terms of order ${V}_{{\rm{A}}}/U$. Moreover, this leads to better agreement with observational data, particularly for the energy-weighted cross helicity ${{\rm{\Sigma }}}_{c}=({Z}^{2}{\sigma }_{c,z}+{W}^{2}{\sigma }_{c,w})/({Z}^{2}+{W}^{2})$, depicted using a red dotted line in Figure 5.