Turbulent Transport in a Three-dimensional Solar Wind

The solar wind MHD model of Shiota et al. (2014) calculates the spatial and temporal variation of the 3D solar wind on the basis of a time series of observed photospheric magnetic field synoptic maps. The inner heliosphere defined in the Heliographic Inertial (HGI) coordinate system is discretized with a Yinyang grid (Kageyama & Sato 2004). In the MHD model (Shiota et al. 2014), we integrate the following MHD equations:

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho }{\partial t}+\mathrm{div}(\rho {\boldsymbol{v}})=0,\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho {\boldsymbol{v}}}{\partial t}+\mathrm{div}(\rho {\boldsymbol{vv}}-{\boldsymbol{BB}}+{p}_{T}{\boldsymbol{I}})=\rho {\boldsymbol{g}},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{B}}}{\partial t}+\mathrm{div}({\boldsymbol{Bv}}-{\boldsymbol{vB}}+\psi {\boldsymbol{I}})=0,\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial e}{\partial t}+\mathrm{div}((e+{p}_{T}){\boldsymbol{v}}-({\boldsymbol{v}}\cdot {\boldsymbol{B}}){\boldsymbol{B}})=\rho {\boldsymbol{v}}\cdot {\boldsymbol{g}}+H,\end{eqnarray} \tag{ 4 }$$

where $e=\tfrac{\rho {{\boldsymbol{v}}}^{2}}{2}+\tfrac{p}{\gamma -1}+\tfrac{{{\boldsymbol{B}}}^{2}}{2},$ ${\boldsymbol{g}}=-\tfrac{{GM}}{{r}^{3}}{\boldsymbol{r}},$ and ${p}_{T}=p+\tfrac{{{\boldsymbol{B}}}^{2}}{2}.$ These equations are normalized with typical values: ${L}_{0}\,={R}_{\odot }=6.96\times {10}^{5}$ km, ${\tau }_{0}=1\,{\rm{h}}$ = 3600 s, ${\rho }_{0}=1$ cm${}^{-3}\times {m}_{{\rm{H}}}$, ${V}_{0}={L}_{0}/{\tau }_{0}=193.3$ km/s, ${B}_{0}={V}_{0}\,\times \sqrt{4\pi {\rho }_{0}}=8.864$ nT, and ${p}_{0}=6.252\times {10}^{-2}$ nPa. I is the unit matrix in Cartesian coordinates. The additional variable ψ is integrated with the following equation

$$\begin{eqnarray}&&\displaystyle \frac{\partial \psi }{\partial t}+{c}_{h}^{2}\mathrm{div}{\boldsymbol{B}}+{c}_{p}\psi =0,\end{eqnarray} \tag{ 5 }$$

that relieves the numerical deviation of the divergence free (solenoidal) condition of magnetic field (Dedner et al. 2002), where c_h and c_p are the propagation speed of ψ defined as the maximum characteristic speed within the numerical domain and the diffusion coefficient of ψ is defined as ${c}_{p}={c}_{h}/0.18$, respectively (Dedner et al. 2002). ${m}_{{\rm{H}}}$, G, and M are the proton mass, the gravitational constant, and the solar mass, respectively. H is a heating term (Adhikari et al. 2014, 2015a) that is discussed later and is set to be 0 in this study. These equations are solved using a finite volume method with a Harten–Lax–van Leer Discontinuities approximate Riemann solver (Miyoshi & Kusano 2005), which is combined with the third-order Monotone Upstream-centered Schemes for Conservation Laws and second-order Runge–Kutta time integration.

The numerical domain in this study is set as $60{R}_{\odot }\leqslant r\,\leqslant 1290{R}_{\odot }$ (6 au). The inner boundary conditions at $r=60{R}_{\odot }$ are rotating and time-varying solar wind maps obtained from the 3D potential coronal magnetic field (Shiota et al. 2012) and empirical velocities between coronal magnetic field and solar wind parameters (Wang & Sheeley 1990; Totten et al. 1995; Arge & Pizzo 2000; Hayashi et al. 2003).

The solar wind 3D structure is determined from the time series of solar wind and IMF maps specified on the inner boundary of the MHD simulation (Shiota et al. 2014). In order to obtain a simple solar wind structure, we use two simplified steady magnetic field configurations:

• 1.  
  a monopole field structure in which the field is outward-directed and its strength (0.55 G or 0.275 G at r = 2.5 ${R}_{\odot }$) is spherically symmetric; and
• 2.  
  a tilted dipole magnetic field whose magnetic axis is inclined at 30° from the rotational axis of the Sun.

Furthermore, we use two types of solar wind velocity maps:

• 1.  
  a spherically symmetric solar wind flow with a constant speed 400 km s⁻¹ or 600 km s⁻¹; and
• 2.  
  a bimodal (fast and slow) solar wind flow whose speed is determined from the tilted dipole potential field and the Wang–Sheeley–Arge 2000 formula (Wang & Sheeley 1990; Arge & Pizzo 2000).

We ran the simulations for a sufficiently long time, so that a steady-state solution was obtained. We show numerical results with combinations of the four configurations and compare with observations of turbulence in the solar wind.

We combine the theoretical turbulence transport model developed by Zank et al. (2012) with our global MHD model (Shiota et al. 2014). In general, turbulence in the MHD regime can be treated with the MHD equations, Equations (1)–(4), if the spatial resolution of the numerical simulation is sufficiently high to resolve the higher wavenumber turbulence. However, it is difficult, if not impossible currently, to solve a global MHD simulation with such a high spatial resolution because the energy-containing rage of the MHD turbulence is too large. We separate the high frequency fluctuations in the velocity ${\boldsymbol{v}}$ and magnetic field ${\boldsymbol{B}}$ (Zhou & Matthaeus 1990), which cannot be described by the global MHD simulation, and denote the scale-separated fluctuations as ${\boldsymbol{u}}$ and ${\boldsymbol{b}}$. The fluctuations in solar wind velocity and magnetic field associated with oppositely propagating Alfvén modes are described by the Elsässar variables (Elsässer 1950):

$$\begin{eqnarray}&&{{\boldsymbol{z}}}^{\pm }={\boldsymbol{u}}\pm {\boldsymbol{b}}/\sqrt{4\pi \rho }.\end{eqnarray} \tag{ 6 }$$

The transport and modulation of the Elsässar variables are described by scale-separated incompressible MHD equations (Marsch & Tu 1989, 1990; Zhou & Matthaeus 1990):

$$\begin{eqnarray}&&\displaystyle \frac{\partial {{\boldsymbol{z}}}^{\pm }}{\partial t}+({\boldsymbol{U}}\mp {{\boldsymbol{V}}}_{A})\cdot {\rm{\nabla }}{{\boldsymbol{z}}}^{\pm }+\displaystyle \frac{1}{2}{\rm{\nabla }}\cdot ({\boldsymbol{U}}/2\pm {{\boldsymbol{V}}}_{A}){{\boldsymbol{z}}}^{\pm }\\ &&+{{\boldsymbol{z}}}^{\mp }\cdot \left[{\rm{\nabla }}{\boldsymbol{U}}\pm \displaystyle \frac{{\rm{\nabla }}{\boldsymbol{B}}}{\sqrt{4\pi \rho }}-\displaystyle \frac{{\boldsymbol{I}}}{2}{\rm{\nabla }}\cdot ({\boldsymbol{U}}/2\pm {{\boldsymbol{V}}}_{A})\right]\\ &&\,={{\boldsymbol{NL}}}_{\pm }+{{\boldsymbol{S}}}^{\pm }\end{eqnarray} \tag{ 7 }$$

where U is the mean velocity, V_A is the mean Alfvén velocity, ${{\boldsymbol{NL}}}_{\pm }$ are nonlinear dissipation terms, and ${{\boldsymbol{S}}}^{\pm }$ are source terms for the turbulence. The detail of the source terms ${{\boldsymbol{S}}}^{\pm }$ are described in the next subsection.

Equation (7) is derived from the standard incompressible reduction of the MHD Equations (1)–(4). Strictly speaking Equation (7) is therefore valid only for solar wind conditions satisfying ${\beta }_{p}\gg 1$, where ${\beta }_{p}$ is the usual plasma beta (Zank & Matthaeus 1991; Hunana & Zank 2010). Despite this limitation, the turbulence models (Zank et al. 1996, 2012; Matthaeus et al. 1999b; Smith et al. 2001; Adhikari et al. 2014, 2015a) derived from it are quite well validated by in situ observations of low-frequency turbulence in the solar wind (Adhikari et al. 2014, 2015a).

The turbulence transport model of Zank et al. (2012) considered the temporal and spatial evolution of three moments of the small-scale fluctuations (unit: energy per unit mass),

$$\begin{eqnarray}&&f\equiv \langle {{z}^{+}}^{2}\rangle =\langle {{\boldsymbol{z}}}^{+}\cdot {{\boldsymbol{z}}}^{+}\rangle ,\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&g\equiv \langle {{z}^{-}}^{2}\rangle =\langle {{\boldsymbol{z}}}^{-}\cdot {{\boldsymbol{z}}}^{-}\rangle ,\end{eqnarray} \tag{ 9 }$$

and

$$\begin{eqnarray}&&{E}_{D}\equiv \langle {{\boldsymbol{z}}}^{+}\cdot {{\boldsymbol{z}}}^{-}\rangle =\langle {u}^{2}\rangle -\langle {b}^{2}/4\pi \rho \rangle .\end{eqnarray} \tag{ 10 }$$

The transport equations of these variables are derived by taking moments of Equation (7), corresponding to Equations (42)–(44) of Zank et al. (2012):

$$\begin{eqnarray}\displaystyle \frac{\partial f}{\partial t} & = & -\,({\boldsymbol{U}}-{{\boldsymbol{V}}}_{A})\cdot {\rm{\nabla }}f-{\rm{\nabla }}\cdot \left(\displaystyle \frac{{\boldsymbol{U}}}{2}+{{\boldsymbol{V}}}_{A}\right)\,f\\ & & +{M}^{+}{E}_{D}-2\displaystyle \frac{{f}^{2}\sqrt{g}}{{L}^{+}}+2{S}^{+};\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}\displaystyle \frac{\partial g}{\partial t} & = & -({\boldsymbol{U}}+{{\boldsymbol{V}}}_{A})\cdot {\rm{\nabla }}g-{\rm{\nabla }}\cdot \left(\displaystyle \frac{{\boldsymbol{U}}}{2}-{{\boldsymbol{V}}}_{A}\right)g\\ & & +{M}^{-}{E}_{D}-2\displaystyle \frac{{g}^{2}\sqrt{f}}{{L}^{-}}+2{S}^{-};\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}\displaystyle \frac{\partial {E}_{D}}{\partial t} & = & -{\boldsymbol{U}}\cdot {\rm{\nabla }}{E}_{D}-\displaystyle \frac{{\rm{\nabla }}\cdot {\boldsymbol{U}}}{2}{E}_{D}+\displaystyle \frac{{M}^{-}}{2}f\\ & & +\displaystyle \frac{{M}^{+}}{2}g-\displaystyle \frac{1}{2}\left(\sqrt{\displaystyle \frac{f}{g}}{{\boldsymbol{V}}}_{A}\cdot {\rm{\nabla }}g-\sqrt{\displaystyle \frac{g}{f}}{{\boldsymbol{V}}}_{A}\cdot {\rm{\nabla }}f\right)\\ & & -\left[\displaystyle \frac{g\sqrt{f}}{{L}^{-}}-\displaystyle \frac{f\sqrt{g}}{{L}^{+}}\right]{E}_{D}+2{S}_{D},\end{eqnarray} \tag{ 13 }$$

where

$$\begin{eqnarray}{M}^{\pm } & = & {\rm{\nabla }}\cdot \left(\displaystyle \frac{{\boldsymbol{U}}}{2}\pm {{\boldsymbol{V}}}_{A}\right)-2\left(a\left[{\rm{\nabla }}\cdot {\boldsymbol{U}}-\displaystyle \frac{({\boldsymbol{B}}\cdot {\rm{\nabla }}{\boldsymbol{U}})}{{B}^{2}}\cdot {\boldsymbol{B}}\right]\right.\\ & & \left.\pm b\left[{\rm{\nabla }}\cdot {{\boldsymbol{V}}}_{A}-\displaystyle \frac{{\boldsymbol{B}}\cdot {\rm{\nabla }}| {V}_{A}| }{| B| }\right]\right)\end{eqnarray} \tag{ 14 }$$

and ${L}^{\pm }$ and L^D are variables related to the correlation lengths (${\lambda }_{\pm }$ and ${\lambda }_{D}$) of these three moments, respectively. They are defined as ${L}^{+}\equiv f{\lambda }_{+}$, ${L}^{-}\equiv g{\lambda }_{-}$, and ${L}^{D}\equiv {E}_{D}{\lambda }_{D}$. In the turbulence transport model (Equations (11)–(13)), we solve simultaneously for the temporal and spatial evolution of the three variables ${L}^{\pm }$ and L^D,

$$\begin{eqnarray}&&\displaystyle \frac{\partial {L}^{\pm }}{\partial t}=-({\boldsymbol{U}}\mp {{\boldsymbol{V}}}_{A})\cdot {\rm{\nabla }}{L}^{\pm }-{\rm{\nabla }}\cdot \left(\displaystyle \frac{{\boldsymbol{U}}}{2}\pm {{\boldsymbol{V}}}_{A}\right){L}^{\pm }+\displaystyle \frac{{M}^{\pm }}{2}{L}_{D};\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}\displaystyle \frac{\partial {L}^{D}}{\partial t} & = & -\,{\boldsymbol{U}}\cdot {\rm{\nabla }}{L}^{D}-\displaystyle \frac{{\rm{\nabla }}\cdot {\boldsymbol{U}}}{2}{L}^{D}+{M}^{-}{L}^{+}+{M}^{+}{L}^{-}\\ & & -\left(\sqrt{\displaystyle \frac{{L}^{+}}{{L}^{-}}}{{\boldsymbol{V}}}_{A}\cdot {\rm{\nabla }}{L}^{-}-\sqrt{\displaystyle \frac{{L}^{-}}{{L}^{+}}}{{\boldsymbol{V}}}_{A}\cdot {\rm{\nabla }}{L}^{+}\right),\end{eqnarray} \tag{ 16 }$$

also derived from Equation (7). According to the derivation in Appendix A, the parameters a and b in Equation (14) are set to $a=b=1/2$ in the present study, which yields

$$\begin{eqnarray}{M}^{\pm } & = & {\rm{\nabla }}\cdot \left(\displaystyle \frac{{\boldsymbol{U}}}{2}\pm {{\boldsymbol{V}}}_{A}\right)-\left[{\rm{\nabla }}\cdot {\boldsymbol{U}}-\displaystyle \frac{({\boldsymbol{B}}\cdot {\rm{\nabla }}{\boldsymbol{U}})}{{B}^{2}}\cdot {\boldsymbol{B}}\right]\\ & & \mp \left[{\rm{\nabla }}\cdot {{\boldsymbol{V}}}_{A}-\displaystyle \frac{{\boldsymbol{B}}\cdot {\rm{\nabla }}| {V}_{A}| }{| B| }\right]\\ & = & -\,{\rm{\nabla }}\cdot \left(\displaystyle \frac{{\boldsymbol{U}}}{2}\right)+\displaystyle \frac{({\boldsymbol{B}}\cdot {\rm{\nabla }}{\boldsymbol{U}})}{{B}^{2}}\cdot {\boldsymbol{B}}\pm \displaystyle \frac{{\boldsymbol{B}}\cdot {\rm{\nabla }}| {V}_{A}| }{| B| }.\end{eqnarray} \tag{ 17 }$$

The model Equations (11)–(13), (15), and (16) are coupled to our 3D MHD model, utilizing the inhomogeneous solar wind ${\boldsymbol{U}}$ and ${{\boldsymbol{V}}}_{A}$ obtained from the MHD simulation, i.e., we solve the temporal variation of eight MHD variables and six turbulence variables simultaneously in the coupled model.

We note that related approaches have been introduced by Usmanov et al. (2000, 2011, 2014) and Kryukov et al. (2012). However, these authors coupled the much simpler turbulent transport models of Zank et al. (1996) and Breech et al. (2008) to 3D MHD models of the heliosphere. The Usmanov et al. (2014) and Kryukov et al. (2012) 3D MHD models included the back reaction of the turbulence on the global MHD flow by including both heating and Reynolds-averaged turbulence terms in the MHD model. This will be an important next step in the further development of our approach.

The coupling of turbulence and heating by dissipation of the turbulence using a 1D Zank et al. (2012) model was investigated by Adhikari et al. (2015a), The intensities of inward and outward propagating turbulence show a large asymmetry in fast solar wind streams, whereas they are comparable in slow solar wind streams.

Although these are coupled with each other, in this study, we show results in a steady state with coupling in one direction only (no source term, H = 0 in Equation (4)).

The origin of the three turbulence intensities described in the previous section depend on the assumed inner boundary conditions of the simulation and the generation of turbulence by various sources within the numerical domain.

2.3.1. Boundary Condition

We consider the turbulent energy densities g and f to be associated with forward and backward propagating modes, respectively. The propagation direction in the radial direction depends on the sign of the IMF radial component. If we used actual directional inhomogeneous magnetic field data, the boundary condition would have to be specified taking into account the IMF direction. However, the development of such a technique is beyond of the scope of this paper, and we specify instead spherically symmetric values on the inner boundary, which are slightly modified from those used in Adhikari et al. (2015a), as listed in Table 1.

Table 1.  Turbulence Variables on the Inner Boundary (Same as in Adhikari et al. 2015a)

Name	symbol	values	unit
Strength of outward propagating mode	g0	14,000	(km s−1)2
Strength of inward propagating mode	f0	268	(km s−1)2
Strength of Residual energy	ED0	−57.07	(km s−1)2
Correlation length of outward propagating mode	${\lambda }_{-0}$	$1.16\times {10}^{5}$	km
Correlation length of inward propagating mode	${\lambda }_{+0}$	$2.14\times {10}^{5}$	km
Correlation length of residual energy	${\lambda }_{D0}$	$3.05\times {10}^{6}$	km

Download table as:  ASCIITypeset image

To examine if the boundary conditions are inappropriate, we performed simulations with the turbulence boundary conditions with forward and backward propagating modes swapped (g, f, ${\lambda }_{\pm }$). The parameters for each case are tabulated in Table 2. We denote the cases where $g={g}_{0}$, $f={f}_{0}$, ${\lambda }_{+}={\lambda }_{+0}$, and ${\lambda }_{-}={\lambda }_{-0}$, as "proper" in Table 2. The swapped cases, where $g={f}_{0}$, $f={g}_{0}$, ${\lambda }_{+}={\lambda }_{-0}$, and ${\lambda }_{-}={\lambda }_{+0}$, are denoted as "opposite."

Table 2.  Parameter List of the Numerical Results

Case	Boundary	Solar wind		IMF	Turbulence
	condition	Structure	${V}_{\mathrm{SW}}$ [km s−1]	Structure	source term
Case 1	proper	1D	600	Monopole	Yes
Case 2	proper	1D	400	Monopole	Yes
Case 3	proper	1D	600	Monopole, $\times 1/2$	Yes
Case 4a	proper	1D	600	Monopole	No
Case 4b	proper	1D	600	Monopole	Yes, $\times 2$
Case 5	opposite	1D	600	Monopole	Yes
Case 6	proper	Bimodal		Monopole	Yes
Case 7a	proper	Bimodal		Dipole	Yes
Case 7b	opposite	Bimodal		Dipole	Yes
Case 8a	proper	Bimodal		Dipole	Yes, $\times 2$
Case 8b	opposite	Bimodal		Dipole	Yes, $\times 2$

Download table as:  ASCIITypeset image

2.3.2. Source Terms

We applied only one form of turbulence source terms, which is that associated with the generation of shear at the interface of fast and slow solar wind streams (Coleman 1968; Roberts et al. 1987, 1992; Zank et al. 1996; Adhikari et al. 2015a). Roberts et al. (1987) investigated observationally the origin and evolution of fluctuations in the solar wind using the Helios and Voyager spacecraft data sets and observed rapid evolution of the normalized cross-helicity in the inner heliosphere, which is indicative of strong in situ nonlinear dynamical processes, i.e., the shear of the streams. Zank et al. (1996) modeled this process with a velocity shear source term in the form of $C{\rm{\Delta }}U\langle {z}^{2}\rangle /r$ though the generation of turbulence by stream shear should be independent of intensities of the background turbulence. Therefore, we apply the source terms for solar wind turbulence, which are extensively studied elsewhere (Zank et al. 2017). We constructed dimensionally new source terms parameterized based on (velocity)³/length ∼(velocity)²/time. The way we chose the (velocity)³ term was to combine the two characteristic speeds, viz., the Alfvén speed and the speed of fluctuations. We find that the combination ${V}_{A}^{2}{\rm{\Delta }}U$ seems to work best where ${\rm{\Delta }}U\sim 5\mbox{--}10$ km s⁻¹ is of the order of the characteristic speed of velocity fluctuations. We could have written the source term in e.g., Equation (11) as Constant $\times {\rm{\Delta }}{{UV}}_{A}^{2}{r}_{0}/{r}^{2}={C}_{+}{V}_{A}^{2}{r}_{0}/{r}^{2}$. The dimension of ${C}_{+}$ and ${C}_{-}$ is that of velocity. We use the following modified 1D model

$$\begin{eqnarray}&&{S}^{+}=\displaystyle \frac{{C}_{+}{V}_{A0}^{2}{r}_{0}}{{r}^{2}};\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&{S}^{-}=\displaystyle \frac{{C}_{-}{V}_{A0}^{2}{r}_{0}}{{r}^{2}};\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&{S}_{D}=-\displaystyle \frac{{S}^{+}+{S}^{-}}{4},\end{eqnarray} \tag{ 20 }$$

where ${V}_{A0}=174.4$ km s⁻¹, ${r}_{0}=62.35{R}_{\odot }$, ${C}_{+}=2.00{V}_{0}$ km s⁻¹, and ${C}_{-}=0.80{V}_{0}$ km s⁻¹. These forms were chosen to ensure consistency of the model and observed values of the various turbulence quantities. To examine the effect of the source terms, we consider cases in which the sources are set to ${C}_{\pm }=0$ or two times stronger than shown above. These cases are denoted as "No" or "$\times 2$" in the last column of Table 2, respectively.

Since we consider only solar minimum conditions, source terms associated with the generation of turbulence by shock waves (Zank et al. 1996; Adhikari et al. 2015a) are neglected. Similarly since the simulated domain is restricted to within 6 au, pickup ion driven turbulence (Lee & Ip 1987; Williams & Zank 1994; Zank et al. 1996) is neglected too.

First, for the Basic Case (Case 1), we assumed a spherically symmetric solar wind flow of 600 km s⁻¹ at the inner boundary (60 ${R}_{\odot }$) and a monopole IMF with magnetic field strength is 0.55 G at source surface 2.5 ${R}_{\odot }$. The radial field in the magnetic field vector becomes B_r = 7.43 nT at 1 au, which is a comparable value to that observed in situ (Adhikari et al. 2015a).

Figure 1 shows two-dimensional distributions of velocity (panel (a)), magnetic field (panel (b)), magnetic field alignment to the radial direction (panel (c)), and turbulence variables (panel (d)–(i)) in the meridional plane in the HGI coordinate. In spite of the spherically symmetric boundary conditions, the steady state shows latitudinal distributions resulting from the rotation of the Sun (the inner boundary), i.e., forming the Parker spiral IMF (Parker 1958; see Figure 1(c)). V_r and B_r have spherically symmetric distributions, while the magnetic field alignment $| {B}_{r}| /| B|$ has a cylindrically symmetric distribution. Because of the inhomogeneous distribution of the background solar wind and IMF, the distributions of turbulence intensities and correlation lengths have latitudinal dependences (Figures 1(d)–(i)).

Zoom In Zoom Out Reset image size

For a quantitative comparison between the four different latitudes, we use the 1D radial plots of the background solar wind, IMF, and the turbulent variables in Figure 2. The different colors in each radial plot correspond to different latitudes (0° (black), 30° (blue), 60° (green), and 90° (red)). Figures 2(a) and (b) show the solar wind speed and IMF along the four different latitude lines, respectively. The profiles of the solar wind speed and the IMF radial component are identical due to the spherically symmetric distribution seen in Figure 1, while only the IMF azimuthal component has different distribution.

Zoom In Zoom Out Reset image size

Figure 2(e) shows radial plots of the three turbulence strengths: the solid lines denote g, the dashed lines f, and the dotted lines $-{E}_{D}$. For the outward propagating intensity g, the distributions are similar and small differences can be seen beyond 2 au. In contrast, the profiles of the other two intensities show a stronger dependence with latitude, with the lower latitude profiles showing a greater decrease at larger distances. From Figure 1(d), since the contours of the distribution f seem to be more vertically oblate with distance from the Sun, the distribution will approach a cylindrical distribution, and the transport of the inward propagating intensity f is likely to depend on the direction of IMF (Figure 1 (c)) in these regions.

From the turbulence transport Equations (11)–(13), (15), and (16), the dependence on the magnetic field alignment $| {B}_{r}| /| B|$ appears only in the mixing term as shown in Equation (17). For large distances, Equation (17) can be reduced as

$$\begin{eqnarray}&&{M}^{\pm }\sim -{\rm{\nabla }}\cdot \left(\displaystyle \frac{{\boldsymbol{U}}}{2}\right)+\displaystyle \frac{{B}_{r}^{2}}{{B}^{2}}\displaystyle \frac{\partial {U}_{r}}{\partial r}+\displaystyle \frac{{B}_{\theta }^{2}+{B}_{\phi }^{2}}{{B}^{2}}\displaystyle \frac{{U}_{r}}{r}\sim -\displaystyle \frac{{B}_{r}^{2}}{{B}^{2}}\displaystyle \frac{{U}_{r}}{r},\end{eqnarray} \tag{ 21 }$$

assuming that ${\boldsymbol{U}}\sim {U}_{r}{{\boldsymbol{e}}}_{r}$ and ${V}_{A}\ll {U}_{r}$. This indicates that mixing of the inward and outward propagating modes and the residual energy can be the relatively dominant factor in determining the backward propagating mode intensity in the distant region. The terms that include the mixing in Equations (11) and (12) are positive because E_D is negative. In the low latitude region, the term approaches 0 and hence the absence of mixing results in lower intensities.

Figures 2(c) and (d) show the normalized cross-helicity (solid), residual energy (dashed), and Alfvén ratio (Bruno et al. 2007), respectively,

$$\begin{eqnarray}&&{\sigma }_{C}=\displaystyle \frac{g-f}{g+f}\times \displaystyle \frac{{B}_{r}}{| {B}_{r}| };\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&{\sigma }_{D}=\displaystyle \frac{2{E}_{D}}{g+f};\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}&&{r}_{A}=\displaystyle \frac{g+f+2{E}_{D}}{g+f-2{E}_{D}}.\end{eqnarray} \tag{ 24 }$$

Corresponding to the latitudinal dependence seen in Figure 2(e), these three parameters reflect the same latitudinal dependence.

Notice too that the decay of the normalized residual energy σ_D indicates that the turbulence is becoming increasingly dominated by the magnetic field fluctuations (see Adhikari et al. 2015a, 2015b) rather than the kinetic energy. In the inner region, the dominance of the outward propagating turbulence intensity relative to the other two intensities (Figure 2(e)) reflects that the normalized cross-helicity is nearly equal to +1 and that the normalized residual energy is ∼0. The normalized residual energy approaches approximately −0.5 as the intensity of the residual energy becomes half of the average of the inward and outward propagating turbulence in the outer regions (Figure 2(e)). This is reflected in the Alfvén ratio plot r_A, which approaches ∼0.3 (Figure 2(d)). The normalized cross-helicity approaches ∼0.5 as the intensity of the inward propagating turbulence becomes half that of the outward propagating turbulence. The reductions in the intensities of the inward propagating turbulence and the residual energy in lower latitudes results in the difference in the radial profiles of ${\sigma }_{C}$, ${\sigma }_{D}$, and r_A (Figures 2(c) and (d)).

Figure 2(f) shows radial plots of the three turbulence correlation lengths: the solid lines denote ${\lambda }_{-}$, the dashed lines ${\lambda }_{+}$, and the dotted lines ${\lambda }_{D}$. The distributions of all the correlation lengths show a similar latitudinal dependence.

The line in each panel of Figure 3 show the radial profiles in the equatorial plane (solid lines) and over the poles (dotted lines) for all the variables. Observed values (Adhikari et al. 2015a) for the three strengths are overlaid using the cross symbols in Panels (d), (e), (g), (h), and (i) for reference. A word of caution is appropriate about the observational data overplotted on the panels of the radial plots. As discussed at length in Adhikari et al. (2015a, 2015b), the observations combine Helios data from 0.29 au to 0.9 au and Ulysses data 1.3–4.8 au. The data, taken during solar minimum conditions, is nonetheless from different decades, used different instrumentation, and is a combination of in-ecliptic and out-of-ecliptic solar wind data. During the periods of observation, the IMF was approximately radial. Despite these caveats, the observed turbulence properties appear to provide a reasonable quantification of the evolution of various turbulence parameters from 0.29 to ∼5 au, and is in reasonable in accordance with the theoretical model of Adhikari et al. (2015a, 2015b).

Zoom In Zoom Out Reset image size

Case 2 differs from Case 1 in that the solar wind speed is reduced to 400 km s⁻¹ from the 600 km s⁻¹ of Case 1. Since the behavior of the global distributions is similar to those in Case 1, only radial profiles along the equatorial line and over the poles are compared, as shown with red lines in Figure 3 to illustrate the latitudinal dependence.

The main differences are in the speed (Figure 3(a)) and hence in the IMF because of the Parker spiral, as shown in Figure 3(c). The azimuthal component in Case 2 at the equator is 1.5 times larger than that in Case 1.

A comparison of the turbulence intensities is shown in Figures 3(g)–(i). The profiles of the three turbulence intensities along the pole (dotted lines in the panels) are almost identical. This means that the advection of the turbulence in the completely radial magnetic field does not introduce any differences.

In the equatorial plane, the inward propagating mode intensity is smaller than that of Case 1, intensity is larger than that of Case 1, while that of outward propagating modes is almost identical. These different trends can be understood from the terms related to the Alfvén velocity in Equation (17), i.e., the term including b in Equation (14), since the relative amplitude of ${{\boldsymbol{V}}}_{A}$ and its gradient relative to ${\boldsymbol{U}}$ are larger in Case 2 than in Case 1. Therefore, the difference between inward and outward propagating modes, i.e., cross-helicity, tends to be amplified (see Figure 3(d)) with the heliocentric distance. This is a consequence of the difference in the Parker spiral of IMF between the two cases, since the magnetic field alignment $| {B}_{r}| /| B|$ decreases faster in Case 2. This effect can suppress the mixing of turbulence in the inner region compared to Case 1.

From observations, it is known that for turbulence in the slow solar wind is balanced, having the same amplitude for inward and outward propagating modes. Thus, cross-helicity should be around zero. However, the results of our model with slow wind show a non-zero cross-helicity distribution. This can be caused by the use of inappropriate boundary values for the three turbulence intensities and corresponding correlation lengths. Although an improvement of the values for the slow solar wind is needed for more realistic studies, this is beyond the scope of the current study.

Case 3 differs from Case 1 in that a different IMF strength is assumed, where the strength is reduced to half that of Case 1. The comparison of the radial profiles to Case 1 is shown in Figure 4. Because the solar wind speed is the same as in Case 1, the shape of the Parker spiral is identical (the ratio between the azimuthal and radial components seen in Figure 4(c)). Hence the mixing terms related to the magnetic field alignment are the same while the terms related to the amplitude of Alfvén velocity (namely ${\rm{\nabla }}\cdot {{\boldsymbol{V}}}_{A}$) are halved.

Zoom In Zoom Out Reset image size

The difference in the turbulence intensities is slight, though they can be seen at both the equator and the pole. This can be understood as the effect of an Alfvén velocity that is relatively small in the high-speed solar wind (Case 3).

Cases 4a (Figure 5) and 4b (Figure 6) are comparative cases to show that the source term (Zank et al. 1996; Adhikari et al. 2015a) is necessary to reproduce quantitatively the radial profiles of both inward and outward propagating turbulence. Case 4a neglects the source terms and Case 4b assumes source terms two times larger than those of Case 1, while the solar wind and IMF are identical to those in Case 1.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

All three turbulence intensities in Case 4a (Figure 5) are reduced compared to those of Case 1 and the reduction in amplitudes is approximately similar (> a factor of five smaller) at both the equator and the pole. The radial dependence of the intensities also differs from those observed. Corresponding to the reduction in intensity, all the correlation lengths are ∼5 times longer. Because the relative amplitudes are different in their intensities, the difference in normalized cross-helicity ${\sigma }_{C}$ and residual energy ${\sigma }_{D}$ is not so different compared to the difference in turbulent intensities.

In contrast to Case 4a, all the turbulence intensities in Case 4b (Figure 6) are increased by a factor of about 1.5 compared to those of Case 1. Corresponding to the increase in the intensities, the correlation lengths are reduced by the same factor. That increase by a factor of ∼1.5, which is smaller than the factor of the source terms, can be understood as the increase in the nonlinear dissipation caused by the shorter correlation lengths. The radial profiles of the turbulence intensities in Case 4b show good agreement with the observed radial profiles (Figures 6(g), (h), and (i)).

In both Cases 4a and 4b, the increase (decrease) in the turbulence intensities is inversely related to the corresponding correlation lengths. This can be interpreted as follows. In our numerical models, the time variation of the correlation variables ${L}^{\pm }$ and L^D are solved using Equations (15) and (16). Because the background solar wind and IMF and the boundary conditions of the correlation variables ${L}^{\pm }$ and L^D are identical to those in Case 1, the distributions for ${L}^{\pm }$ and L^D are also identical. In contrast to ${L}^{\pm }$ and L^D, the turbulence intensities are different reflecting the difference in the source terms. Because the correlation lengths are obtained by dividing the correlation variables by the intensities (i.e., ${\lambda }^{+}={L}^{+}/f$, ${\lambda }^{-}={L}^{-}/g$, ${\lambda }^{D}={L}^{D}/{E}_{D}$) , they are inversely related to the turbulence intensities.

Case 5 (Figure 7) is one for which the inner boundary conditions of the inward and outward propagating modes (f and g) and corresponding correlation variables (${L}^{\pm }$) are swapped as well as the source term (S_±). Because of the dominance of the inward propagating mode, the cross-helicity has negative values. The increase of residual energy in vicinity of the inner boundary is suppressed at both the equator and the pole. This indicates that the energy density in the inward and outward propagating modes in Case 5 is more inefficiently converted to the residual energy than in Case 1, i.e., more energy is dissipated through the nonlinear terms. The turbulence intensity of each dominant mode (f in Case 5) at r = 6 au is smaller by a factor of $\gt 1.5$ than g in Case 1 in both the equatorial plane and over the pole.

Zoom In Zoom Out Reset image size

The radial intensity profiles of inward (f) and outward (g) propagating modes and the residual energy (E_D) finally converge to the same value of ∼100 km² s⁻² at the poles (red dotted lines in Figures 7(g)–(i)) at r = 6 au. This indicates that the turbulence over the pole is similar isotopically to that observed in the slow solar wind. In contrast to the pole, in the equatorial plane, the intensity of the dominant (inward propagating) mode f is a factor of 1.5 times higher than that of the outward propagating mode g (Figure 7(h)). This indicates that the radial dependence of the turbulence intensities is strongly dependent on the values specified at the inner boundary. Around lower latitudes, the mixing between f and g is less efficient and their dissipation is stronger than at the pole.

The result shows the importance of the choice in boundary values for turbulence intensities and indicates that the highly asymmetric nature of turbulence in the fast solar wind is very likely generated in the inner heliosphere, i.e., within the inner boundary of the simulation in this study.

For Case 6, a tilted bimodal solar wind is introduced with the identical IMF boundary used condition in Case 1. Figures 8 and 9 show distributions of the background solar wind and turbulence variables in the meridional and equatorial planes, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In spite of the spherically symmetric radial magnetic field on the inner boundary, the magnetic field has a 3D distribution (Figures 8(b), (c), 9(b), and (c)). The radial distributions of the variables along the solid lines in Figures 8(a) and 9(a) are displayed with solid lines in each panel of Figure 10. The smallest radial flow velocity along the solid line is located around 2 au (Figure 10(a)). The magnetic field strengths of both components are enhanced behind this region and reduced ahead of it (Figures 10(b) and (c)). Before the slowest radial flow velocity, the radial gradient of the flow speed is negative, and therefore the divergence of the flow is negative, i.e., corresponding to compression (Equations (1) and (3)). Beyond the radial flow minimum, there is a positive gradient, i.e., expansion. A slight enhancement of the IMF results from the expansion and compression of the flow due to its interaction with the solar wind distribution (Figures 8(a) and 9(a)). The slightly enhanced field region ahead of fast solar wind streams is called a stream interaction region, or co-rotating interaction region (CIR) in general.

Zoom In Zoom Out Reset image size

Corresponding to the IMF enhancement, the turbulence intensities of the inwardly propagating modes and the residual energy are enhanced and reduced at the same locations when considered along the equator line (the black solid lines). In contrast to the case without the flow inhomogeneity (Case 1, Figure 2), the difference in inward propagating mode intensity between on the equator (black) and 30° (blue) is within a factor of 0.1 in logarithmic scale. For Case 6, the enhancement is more than a factor of 0.2 in the equatorial plane behind the velocity minimum. In addition to the divergence and convergence effects, the solar wind speed itself affects the amplitude of the turbulence intensities as shown in Section 3.2. The slow speed effect and the expansion can explain the strong reduction in the turbulence intensities ahead of the radial flow velocity minimum (Figures 10(g) and (i)).

A consequence of the decrease in the intensities of the inward propagating mode is the rapid decrease of the cross-helicity (Figure 10(d)). In contrast, ahead of the flow speed minimum, the cross-helicity increases. The decrease ahead of the radial flow minimum can be seen in the residual energy. This results in a rapid increase in the normalized residual energy (Figure 10(e)) and Alfvén ratio (Figure 10(f)).

The various enhancements discussed above (i.e., strong inhomogeneity of the solar wind speed) can be seen mainly in the radial profiles in the equatorial plane shown in Figure 10. Because the tilt of the bipolar flow on the inner boundary is 30°, the slow wind streams located along the parallel plane are distributed in at latitudes below 30°. Spiky increases and decreases in the turbulence variables are prominent at the equator. That is interpreted as a consequence of shock formation occurring due to the nonlinear steepening of the global MHD structure.

For Case 7a, the magnetic field boundary condition is replaced by a tilted dipole that aligns with the bipolar flow specified in Case 6. The 2D distributions in the equatorial plane are shown in Figure 11 and radial plots along a line in the equatorial plane and over both poles are shown in Figure 12. Profiles for the same locations as in Case 6 are overplotted in Figure 12 with black lines.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

A comparison between Case 6 (Figure 9) and Case 7a (Figure 11) shows the effect of the inhomogeneous background IMF. Note that we specified a tilted dipole field whose strength is not spherically symmetric. The difference in the field strength can cause a slight difference in the solar wind speed profile (Figure 12(a)) compared to that of Case 6. The sign of magnetic field changes at the line perpendicular to the axis of the dipole, at which the solar wind is slowest in this model. The field strength near the magnetic neutral line is close to 0. As the solar wind propagates, the compressions and expansions described in Section 4.1 developed similarly to Case 6. However, there is no enhancement around the converging area in Case 7a. This can be understood since the compression toward the magnetic neutral line (surface) can cause the enhancement of the weak field strength region and might cause magnetic reconnection. The flow due to the reconnection can contribute to a change in the solar wind speed profiles.

Because we did not take into account the effect of the magnetic field polarity on the assumed boundary values of the turbulence, over part of the inner hemisphere (north in Case 7a and south in Case 7b, which is mentioned in the next section), the turbulence boundary values and magnetic field polarity are inconsistent. The source terms also depend on the sign of the IMF. As discussed in Section 3.5, when inappropriate turbulent intensities are specified on the boundary (Case 5) the intensity of the dominant mode decrease more rapidly than the case with appropriate boundary conditions (Case 1). Along a line on the equatorial plane in Case 7a, the polarity of the IMF is opposite before and after $r\sim 2$ au (Figure 12(b)). The IMF in the inner region ($r\lt 2$ au) is the same polarity as in Case 6 and that in the outer region is opposite. Comparing the turbulence intensities in the outer region in Case 6 and Case 7a shows that the influence of the boundary condition appears as a decrease in the dominant propagating mode g. In this region, g seems to be a factor of 1.3 smaller than that in Case 6. This is consistent with the results of Case 5 shown in Section 3.5. The other two intensities do not change as much. The correlation lengths in this region are shorter for the outward propagating mode intensity f here but the other two are not different from those of Case 6.