DIRECT AND INVERSE CASCADES IN THE ACCELERATION REGION OF THE FAST SOLAR WIND

The fast solar wind emanating from coronal holes is believed to be driven by Alfvén waves that propagate outward along the open field lines (e.g., Parker 1965; Heinemann & Olbert 1980; Velli 1993; Dmitruk et al. 2002; Suzuki & Inutsuka 2005; Cranmer et al. 2015). In situ observations in the heliosphere indicate that the waves are in a turbulent state with a broad spectrum of wavenumbers and frequencies (e.g., Coleman 1968; Belcher & Davis 1971; Hollweg 1986; Matthaeus et al. 1990; Bale et al. 2005; Borovsky 2012). Alfvén waves have also been detected by remote sensing of the solar atmosphere (De Pontieu et al. 2007; Tomczyk et al. 2007; Tomczyk & McIntosh 2009; Tian et al. 2011, 2014; Threlfall et al. 2013; Liu et al. 2015; Morton et al. 2015). Alfvén waves are a prime candidate for heating and accelerating the fast wind because they have the ability to transport energy over large distances in the corona (Barnes 1966; Belcher 1971; Hollweg 1973; Jacques 1977; Velli et al. 1989; Marsch & Tu 1997; Matthaeus et al. 1999; Suzuki & Inutsuka 2006; Cranmer et al. 2007; Chandran et al. 2011). The slow solar wind may also be driven by Alfvén waves (Oran et al. 2015).

A turbulent cascade has long been considered a promising mechanism for dissipation of Alfvén waves in the solar wind (e.g., Hollweg et al. 1982; Hollweg 1986; Velli et al. 1989). Nonlinear interactions between counterpropagating Alfvén waves are known to produce turbulence (Iroshnikov 1963; Kraichnan 1965; Shebalin et al. 1983; Goldreich & Sridhar 1995, 1997; Bhattacharjee & Ng 2001; Maron & Goldreich 2001; Cho et al. 2002). The turbulence can be described in terms of Elsasser variables, ${{\boldsymbol{z}}}_{\pm }={{\boldsymbol{v}}}_{1}\mp {{\boldsymbol{B}}}_{1}/\sqrt{4\pi {\rho }_{0}}$, where ${{\boldsymbol{B}}}_{1}$ and ${{\boldsymbol{v}}}_{1}$ are the magnetic and velocity fluctuations of the waves, and ${\rho }_{0}$ is the mean plasma density (Elsasser 1950). The ${{\boldsymbol{z}}}_{+}$ and ${{\boldsymbol{z}}}_{-}$ waves are linearly coupled because of radial gradients in plasma density and magnetic field strength (Heinemann & Olbert 1980; Velli 1993; Hollweg & Isenberg 2007), a process often described as wave "reflection." Hence, the dominant ${{\boldsymbol{z}}}_{+}$ waves produce a lower level of ${{\boldsymbol{z}}}_{-}$ waves, which we refer to as the "minority" waves. In the solar wind, the dominant waves are outward-propagating, but the minority waves can have both inward- and outward-propagating components (e.g., Velli et al. 1989; Verdini et al. 2009; Perez & Chandran 2013). Matthaeus et al. (1999) proposed that the corona may be heated by Alfvén wave turbulence driven by nonlinear interactions between the ${{\boldsymbol{z}}}_{+}$ and ${{\boldsymbol{z}}}_{-}$ waves. Detailed models of the solar wind based on these ideas since have been developed (e.g., Cranmer & van Ballegooijen 2005; Cranmer et al. 2007; Verdini & Velli 2007; Verdini et al. 2010; Chandran et al. 2011; Oran et al. 2013; Sokolov et al. 2013; Lionello et al. 2014; van der Holst et al. 2014). Woolsey & Cranmer (2015) have shown that the turbulent heating varies strongly in time and space. However, it should be kept in mind that the turbulent cascade is not the only mechanism for dissipating Alfvén waves in the corona. Nonlinear coupling between Alfvén and compressive waves may also play an important role (Kudoh & Shibata 1999; Moriyasu et al. 2004; Chandran 2005; Suzuki & Inutsuka 2005, 2006; Matsumoto & Shibata 2010; Cranmer & van Ballegooijen 2012).

Models of solar wind turbulence have been developed by several authors. One approach is to use the "shell" model, which simplifies the nonlinear interactions by reducing the number of wave modes that are allowed to interact (Velli et al. 1989; Buchlin & Velli 2007). This model has the great advantange that very high perpendicular wavenumbers can be reached. Verdini et al. (2009) and Verdini et al. (2012) used the shell model to study the formation and evolution of a turbulent spectrum of Alfvén waves produced by linear and nonlinear wave couplings. Another approach to turbulence modeling is to perform direct numerical simulations using the reduced magnetohydrodynamic (RMHD) equations (e.g., Dmitruk et al. 2002; Dmitruk & Matthaeus 2003; Perez & Chandran 2013). These equations include the nonlinear couplings that lead to turbulence, but omit other effects such as the coupling between Alfvén and compressive waves. Dmitruk et al. (2002) argued that reflection-driven turbulence provides a robust heating mechanism that can explain the observed temperatures in the region below $2\,{R}_{\odot }$. Dmitruk & Matthaeus (2003) showed that a high dissipation efficiency can be obtained when the nonlinear timescale of the turbulence is less than the Alfvén crossing time. Perez & Chandran (2013) were the first to include the effects of the solar wind flow on wave propagation in the RMHD model. They found that up to one-third of the wave energy launched at the coronal base is dissipated in the corona below the Alfvén critical point, and another third goes into doing work on the solar wind outflow.

In a previous paper (van Ballegooijen & Asgari-Targhi 2016, hereafter paper I), we presented RMHD simulations of Alfvén wave turbulence for the fast solar wind emanating from a polar coronal hole. This modeling is an extension of our earlier work on turbulence in coronal loops (e.g., van Ballegooijen et al. 2011, 2014; Asgari-Targhi & van Ballegooijen 2012; Asgari-Targhi et al. 2013, 2014). Paper I includes the effects of the solar wind flow on Alfvén-wave propagation. Two models of the solar wind are considered. In the first model, the plasma density and Alfvén speed vary smoothly with height along the modeled flux tube. We find that for this "smooth" model the linear wave coupling is relatively weak, producing only a low level of minority waves. Therefore, the energy dissipation rate of the turbulence is insufficient to maintain the temperature of the background atmosphere. We also present a second model with additional, random density variations that approximate the effects of compressive MHD waves in the solar wind. We find that such spatial variations in density can significantly enhance the minority waves and thereby the turbulent dissipation rates.

The results of paper I led us to conclude that interactions between Alfvén and compressive waves may play an important role in the turbulent heating of the fast solar wind. However, this conclusion is somewhat premature because the reason(s) for the low dissipation rates are not yet well understood. In particular, we do not know whether the low rates are a real physical effect or a numerical artifact. In the present paper we further improve our numerical model, and we compute for the first time the cascade rates ${\epsilon }_{\pm }$ as functions of perpendicular wavenumber and radial distance from the Sun center. Cascade rates have also been measured in the solar wind, using third-order structure functions (e.g., Stawarz et al. 2009; Coburn et al. 2014). We find that an analysis of the cascade rates can shed light on the question why "smooth" solar wind models produce relatively low dissipation rates.

The RMHD equations describe the nonlinear interactions between counterpropagating Alfvén waves (e.g., Strauss 1976, 1997). The waves are assumed to propagate in a medium with a fixed background magnetic field ${{\boldsymbol{B}}}_{0}({\boldsymbol{r}})$ and density ${\rho }_{0}({\boldsymbol{r}})$, where ${\boldsymbol{r}}$ denotes the position. Specifically, we consider a thin, open magnetic flux tube extending radially from the Sun inside a coronal hole, so the field strength B₀ and density ${\rho }_{0}$ are functions of radial distance r from the Sun center. The waves can be described in terms of Elsasser variables, ${{\boldsymbol{z}}}_{\pm }(x,y,r,t)\equiv {{\boldsymbol{v}}}_{1}\mp {{\boldsymbol{B}}}_{1}/\sqrt{4\pi {\rho }_{0}}$, where ${{\boldsymbol{B}}}_{1}$ and ${{\boldsymbol{v}}}_{1}$ are the magnetic and velocity fluctuations of the waves (Elsasser 1950), x and y are the coordinates perpendicular to the flux tube axis, and t is the time. Oughton et al. (2001) and Dmitruk & Matthaeus (2003) were the first to present RMHD simulations for such an open flux tube, and they found that reflection-driven turbulence can be maintained in this environment despite the fact that the waves can escape into the heliosphere. Perez & Chandran (2013) included the effects of the solar wind outflow on the waves and presented a variety of RMHD models with different perpendicular correlation lengths and correlation times of the imposed footpoint motions. In paper I we investigated whether RMHD models of wave turbulence can explain the observed heating of the fast solar wind. In the present work, we continue this investigation with a more detailed analysis of the cascade processes. We shall refer to the ${{\boldsymbol{z}}}_{+}$ and ${{\boldsymbol{z}}}_{-}$ waves as the dominant and minority waves, respectively.

The Elsasser variables are nearly incompressible velocity fields and can be written as ${{\boldsymbol{z}}}_{\pm }={{\rm{\nabla }}}_{\perp }{f}_{\pm }\times \,{\hat{{\boldsymbol{B}}}}_{0}$, where ${{\rm{\nabla }}}_{\perp }$ is the spatial derivative in the x and y directions, ${f}_{\pm }({\boldsymbol{r}},t)$ are the velocity stream functions, and ${\hat{{\boldsymbol{B}}}}_{0}(x,y,r)$ is the unit vector along the background field. The RMHD equations can be written as

$$\begin{eqnarray}\displaystyle \frac{\partial {\omega }_{\pm }}{\partial t} & = & -({u}_{0}\pm {v}_{{\rm{A}}})\displaystyle \frac{\partial {\omega }_{\pm }}{\partial r}\\ & & +\displaystyle \frac{1}{2}\left(\displaystyle \frac{{{dv}}_{{\rm{A}}}}{{dr}}\pm \displaystyle \frac{{u}_{0}}{2{H}_{\rho }}\right)({\omega }_{+}-{\omega }_{-})+\displaystyle \frac{{u}_{0}}{2{H}_{{\rm{B}}}}({\omega }_{+}+{\omega }_{-})\\ & & -\displaystyle \frac{1}{2}[{\omega }_{+},{f}_{-}]-\displaystyle \frac{1}{2}[{\omega }_{-},{f}_{+}]\pm {{\rm{\nabla }}}_{\perp }^{2}\left(\displaystyle \frac{1}{2}[{f}_{+},\,{f}_{-}]\right),\end{eqnarray} \tag{ 1 }$$

where ${\omega }_{\pm }\equiv -{{\rm{\nabla }}}_{\perp }^{2}{f}_{\pm }$ are the vorticities associated with the dominant and minority waves, ${u}_{0}(r)$ is the outflow velocity of the wind, ${v}_{{\rm{A}}}(r)$ is the Alfvén speed, ${H}_{B}(r)\equiv {B}_{0}/({{dB}}_{0}/{dr})$ is the magnetic scale length, and ${H}_{\rho }(r)\equiv {\rho }_{0}/(d{\rho }_{0}/{dr})$ is the density scale length. Equation (1) can be derived from the expressions given in Perez & Chandran (2013) and paper I. The first term on the right-hand side of this equation describes the effects of wave propagation, and the second and third terms describe the linear couplings between the dominant and minority waves. The bracket operator [⋯ , ⋯] is defined by

$$\begin{eqnarray}&&b(x,y,r,t)=[f,g]\equiv \displaystyle \frac{\partial f}{\partial x}\displaystyle \frac{\partial g}{\partial y}-\displaystyle \frac{\partial f}{\partial y}\displaystyle \frac{\partial g}{\partial x},\end{eqnarray} \tag{ 2 }$$

where $f(x,y,r,t)$ and $g(x,y,r,t)$ are two arbitrary functions. All nonlinearities of the RMHD model are contained within such bracket terms. In Equation (1) we omit the dissipative terms, which will be described in more detail below.

Most RMHD models use a spectral method in which all functions of x and y are written in terms of a set of normalized basis functions ${\tilde{F}}_{k}(\tilde{x},\tilde{y})$. Here $\tilde{x}$ and $\tilde{y}$ are dimensionless coordinates, and index k is in the range $k=1,\ldots ,{k}_{\max }$, where ${k}_{\max }$ is the total number of modes. Then an arbitrary function $f(x,y,r,t)$ can be written as

$$\begin{eqnarray}&&f(x,y,r,t)=\displaystyle \sum _{k}\,{f}_{k}(r,t){\tilde{F}}_{k}(\tilde{x},\tilde{y}),\end{eqnarray} \tag{ 3 }$$

where ${f}_{k}(r,t)$ is the amplitude of the mode with index k. The basis functions depend on the dimensionless perpendicular coordinates $\tilde{x}\equiv x/R(r)$ and $\tilde{y}\equiv y/R(r)$, where R(r) is the radius of the cross section of the flux tube. In paper I we assumed a circular cross section, ${\tilde{x}}^{2}+{\tilde{y}}^{2}\leqslant 1$, but in the present work we follow Perez & Chandran (2013) by assuming a square cross section, $-1\leqslant \tilde{x}\leqslant +1$ and $-1\leqslant \tilde{y}\leqslant +1$, and we use periodic boundary conditions on this square domain. Then the basis functions are products of $\tilde{x}$- and $\tilde{y}$-dependent parts, each of which are sine or cosine functions with periods ${\rm{\Delta }}\tilde{x}={\rm{\Delta }}\tilde{y}=2$. In this case, Equation (3) is essentially the Fourier transform written in a compact form. The width of the computational domain in dimensional units is ${\rm{\Delta }}x={\rm{\Delta }}y=2R(r)$, which increases with radial distance r from the Sun center. The basis functions have well-defined dimensionless perpendicular wavenumbers ${a}_{x,k}=\pi {n}_{x,k}$ and ${a}_{y,k}=\pi {n}_{y,k}$, where ${n}_{x,k}$ and ${n}_{y,k}$ are integers. The total dimensionless wavenumber is ${a}_{k}\equiv \sqrt{{a}_{x,k}^{2}+{a}_{y,k}^{2}}$, and the actual wavenumber in physical units is ${k}_{\perp }={a}_{k}/R(r)$. Inserting Equation (3) into Equation (2), we find the following for the mode amplitudes of the function $b(x,y,r,t)$:

$$\begin{eqnarray}&&{b}_{k}(r,t)=\displaystyle \frac{1}{{R}^{2}(r)}\displaystyle \sum _{j}\,\displaystyle \sum _{i}{M}_{{kji}}{f}_{j}(r,t){g}_{i}(r,t),\end{eqnarray} \tag{ 4 }$$

where ${f}_{j}(r,t)$ and ${g}_{i}(r,t)$ are the mode amplitudes of the arbitrary functions $f(x,y,r,t)$ and $g(x,y,r,t)$, and M_kji is a sparse, dimensionless matrix describing the nonlinear coupling between certain mode triples $(i,j,k)$. For the present case of a square cross section, we obtain

$$\begin{eqnarray}&&{M}_{{kji}}=\displaystyle \frac{1}{4}{\int }_{-1}^{+1}{\int }_{-1}^{+1}{\tilde{F}}_{k}(\tilde{x},\tilde{y})\left(\displaystyle \frac{\partial {\tilde{F}}_{j}}{\partial \tilde{x}}\displaystyle \frac{\partial {\tilde{F}}_{i}}{\partial \tilde{y}}-\displaystyle \frac{\partial {\tilde{F}}_{j}}{\partial \tilde{y}}\displaystyle \frac{\partial {\tilde{F}}_{i}}{\partial \tilde{x}}\right)\,d\tilde{x}\,d\tilde{y}.\end{eqnarray} \tag{ 5 }$$

In general, the details of the M_kji matrix depend on whether the flux tube has a circular or square cross section and on the type of boundary condition used, but the matrix is always fully antisymmetric in its indices, as was shown for the circular case in Appendix B of van Ballegooijen et al. (2011). Using Equation (3), the RMHD equations can be written as

$$\begin{eqnarray}\displaystyle \frac{\partial {\omega }_{\pm ,k}}{\partial t} & = & -({u}_{0}\pm {v}_{{\rm{A}}})\displaystyle \frac{\partial {\omega }_{\pm ,k}}{\partial r}\\ & & +\displaystyle \frac{1}{2}\left(\displaystyle \frac{{{dv}}_{{\rm{A}}}}{{dr}}\pm \displaystyle \frac{{u}_{0}}{2{H}_{\rho }}\right)({\omega }_{+,k}-{\omega }_{-,k})\\ & & +\displaystyle \frac{{u}_{0}}{2{H}_{{\rm{B}}}}({\omega }_{+,k}+{\omega }_{-,k})\\ & & +\displaystyle \frac{1}{2{R}^{4}}\displaystyle \sum _{j}\displaystyle \sum _{i}{M}_{{kji}}({a}_{i}^{2}-{a}_{j}^{2}-{a}_{k}^{2}){f}_{\pm ,j}{f}_{\mp ,i}\\ & & -{\nu }_{\pm ,k}\,{\omega }_{\pm ,k},\end{eqnarray} \tag{ 6 }$$

where ${\nu }_{\pm ,k}$ are artificial damping rates for dominant and minority waves. The damping model will be described in more detail in Section 3.

Multiplying Equation (6) by $\tfrac{1}{2}{\rho }_{0}{f}_{\pm ,k}$ and summing over modes, we obtain the wave energy equations for the dominant and minority waves:

$$\begin{eqnarray}&&\displaystyle \frac{\partial {U}_{\pm }}{\partial t}+{B}_{0}\displaystyle \frac{\partial }{\partial r}\left(\displaystyle \frac{{F}_{\pm }}{{B}_{0}}\right)=-{u}_{0}{D}_{\pm }\mp \displaystyle \frac{1}{2}\displaystyle \frac{{{dv}}_{{\rm{A}}}}{{dr}}{U}_{{\rm{R}}}-{Q}_{\pm },\end{eqnarray} \tag{ 7 }$$

where

$$\begin{eqnarray}&&{U}_{\pm }=\displaystyle \frac{{\rho }_{0}}{4{R}^{2}}\displaystyle \sum _{k}{a}_{k}^{2}{f}_{\pm ,k}^{2},\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{U}_{{\rm{R}}}=\displaystyle \frac{{\rho }_{0}}{2{R}^{2}}\displaystyle \sum _{k}{a}_{k}^{2}{f}_{+,k}{f}_{-,k},\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{Q}_{\pm }=\displaystyle \frac{{\rho }_{0}}{2{R}^{2}}\displaystyle \sum _{k}{\nu }_{\pm ,k}{a}_{k}^{2}{f}_{\pm ,k}^{2},\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{F}_{\pm }=({u}_{0}\pm {v}_{{\rm{A}}}){U}_{\pm }+\displaystyle \frac{1}{2}{u}_{0}\left({U}_{\pm }-\displaystyle \frac{{U}_{{\rm{R}}}}{2}\right),\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{D}_{\pm }=-\displaystyle \frac{1}{2}\displaystyle \frac{\partial }{\partial r}\left({U}_{\pm }-\displaystyle \frac{{U}_{{\rm{R}}}}{2}\right)-\displaystyle \frac{{U}_{{\rm{R}}}}{2{H}_{B}},\end{eqnarray} \tag{ 12 }$$

and we use mass conservation (${\rho }_{0}{u}_{0}/{B}_{0}$ = constant). Here ${U}_{\pm }(r,t)$ are the wave energy densities, ${U}_{{\rm{R}}}(r,t)$ is the "residual" energy density (Grappin et al. 1982, 1983), ${Q}_{\pm }(r,t)$ are the wave dissipation rates, ${F}_{\pm }(r,t)$ are the energy fluxes, and ${D}_{\pm }(r,t)$ are the contributions to the wave pressure force. The total wave energy density is given by ${U}_{\mathrm{tot}}={U}_{+}+{U}_{-}$, and the contributions from magnetic and kinetic energy are given by ${U}_{\mathrm{mag}}=({U}_{\mathrm{tot}}-{U}_{{\rm{R}}})/2$ and ${U}_{\mathrm{kin}}=({U}_{\mathrm{tot}}+{U}_{{\rm{R}}})/2$. Similarly, the total dissipation rate ${Q}_{\mathrm{tot}}={Q}_{+}+{Q}_{-}$, the total energy flux ${F}_{\mathrm{tot}}={F}_{+}+{F}_{-}$, and the total wave pressure force ${D}_{\mathrm{wp}}={D}_{+}+{D}_{-}$. Note that the nonlinear terms in the RMHD equations drop out in the energy Equations (7) (also see Appendix C of paper I). The terms ${u}_{0}{D}_{\pm }$ in the energy equations represent the work done by the wave pressure forces on the background flow.

We now derive an expression for the energy cascade rates. For a given value ${a}_{\perp }$ of the dimensionless perpendicular wavenumber, the basis functions can be split into two sets, a low-wavenumber set, $L\equiv \{k\,| \,{a}_{k}\lt {a}_{\perp }\}$, and a high-wavenumber set, $H=\{k\,| \,{a}_{k}\gt {a}_{\perp }\}$. All wave-related quantities have contributions from both low- and high-wavenumber sets. For example, the wave energy densities can be written as ${U}_{\pm }={U}_{L,\pm }+{U}_{H,\pm }$, where the subscripts L and H refer to the two subsets:

$$\begin{eqnarray}&&{U}_{L,\pm }=\displaystyle \frac{{\rho }_{0}}{4{R}^{2}}\displaystyle \sum _{k\in L}{a}_{k}^{2}{f}_{\pm ,k}^{2},\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{U}_{H,\pm }=\displaystyle \frac{{\rho }_{0}}{4{R}^{2}}\displaystyle \sum _{k\in H}{a}_{k}^{2}{f}_{\pm ,k}^{2}.\end{eqnarray} \tag{ 14 }$$

Similar expressions can be written for the other quantities listed in Equations (9) through (12). Taking the time derivatives of ${U}_{L,\pm }$ and ${U}_{H,\pm }$, and using Equation (6), we can derive separate energy equations for the low- and high-wavenumber sets:

$$\begin{eqnarray}\displaystyle \frac{\partial {U}_{L,\pm }}{\partial t}+{B}_{0}\displaystyle \frac{\partial }{\partial r}\left(\displaystyle \frac{{F}_{L,\pm }}{{B}_{0}}\right) & = & -{u}_{0}{D}_{L,\pm }\\ & & \mp \tfrac{1}{2}\displaystyle \frac{{{dv}}_{{\rm{A}}}}{{dr}}{U}_{L,{\rm{R}}}-{Q}_{L,\pm }-{\epsilon }_{\pm },\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}\displaystyle \frac{\partial {U}_{H,\pm }}{\partial t}+{B}_{0}\displaystyle \frac{\partial }{\partial r}\left(\displaystyle \frac{{F}_{H,\pm }}{{B}_{0}}\right) & = & -{u}_{0}{D}_{H,\pm }\\ & & \mp \tfrac{1}{2}\displaystyle \frac{{{dv}}_{{\rm{A}}}}{{dr}}{U}_{H,{\rm{R}}}-{Q}_{H,\pm }+{\epsilon }_{\pm }.\end{eqnarray} \tag{ 16 }$$

Here ${\epsilon }_{\pm }$ are the rates at which energy is transferred from set L to set H by nonlinear coupling:

$$\begin{eqnarray}&&{\epsilon }_{\pm }({a}_{\perp },r,t)=\displaystyle \frac{{\rho }_{0}}{4{R}^{4}}\displaystyle \sum _{k\in H}\displaystyle \sum _{j\in L}\,\displaystyle \sum _{i}\,{M}_{{kji}}({a}_{i}^{2}-{a}_{j}^{2}-{a}_{k}^{2})\\ &&{f}_{\pm ,k}{f}_{\pm ,j}{f}_{\mp ,i},\end{eqnarray} \tag{ 17 }$$

where the sum over k is restricted to set H, the sum over j can be restricted to the set L (because the contributions from $j\in H$ cancel each other), and the sum over i includes all modes. These rates are functions of dimensionless perpendicular wavenumber ${a}_{\perp }$, position r along the flux tube, and time t. The time-averaged cascade rates are given by

$$\begin{eqnarray}{\epsilon }_{\pm }({a}_{\perp },r) & = & \displaystyle \frac{{\rho }_{0}}{4{R}^{4}}\displaystyle \sum _{k\in H}\displaystyle \sum _{j\in L}\,\displaystyle \sum _{i}{M}_{{kji}}({a}_{i}^{2}-{a}_{j}^{2}-{a}_{k}^{2})\\ & & \times \langle {f}_{\pm ,k}{f}_{\pm ,j}{f}_{\mp ,i}\rangle ,\end{eqnarray} \tag{ 18 }$$

where $\langle \cdots \rangle$ denotes a time average, and the turbulence is assumed to be in a statistically stationary state. Note that the indices i, j, and k refer to three distinct modes ($i\ne j\ne k$). Also, the cascade rate ${\epsilon }_{+}$ for the dominant waves depends linearly on the amplitude ${f}_{-,i}$ of the minority waves, and conversely, the cascade rate ${\epsilon }_{-}$ for the minority waves depends linearly on ${f}_{+,i}$. Therefore, the minority waves play an important role in the cascade of the dominant waves, and vice versa (e.g., Matthaeus et al. 1999; Chandran et al. 2009).

According to Equations (17) and (18), the cascade rates ${\epsilon }_{\pm }$ at wavenumber ${a}_{\perp }$ have contributions from all mode triples $(i,j,k)$ that straddle the chosen wavenumber. In general, there are a large number of triples contributing to the overall cascade rate. Therefore, each term in the above equations has only a small contribution to the overall cascade rate. This means that the random variables ${f}_{\pm ,k}(r,t)$, ${f}_{\pm ,j}(r,t)$, and ${f}_{\mp ,i}(r,t)$ are only weakly correlated, and the triple correlation $\langle {f}_{\pm ,k}{f}_{\pm ,j}{f}_{\mp ,i}\rangle$ is small compared to the product of the typical values of the three mode amplitudes. This makes it difficult to develop a statistical theory of reflection-driven wave turbulence in an inhomogeneous atmosphere. In this paper we avoid this problem by taking the mode amplitudes from numerical simulations. We use Equation (17) to compute the time-dependent cascade rates ${\epsilon }_{\pm }({a}_{\perp },r,t)$, and we then average the results to obtain ${\epsilon }_{\pm }({a}_{\perp },r)$. This avoids any assumptions about the statistical properties of the turbulence.

In paper I we developed an RMHD model for wave turbulence in a thin flux tube inside a coronal hole. The flux tube extends from the coronal base (${r}_{\mathrm{base}}=1.003\,{R}_{\odot }$) to r = 20 ${R}_{\odot }$, well into the super-Alfvenic part of the wind. In the present version of the model, the flux tube has a square cross section of size $2R(r)$, but we still refer to R(r) as the tube radius. The method for constructing the background atmosphere is described in Appendix D of paper I. Briefly, the field strength ${B}_{0}(r)$ and temperature ${T}_{0}(r)$ are specified analytically, and the outflow velocity ${u}_{0}(r)$ is computed by solving the wind equation, including the effects of wave pressure forces on the background medium. Three models for the background atmosphere will be considered: in Models A and B the plasma density ${\rho }_{0}(r)$ and Alfvén speed ${v}_{{\rm{A}}}(r)$ vary smoothly with position along the flux tube, and in Model C there are additional density variations that simulate the effect of compressive MHD waves. The background model includes the effects of the wave pressure force ${D}_{\mathrm{wp}}$ on the outflowing plasma. For Models A and C, the rms velocity of the waves at the coronal base is assumed to be ${v}_{\mathrm{rms},\odot }=40\,\mathrm{km}\,{{\rm{s}}}^{-1}$, and for Model B we use ${v}_{\mathrm{rms},\odot }=30\,\mathrm{km}\,{{\rm{s}}}^{-1}$.

Figure 1 shows the background atmosphere for Model A, which is nearly identical to the model described in Section 3 of paper I. Figure 1(a) shows the flux tube radius R(r), which increases from 6 Mm at the coronal base to 360 Mm at r = 20 ${R}_{\odot }$. Figure 1(b) shows the assumed temperature ${T}_{0}(r)$, which is given by Equation (59) of paper I with parameter values ${C}_{0}=0.35$, ${C}_{1}=2$, m = 0.3, and k = 8. Figure 1(c) shows the plasma heating rate ${Q}_{{\rm{A}}}(r)$ necessary to maintain this temperature. The heating is assumed to be balanced by cooling due to the expansion of the outflowing plasma (${Q}_{\mathrm{adv}}$), radiative losses (${Q}_{\mathrm{rad}}$), and conductive losses (${Q}_{\mathrm{cond}}$); these contributions are shown by the colored curves in Figure 1(c). The dashed red curves indicate regions where ${Q}_{\mathrm{cond}}\lt 0$, that is, the con-vergence of the conductive flux is heating the plasma. Figure 1(d) shows the outflow velocity ${u}_{0}(r)$ and the Alfvén speed ${v}_{{\rm{A}}}(r)$. Note that the Alfvén critical point is located at $r\approx 7.2\,{R}_{\odot }$.

Zoom In Zoom Out Reset image size

On the real Sun, Alfvén or kink waves may be produced by interactions of photospheric magnetic elements with granule-scale convective flows (e.g., Spruit 1982; Edwin & Roberts 1983; Morton et al. 2013). Due to the density stratification of the lower atmosphere, the waves are significantly amplified on their way to the corona. However, the present model does not include the lower atmosphere, but starts at the coronal base. The waves are launched by imposing random "footpoint" motions on the plasma and magnetic field at the coronal base. The footpoint velocity is given by ${\boldsymbol{v}}={{\rm{\nabla }}}_{\perp }f\times {\hat{{\boldsymbol{B}}}}_{0}$, where $f(x,y,{r}_{\mathrm{base}},t)$ is the velocity stream function at the coronal base. The latter is written as a sum over basis functions:

$$\begin{eqnarray}&&f(x,y,{r}_{\mathrm{base}},t)=\displaystyle \sum _{k\in D}\,{f}_{k}({r}_{\mathrm{base}},t){\tilde{F}}_{k}(\tilde{x},\tilde{y}),\end{eqnarray} \tag{ 19 }$$

where D is a set of "driver" modes with dimensionless wavenumbers in the range $3.5\pi \lt {a}_{k}\lt 5.5\pi$. The amplitudes ${f}_{k}({r}_{\mathrm{base}},t)$ of the driver modes vary randomly with time t in the simulation. For each mode, we first create a normally distributed random sequence f(t) on a grid of times covering the entire simulation (${t}_{\max }\,=$ 30,000 s). Then the sequence is Fourier filtered using a Gaussian function $G{(\tilde{\nu })=\exp [-({\tau }_{0}\tilde{\nu })}^{2}]$, where $\tilde{\nu }$ is the temporal frequency (in Hz) and ${\tau }_{0}$ is a specified parameter. In the present work we use ${\tau }_{0}=120\,{\rm{s}}$, which corresponds to a correlation time ${\tau }_{{\rm{c}}}={\tau }_{0}/\sqrt{2\pi }\approx 48\,{\rm{s}};$ this value was chosen to be comparable to the timescale of the solar granulation. The filtered sequences are renormalized such that each driver mode has an equal contribution to the square of the velocity:

$$\begin{eqnarray}&&\displaystyle \frac{{a}_{k}^{2}}{{R}_{\mathrm{base}}^{2}}\langle {f}_{k}^{2}\rangle =\displaystyle \frac{{v}_{\mathrm{rms},\odot }^{2}}{{N}_{D}},\end{eqnarray} \tag{ 20 }$$

where $\langle \cdots \rangle$ denotes a statistical average, and N_D is the number of driver modes (N_D = 60). The driver modes are assumed to be uncorrelated, $\langle {f}_{l}{f}_{k}\rangle =0$ for $l\ne k$. We assume ${v}_{\mathrm{rms},\odot }=40\,\mathrm{km}\,{{\rm{s}}}^{-1}$, consistent with the value used in the setup of the background model. This value is also consistent with the observed spectral line widths and nonthermal velocities in coronal holes (Wilhelm et al. 1998; McIntosh et al. 2008; Banerjee et al. 2009; Landi & Cranmer 2009; Singh et al. 2011; Bemporad & Abbo 2012; Hahn et al. 2012). Note that the dynamical time of the footpoint motions is comparable to the correlation time, ${\tau }_{\mathrm{dyn}}=2{\lambda }_{\perp ,\odot }/{v}_{\mathrm{rms},\odot }=50\,{\rm{s}}$. The normalized autocorrelation function for the x component of velocity is given by

$$\begin{eqnarray}{C}_{x}({\rm{\Delta }}x,{\rm{\Delta }}y) & \equiv & \displaystyle \frac{\langle {v}_{x}(x+{\rm{\Delta }}x,y+{\rm{\Delta }}y,t){v}_{x}(x,y,t)\rangle }{\langle {v}_{x}^{2}\rangle }\\ & = & \displaystyle \frac{2}{{N}_{D}}\displaystyle \sum _{k\in D}\tfrac{{a}_{y,k}^{2}}{{a}_{k}^{2}}\cos ({a}_{x,k}{\rm{\Delta }}\tilde{x})\cos ({a}_{y,k}{\rm{\Delta }}\tilde{y}),\end{eqnarray} \tag{ 21 }$$

where $({\rm{\Delta }}x,{\rm{\Delta }}y)$ are spatial offsets, and $({\rm{\Delta }}\tilde{x},{\rm{\Delta }}\tilde{y})$ are dimensionless values of these offsets (normalized by ${R}_{\mathrm{base}}$). Figure 2(a) shows this correlation function as a color-scale plot. The anisotropy of the distribution is due to the fact that we consider here only the x component of the velocity; the correlation function for ${v}_{y}(x,y,t)$ would be rotated by 90 degrees. Figure 2(b) shows two cross sections of the correlation function, ${C}_{x}({\rm{\Delta }}x,0)$ (red curve) and ${C}_{x}(0,{\rm{\Delta }}y)$ (green curve). The circle in Figure 2(a) indicates the region where the strongest (positive or negative) correlations are found. We use the radius of this circle as our definition of the autocorrelation length ${\lambda }_{\perp ,\odot }$ of the footpoint motions. Note that ${\lambda }_{\perp ,\odot }=1$ Mm, significantly smaller than the domain half-width, ${R}_{\mathrm{base}}=6$ Mm.

Zoom In Zoom Out Reset image size

The present model neglects all details of the collisionless processes by which the waves are dissipated at small spatial scales. The "dissipation range" of the turbulence is defined as the region in wavenumber space where the simulated waves are dissipated. In the present model, the dissipation range is given by ${a}_{k}\gt (2/3){a}_{\max }$, where ${a}_{\max }$ is the maximum dimensionless wavenumber in the model. In the region below the dissipation range we set ${\nu }_{\pm ,k}=0$ so that these waves can propagate over a long distance without significant dissipation. Inside the dissipation range we use ${\nu }_{\pm ,k}={\nu }_{\pm }$, the same for all modes. The damping rate ${\nu }_{\pm }$ is twice the nonlinear cascade rate for those waves:

$$\begin{eqnarray}&&{\nu }_{\pm }=2{k}_{{\rm{d}}}{Z}_{\mp ,{\rm{d}}},\end{eqnarray} \tag{ 22 }$$

where ${k}_{{\rm{d}}}=(2/3){a}_{\max }/R$ is the wavenumber at the start of the dissipation range, and ${Z}_{\mp ,{\rm{d}}}$ is the Elsasser variable just below this range. The latter is given by

$$\begin{eqnarray}&&{Z}_{\mp ,{\rm{d}}}\equiv \sqrt{\displaystyle \sum _{i}{({a}_{i}/R)}^{2}{f}_{\mp ,i}^{2}},\end{eqnarray} \tag{ 23 }$$

where the sum is taken over all modes with wavenumbers in the range $(1/2){a}_{\max }\lt {a}_{i}\lt (2/3){a}_{\max }$. Note that the damping rate ${\nu }_{+}$ for the dominant waves depends on the Elsasser variable ${Z}_{-,{\rm{d}}}$ of the minority waves, and vice versa, so the damping rates satisfy ${\nu }_{+}\ll {\nu }_{-}$. This differs from the approach used in paper I, where we assumed ${\nu }_{+}={\nu }_{-}$. The present method has the advantage that the waves are dissipated at a rate comparable to the rate at which energy is injected into the dissipation range by the cascade. We find that this approach produces wave energy spectra that are only weakly affected by the "bottleneck" effect (e.g., Beresnyak & Lazarian 2009b).

The numerical methods for solving the RMHD Equation (6) are mostly described in Appendix B of paper I, but there are some important differences. In paper I we assumed a circular cross section, and the nonlinear terms were evaluated by summing over the nonzero elements of M_kji. As already mentioned, we now use a square domain and periodic boundary conditions. We use a set of modes with ${n}_{x,k}$ and ${n}_{y,k}$ in the range 0–21, so the maximum dimensionless wavenumber ${a}_{\max }=21\pi =65.97$, which is higher than the value ${a}_{\max }=30$ used in paper I. In the present case, the total number of modes ${k}_{\max }=1848$, and there are about 7 million matrix elements with ${M}_{{kji}}\ne 0$. Instead of summing over mode triples, we evaluate the nonlinear terms using the fast Fourier transform (FFT) method (also see Perez & Chandran 2013). In essence, the arrays ${f}_{\pm ,k}$ are transformed into functions in real space $(\tilde{x},\tilde{y})$, and the brackets are then computed as described in Equation (2). We verified that the FFT method for computing the brackets produces exactly the same result as explicitly summing over mode triples. The half-width ${R}_{\mathrm{base}}$ of the computational domain at the coronal base is 6 Mm, significantly larger than the correlation length of the footpoint motions, ${\lambda }_{\perp ,\odot }=1$ Mm. The RMHD equations are still integrated with a time step ${\rm{\Delta }}{t}_{0}=1\,{\rm{s}}$. The code is parallelized using OPENMP.

In this section we first describe RMHD simulations for Model A, which has a smooth background atmosphere (see Figure 1). The Alfvén waves launched at the coronal base produce reflection-driven turbulence at larger heights. The outward-propagating waves first reach the outer boundary of the model (r = 20 ${R}_{\odot }$) after about 10,859 s, and we simulate the turbulence for a period of 30,000 s. Figure 3 shows wave velocity patterns in cross sections of the flux tube at the end of the simulation. The first and second rows show the velocity stream functions ${f}_{\pm }(x,y)$, and the third and fourth rows show the vorticities ${\omega }_{\pm }(x,y)$. The different columns correspond to different positions along the tube and are labeled with the radial distance $r/{R}_{\odot }$. Each panel is normalized, so Figure 3 does not provide any quantitative information on the amplitude of the waves.

Zoom In Zoom Out Reset image size

Velli et al. (1989) predicted that the minority waves have both an inward-propagating "classical" component and an outward-propagating "anomalous" component. Figure 4 shows the vorticities ${\omega }_{\pm ,k}(r,t)$ of the simulated waves plotted as a function of radial distance r and time t for three different wave modes k. The selected modes have basis functions of the form ${\tilde{F}}_{k}(\tilde{x},\tilde{y})=2\cos (\pi {n}_{x,k}\tilde{x})\cos (\pi {n}_{y,k}\tilde{y})$, and the values of ${n}_{x,k}$ and ${n}_{y,k}$ are given at the top of each column of Figure 4. The upper panels show the dominant waves ${\omega }_{+,k}$, and the lower panels show the corresponding minority waves ${\omega }_{-,k}$. Note that the velocity patterns in the lower panels have the same positive slopes as the patterns in the upper panels. Therefore, the minority waves travel radially outward with the same velocity (${u}_{0}+{v}_{{\rm{A}}}$) as the dominant waves. There is no evidence in these diagrams for inward-propagating waves with negative slopes, and the same is true for all wave modes in our simulation. This means that the minority waves are dominated by the "anomalous" component (also see Perez & Chandran 2013, and paper I). Therefore, it is not correct to think of the minority waves as inward-propagating waves.

Zoom In Zoom Out Reset image size

Figure 5 shows various wave-related quantities averaged over the cross section of the flux tube and over time. Each quantity is averaged over the time interval ${t}_{0}(r)+300\leqslant t\leqslant {\rm{30,000}}$ (in seconds), where ${t}_{0}(r)$ is the time for an outward-propagating wave to reach a certain height:

$$\begin{eqnarray}&&{t}_{0}(r)\equiv {\int }_{{r}_{\mathrm{base}}}^{r}\displaystyle \frac{{{dr}}^{\prime }}{{u}_{0}({r}^{\prime })+{v}_{{\rm{A}}}({r}^{\prime })}.\end{eqnarray} \tag{ 24 }$$

The black curve in Figure 5(a) shows the rms velocity amplitude of the waves, ${v}_{\mathrm{rms}}(r)$. The solid red and green curves in Figure 5(a) show the rms values of the Elsasser variables, ${Z}_{\pm }(r)=\sqrt{\langle | {{\boldsymbol{z}}}_{\pm }{| }^{2}\rangle }$. Note that the minority waves are much weaker than the dominant waves; at $r\gt 5\,{R}_{\odot }$ the ratio ${Z}_{-}/{Z}_{+}\approx 0.016$. The function ${Z}_{-}(r)$ has a sharp minimum at $r\approx 1.3\,{R}_{\odot }$, which is due to the fact that the Alfvén speed has a maximum near that height; see Figure 1(d). Figure 5(b) shows the rms vorticity of the waves, which is dominated by waves with high perpendicular wavenumbers and is therefore more sensitive to the spatial resolution of the model. Figure 5(c) shows the rms value of the magnetic fluctuations.

Zoom In Zoom Out Reset image size

The numerical results for the Elsasser variables can be compared with predictions from a turbulence model that uses a simple phenomenology for the cascade and dissipation of waves (Chandran & Hollweg 2009). This analytical model gives the following estimates (also see Chandran et al. 2011):

$$\begin{eqnarray}&&{Z}_{+,\mathrm{est}}=\displaystyle \frac{2{v}_{\mathrm{rms},\odot }}{1+{M}_{{\rm{A}}}}{\left(\displaystyle \frac{{\rho }_{0}}{{\rho }_{0,\mathrm{base}}}\right)}^{-1/4},\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}&&{Z}_{-,\mathrm{est}}=(1+{M}_{{\rm{A}}}){\lambda }_{\perp }\left|\displaystyle \frac{{{dv}}_{{\rm{A}}}}{{dr}}\right|,\end{eqnarray} \tag{ 26 }$$

where ${v}_{\mathrm{rms},\odot }$ is the velocity amplitude at the coronal base, ${\rho }_{0,\mathrm{base}}$ is the base density, ${M}_{{\rm{A}}}(r)\equiv {u}_{0}/{v}_{{\rm{A}}}$ is the Alfvén Mach number, and ${\lambda }_{\perp }(r)$ is the perpendicular correlation length of the turbulence. The latter is estimated by extrapolation from the coronal base: ${\lambda }_{\perp }(r)={\lambda }_{\perp ,\odot }{[{B}_{0}(r)/{B}_{\mathrm{base}}]}^{-1/2}$, where ${\lambda }_{\perp ,\odot }=1$ Mm is the autocorrelation length of the footpoint velocity at the base, and ${B}_{\mathrm{base}}=10$ G is the field strength at the base. We further impose a minimum value on the Elsasser variable for the minority waves: ${Z}_{-,\mathrm{est}}\gt 2\,\mathrm{km}\,{{\rm{s}}}^{-1}$. The quantities ${Z}_{+,\mathrm{est}}$ and ${Z}_{-,\mathrm{est}}$ are plotted in Figure 5(a) as dashed red and green curves. Note that these estimates are accurate to about a factor of 2. Therefore, the model by Chandran & Hollweg (2009) indeed provides an approximate description of the Elsasser variables in the acceleration region of the wind.

Figure 5(d) shows the total energy density ${U}_{\mathrm{tot}}$ of the simulated waves (solid black curve), together with the contributions from the kinetic energy ${U}_{\mathrm{kin}}$ (red curve) and magnetic energy ${U}_{\mathrm{mag}}$ (green curve). The dashed curve shows the wave energy density ${U}_{{\rm{A}}}$ used in the setup of the background model. We see that ${U}_{\mathrm{kin}}\approx {U}_{\mathrm{mag}}$ and ${U}_{\mathrm{tot}}\approx {U}_{{\rm{A}}}$, consistent with the assumptions made in the model setup (see paper I).

Figure 5(e) shows the total energy dissipation rate ${Q}_{\mathrm{tot}}(r)$ of the simulated turbulence (solid black curve). Unlike for the model of paper I, this rate is now dominated by the contribution from damping at high perpendicular wavenumbers, and the contribution from damping at high parallel wavenumbers is no longer significant (but still included). The dissipation rate ${Q}_{\mathrm{tot}}$ is higher than that found for the smooth model in paper I, even though the background atmospheres are nearly identical. This indicates that the results of paper I are to some degree affected by a numerical artifact, namely, a "bottleneck" effect that flattens the power spectrum (see Figure 3(a) of paper I) and reduces the cascade rate for the dominant waves. The dashed black curve in Figure 5(e) shows the plasma heating rate ${Q}_{{\rm{A}}}(r)$ used in the model setup. Note that ${Q}_{\mathrm{tot}}\lt {Q}_{{\rm{A}}}$ over a significant height range in the model. Therefore, the wave dissipation rate is still smaller than the plasma heating rate needed to sustain the background atmosphere, and the model is not in thermal equilibrium. Figure 5(f) shows the same wave dissipation and plasma heating rates per unit mass. We also compare our results with predictions from a "phenomenological" turbulence model (e.g., Zhou & Matthaeus 1990; Hossain et al. 1995; Matthaeus et al. 1999; Dmitruk et al. 2001), which predicts that the dissipation rate ${Q}_{\mathrm{phen}}$ is given by

$$\begin{eqnarray}&&{Q}_{\mathrm{phen}}={c}_{{\rm{d}}}{\rho }_{0}\displaystyle \frac{{Z}_{+}^{2}{Z}_{-}+{Z}_{-}^{2}{Z}_{+}}{4{\lambda }_{\perp }}.\end{eqnarray} \tag{ 27 }$$

Here ${c}_{{\rm{d}}}$ is a dimensionless factor of order unity. The blue curve in Figure 5(f) shows the quantity ${Q}_{\mathrm{phen}}/{\rho }_{0}$ as a function of radial distance for ${c}_{{\rm{d}}}=0.1$. This value was chosen to obtain a crude fit to the actual dissipation rate ${Q}_{\mathrm{tot}}(r)$ predicted by the RMHD simulation. Without this correction factor, the above expression would significantly overestimate the dissipation rate.

Figures 6(a) and 6(b) show power spectra for the Elsasser variables as a function of dimensionless perpendicular wavenumber ${a}_{\perp }$ for four different heights in the model. For each height, we compute the wave power in individual modes with wavenumbers a_k, and then collect the results into bins in wavenumber space with ${\rm{\Delta }}{a}_{\perp }=2$ (for details see van Ballegooijen et al. 2011). These results are derived from the last 800 s of the simulation. Figure 6(a) shows the power spectra for the dominant waves. The waves are injected at ${a}_{\perp }\sim 15$, which corresponds to the correlation length ${\lambda }_{\perp ,\odot }$ of the footpoint motions. Figure 6(b) shows similar spectra for the minority waves. In the present work both spectra have approximately the same slopes, which are similar to those found in high-resolution turbulence simulations (e.g., Beresnyak & Lazarian 2008, 2009a, 2009b; Perez & Boldyrev 2009; Perez et al. 2012; Perez & Chandran 2013). We also compute temporal power spectra of dominant and minority waves, and we derive the average wave frequency ${\tilde{\omega }}_{\pm }$ as a function of dimensionless perpendicular wavenumber ${a}_{\perp }$. The results are shown in Figures 6(c) and (d) for four different heights in the model.

Zoom In Zoom Out Reset image size

Figure 7 shows the time-averaged energy cascade rates per unit mass, ${\epsilon }_{\pm }/{\rho }_{0}$. These rates are functions of dimensionless perpendicular wavenumber ${a}_{\perp }$ and radial distance r. Figures 7(a) and (b) show color-scale plots for the dominant and minority waves, respectively. Note that each plot has its own color bar, and that the cascade rates for the dominant waves are much larger than those for the minority waves. Red and blue colors indicate direct (${\epsilon }_{\pm }\gt 0$) and inverse (${\epsilon }_{\pm }\lt 0$) cascades, respectively. Figure 7(a) shows that at low heights the dominant waves have a direct cascade, but ${\epsilon }_{+}$ changes sign at r = 1.4 ${R}_{\odot }$ for wavenumbers in the range $15\lt {a}_{\perp }\lt 44$. The lower boundary of this range (${a}_{\perp }=15$) is approximately where energy is injected into the turbulence by the driver waves, and the upper boundary (${a}_{\perp }=44$) is where the waves start to be dissipated. The inverse cascade continues up to r = 2.5 ${R}_{\odot }$, but in a narrowing wavenumber range. These negative cascade rates reduce the amount of energy that can cascade into the dissipation range (${a}_{\perp }\gt 44$) and therefore affect the overall wave dissipation rate between 1.4 and 2.5 ${R}_{\odot }$. There is also a further extension of the region of the inverse cascade to larger heights ($r\gt 2.5\,{R}_{\odot }$) and low perpendicular wavenumbers (${a}_{\perp }\lt 15$), but this feature is relatively weak and does not seem to have a strong effect on the dissipation rates. In contrast, Figure 7(b) shows that the minority waves have a direct cascade at all heights.

Zoom In Zoom Out Reset image size

Figures 7(c) and (d) show plots of the cascade rates, ${\epsilon }_{\pm }/{\rho }_{0}$, as a function of perpendicular wavenumber for four different heights. These heights were chosen to represent the region of direct cascade near the coronal base (r = 1.11 ${R}_{\odot }$, red curves), the region of inverse cascade at intermediate heights (r = 2.04 ${R}_{\odot }$, green curves), and the direct cascades at large heights (r = 3.94 ${R}_{\odot }$, blue curves; r = 8.0 ${R}_{\odot }$, magenta curves). The colored squares give the wave dissipation rates ${Q}_{\pm }/{\rho }_{0}$. The dissipation rates are approximately equal to the cascade rates at the start of the dissipation range (${a}_{\perp }=44$), where the squares are plotted. Hence, the energy that cascades into the dissipation range is indeed dissipated shortly afterward, as expected. Note that at most heights the cascade rates vary strongly with wavenumber ${a}_{\perp }$, even in the "inertial" range ($15\lt {a}_{\perp }\lt 44$) where no injection or dissipation of the energy occurs. This indicates that the transport of wave energy is a highly nonlocal process: as the waves cascade, they also propagate upward in height over a significant distance. This is because the nonlinear cascade time for the dominant waves is comparable to the wave propagation time (see paper I).

What is the cause of the inverse cascade for the dominant waves? To answer this question, we first consider the production of minority waves by "reflection" of dominant waves. Converting vorticities to mode amplitudes ${f}_{\pm ,k}$ and using the fact that the minority waves are weak ($| {f}_{-,k}| \ll | {f}_{+,k}|$), Equation (6) can be written as

$$\begin{eqnarray}&&\displaystyle \frac{\partial {f}_{-,k}}{\partial t}+({u}_{0}-{v}_{{\rm{A}}})\displaystyle \frac{\partial {f}_{-,k}}{\partial r}\approx \displaystyle \frac{1}{2}(1+{M}_{{\rm{A}}})\displaystyle \frac{{{dv}}_{{\rm{A}}}}{{dr}}{f}_{+,k}\\ &&\quad +\displaystyle \frac{1}{2{R}^{2}{a}_{k}^{2}}\displaystyle \sum _{j}\,\displaystyle \sum _{i}{M}_{{kji}}({a}_{i}^{2}-{a}_{j}^{2}-{a}_{k}^{2}){f}_{-,j}{f}_{+,i},\end{eqnarray} \tag{ 28 }$$

where we omit the wave damping terms. Note that the "nonlinear" term in Equation (28) is in fact linear in the amplitudes ${f}_{-,j}(r,t)$ of the minority waves. Also, the production of minority waves is proportional to the Alfvén speed gradient ${{dv}}_{{\rm{A}}}/{dr}$, which is positive at heights below the peak in Alfvén speed ($r\lt 1.4\,{R}_{\odot }$) and negative above the peak ($r\gt 1.4\,{R}_{\odot }$). The dominant waves have a long cascade time and evolve only gradually with height. However, the minority waves have a short cascade time (about 10 s at the outer scale of the turbulence; see paper I) and respond much more rapidly to the changes in ${{dv}}_{{\rm{A}}}/{dr}$ with height. Therefore, as the waves propagate outward through the region around the peak in Alfvén speed, the dominant waves ${f}_{+,k}(r,t)$ remain more or less unchanged, while the minority waves ${f}_{-,k}(r,t)$ change sign. This reversal of the minority waves at $r\approx 1.4\,{R}_{\odot }$ occurs for most modes and can be seen in diagrams such as Figure 4. The reversal occurs even when the dominant and minority waves are not well correlated with each other, as is the case at higher wavenumbers. We now consider two heights, one just below the peak in Alfvén speed ($r={r}_{1}$) and another just above it ($r={r}_{2}$), such that the magnitudes of the gradients are the same at the two heights: ${({{dv}}_{{\rm{A}}}/{dr})}_{2}=-{({{dv}}_{{\rm{A}}}/{dr})}_{1}$. The dominant and minority waves at these heights are approximately related by

$$\begin{eqnarray}&&{f}_{+,k}({r}_{2},t)\approx {f}_{+,k}({r}_{1},t-{\rm{\Delta }}{t}_{12}),\end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray}&&{f}_{-,k}({r}_{2},t)\approx -{c}_{0}\,{f}_{-,k}({r}_{1},t-{\rm{\Delta }}{t}_{12}),\end{eqnarray} \tag{ 30 }$$

where ${c}_{0}\approx 1$ (independent of mode index k), and ${\rm{\Delta }}{t}_{12}$ is the wave propagation time between the two heights. These relationships follow from a symmetry of Equation (28): the equation remains valid when the signs of ${{dv}}_{{\rm{A}}}/{dr}$ and all ${f}_{-,k}$ are reversed but ${f}_{+,k}$ is unchanged. Inserting expressions (29) and (30) into Equation (18), we find the following relationships between the cascade rates at the two heights:

$$\begin{eqnarray}&&{\epsilon }_{+}({a}_{\perp },{r}_{2})\approx -{c}_{0}\,{\epsilon }_{+}({a}_{\perp },{r}_{1}),\end{eqnarray} \tag{ 31 }$$

$$\begin{eqnarray}&&{\epsilon }_{-}({a}_{\perp },{r}_{2})\approx {c}_{0}^{2}\,{\epsilon }_{-}({a}_{\perp },{r}_{1}).\end{eqnarray} \tag{ 32 }$$

Figure 7 indicates that at low heights there is a direct cascade for both wave types, ${\epsilon }_{\pm }({a}_{\perp },{r}_{1})\gt 0$. Then Equations (31) and (32) predict that at $r={r}_{2}$ there is an inverse cascade for the dominant waves and a direct cascade for the minority waves, ${\epsilon }_{+}({a}_{\perp },{r}_{2})\lt 0$ and ${\epsilon }_{-}({a}_{\perp },{r}_{2})\gt 0$. The inverse cascade of the dominant waves occurs at all wavenumbers ${a}_{\perp }$ because the reversal in sign of ${f}_{-,k}$ is rapidly transmitted to larger wavenumbers by the direct cascade of the minority waves.

As the dominant waves propagate farther out, they gradually adjust to the condition ${{dv}}_{{\rm{A}}}/{dr}\lt 0$, and ${\epsilon }_{+}$ becomes positive again. The adjustment of the dominant waves occurs first at large wavenumbers where the cascade times for the dominant waves are shortest, and later also at smaller wavenumbers. Therefore, in Figure 7(a) the upper boundary of the region with the inverse cascade lies at an angle in the $({a}_{\perp },r)$ plane. We conclude that the inverse cascade in Figure 7(a) is linked to the change of sign of ${{dv}}_{{\rm{A}}}/{dr}$ at r = 1.4 ${R}_{\odot }$, together with the fact that the cascade time for the dominant waves is relatively large and comparable to the wave travel time ${t}_{0}(r)$.

The background atmosphere for Model A was chosen to be the same as that used in paper I so we could directly compare our results and understand why the model of paper I gives such low wave dissipation rates. However, the outflow speed in this model reaches 800 $\mathrm{km}\,{{\rm{s}}}^{-1}$ by 20 ${R}_{\odot }$ (see Figure 1(d)), which is high considering that further acceleration may occur at larger radii. Also, the Alfvén critical point is located at $7.2\,{R}_{\odot }$, which is low compared to other models that rely on Alfvén waves to heat and accelerate the fast solar wind (e.g., Cranmer et al. 2007; Verdini et al. 2010; Chandran et al. 2011). Therefore, we now consider an alternative, Model B, which has a different set of model parameters. Three of the four parameters describing the background temperature were modified (${C}_{0}=0.30$, m = 0.35, k = 12; see Appendix D in paper I), which leads to a reduction of the peak temperature from 1.31 to 1.05 MK. The revised temperature profile is shown in Figure 8(a). We also increased the coronal base pressure from 0.1 to 0.2 $\mathrm{dyne}\,{\mathrm{cm}}^{-2}$ and decreased the wave amplitude to ${v}_{\mathrm{rms},\odot }=30\,\mathrm{km}\,{{\rm{s}}}^{-1}$, which reduces the wave pressure acceleration. Figure 8(b) shows the outflow velocity ${u}_{0}(r)$ and Alfvén speed ${v}_{{\rm{A}}}(r)$ resulting from these changes. In Model B the peak in Alfvén speed occurs at $r\approx 1.6\,{R}_{\odot }$, and the Alfvén critical point is located at r = 9.6 ${R}_{\odot }$, more in line with the values used in earlier models. The plasma heating rate ${Q}_{{\rm{A}}}(r)$ needed to maintain the background temperature is shown by the black curve in Figure 8(c), together with the cooling rates due to radiation (blue curve), thermal conduction (red curve), and solar wind expansion (green curve). A comparison with Figure 1(c) for Model A shows that the heating rate is significantly reduced. In fact, at large radii, ${Q}_{{\rm{A}}}$ becomes slightly negative, which is due to the near balance of conduction heating and expansion cooling at r = 20 ${R}_{\odot }$ in Model B.

Zoom In Zoom Out Reset image size

We simulated the dynamics of the Alfvén waves in Model B for a period of 30,000 s. The imposed footpoint motions are the same as in Model A. The solid black curve in Figure 8(d) shows the time-averaged wave dissipation rate ${Q}_{\mathrm{tot}}(r)$, and the blue curve shows the rate ${Q}_{\mathrm{phen}}(r)$ predicted by the phenomenological turbulence model, Equation (27) with ${c}_{{\rm{d}}}=0.1$. Note that the wave dissipation rate ${Q}_{\mathrm{tot}}$ is reduced compared to Model A and is below the required heating rate ${Q}_{{\rm{A}}}$ (dashed curve) for $1.5\lt r\lt 5\,{R}_{\odot }$. This is again due to the presence of an inverse cascade, in this case in the height range $1.6\lt r\lt 4\,{R}_{\odot }$. The heating produced by the waves is insufficient to maintain the assumed background temperature, so the model is not in thermal equilibrium.

In this section we describe simulation results for Model C, which has a background atmosphere with additional density variations along the flux tube. These variations simulate the effect that compressive MHD waves may have on the propagation and reflection of Alfvén waves. The density variations $\delta {\rho }_{0}(r)$ are assumed to be random in position but constant in time, consistent with our RMHD methodology. The present model is similar to that described in Section 4 of paper I, but the model parameters are different. The magnitude of the density variations has been doubled (to ${\epsilon }_{\mathrm{rms}}=0.2$), and the correlation length of the variations has been increased by a factor of 5 (to ${\lambda }_{{\rm{c}}}=0.2\,{R}_{\odot }$). Therefore, the magnitude of the variations in the Alfvén speed gradient ${{dv}}_{{\rm{A}}}/{dr}$ has decreased by a factor of 2.5, reducing the wave reflection. The main reason for these changes is to bring out more clearly the spatial variations of the cascade rates ${\epsilon }_{\pm }$ of the dominant and minority waves. In reality, there are density fluctuations both along and perpendicular to the field lines, and the latter may actually be much larger than the former (see references in Section 4 of paper I). We believe that perpendicular density fluctuations will have an important effect on the cascade rates, but unfortunately we are unable to simulate this effect with our RMHD code, which assumes constant density over the cross section of the flux tube. To compensate, we use a high value for the magnitude of the density fluctuations along the field lines. This approach is rather artificial, but it is the best we can do right now.

In model C, the Alfvén speed gradient ${{dv}}_{{\rm{A}}}/{dr}$ changes sign multiple times with increasing r, which significantly enhances the wave reflection. Figure 9 shows the vorticities ${\omega }_{\pm ,k}(r,t)$ for three different wave modes, the same modes as in Figure 4. The upper panels of Figure 9 show the dominant waves ${\omega }_{+,k}$, and the lower panels show the corresponding minority waves ${\omega }_{-,k}$. The positive slopes in the upper panels indicate that the dominant waves travel radially outward with velocity ${u}_{0}+{v}_{{\rm{A}}}$. However, the velocity patterns in Figure 9(d) have a negative slope, indicating this low-wavenumber minority wave (a_k = 4.44) is propagating radially inward with velocity ${u}_{0}-{v}_{{\rm{A}}}\lt 0;$ this is the "classical" component of the minority waves described by Velli et al. (1989). In Figures 9(e) and (f) the minority waves have mostly outward-propagating components. All panels show a stationary (vertical) pattern, which is an artifact of our assumption that the density variations $\delta {\rho }_{0}(r)$ are constant in time. At high wavenumbers the minority waves are uncorrelated with the dominant waves (compare Figures 9(f) and (c)). These patterns are quite different from those for the smooth Model A (see Figure 4).

Zoom In Zoom Out Reset image size

Figure 10 shows the time-averaged cascade rates ${\epsilon }_{\pm }({a}_{\perp },r)$ for the model with density variations. Figure 10(a) shows the cascade rate ${\epsilon }_{+}/{\rho }_{0}$ for the dominant waves. Note that with increasing height r the cascade rate changes sign multiple times, going from a direct cascade (red) to an inverse cascade (blue) over short distances. In contrast, the minority waves always have a direct cascade; see Figure 10(b). The changes in ${\epsilon }_{+}$ as a function of r are due to changes in ${{dv}}_{{\rm{A}}}/{dr}$, but the dominant waves try to adjust to these changes, so at heights above $3\,{R}_{\odot }$ the quantities ${\epsilon }_{+}$ and ${{dv}}_{{\rm{A}}}/{dr}$ are only poorly correlated. A comparison of the color bars in Figures 7 and 10 indicates that the magnitudes of the cascade rates in Model C are much larger than those in Model A. Figures 10(c) and (d) show the cascade rates ${\epsilon }_{\pm }/{\rho }_{0}$ as functions of wavenumber for four different heights. The colored squares indicate the wave dissipation rates ${Q}_{\pm }/{\rho }_{0}$. Note that for the dominant waves the cascade rates vary strongly with wavenumber, and the dissipation rates ${Q}_{+}$ are much smaller than the peak cascade rates $| {\epsilon }_{+}{| }_{\max }$, which generally occur at low wavenumber. This indicates that the rapid changes in ${\epsilon }_{+}$ as a function of r prevent the efficient cascade of wave energy to higher wavenumbers and thereby have a negative effect on the dissipation rate. The total dissipation rate ${Q}_{\mathrm{tot}}(r)$ for Model C is slightly larger than that for Model A, but it is still insufficient to maintain the temperature of the background atmosphere. Therefore, Model C is also not in thermal equilibrium.

Zoom In Zoom Out Reset image size

In this work we simulate the dynamics of Alfvén waves for three models of the fast solar wind, two with a smooth background atmosphere (Models A and B) and one with density fluctuations (Model C). These models are improved versions of the models presented in paper I. We compute for the first time the energy cascade rates ${\epsilon }_{\pm }({a}_{\perp },r)$ for dominant and minority waves, and we find that at certain heights and wavenumbers the dominant waves undergo an inverse cascade, ${\epsilon }_{+}\lt 0$. This means that the nonlinear interactions between dominant and minority waves cause energy to be transported from smaller to larger scales, opposite to the direction usually assumed for Alfvén wave turbulence. Inverse cascades are predicted to occur in two-dimensional magnetohydrodynamic systems (Fyfe & Montgomery 1976; Fyfe et al. 1977) and in three-dimensional systems with large magnetic helicity (Frisch et al. 1975; Dmitruk & Matthaeus 2007). In the present case, the inverse cascade appears to be due to a different mechanism, namely, the change in sign of the gradient of the Alfvén speed ${{dv}}_{{\rm{A}}}/{dr}$ as a function of height in the model. In Model A the inverse cascade occurs for radii between 1.4 and 2.5 ${R}_{\odot }$ and wavenumbers in the inertial range ($15\lt {a}_{\perp }\lt 44$); in Model B the inverse cascade occurs between 1.6 and $4\,{R}_{\odot }$. In Model C with density fluctuations, the cascade rate changes sign multiple times as a function of height. In all models the cascade rate ${\epsilon }_{+}({a}_{\perp },r)$ varies significantly with perpendicular wavenumber. This indicates that wave propagation plays an important role in the cascade process. The cascade timescale for the dominant waves is comparable to the wave travel time ${t}_{0}(r)$, so in one cascade time the waves travel a distance comparable to the radial distance r. The energy injected into the cascade by the driver waves at position r is not dissipated locally, but is dissipated only much later at significantly larger heights. Therefore, the cascade of the dominant waves is a highly nonlocal process.

The inverse cascade impedes the efficient transport of wave energy to large perpendicular wavenumbers, and thereby has a negative effect on the wave dissipation rate. Therefore, the low dissipation rates found here (and in paper I) are a real physical effect, not a numerical artifact. In both the smooth models and the model with density fluctuations, there is a significant height range where ${Q}_{\mathrm{tot}}(r)$ is less than the plasma heating rate ${Q}_{{\rm{A}}}(r)$ needed to maintain the temperature of the background atmosphere. Hence, these models are not self-consistent from an energy point of view. Obtaining a model with higher dissipation rates would require that the inverse cascade is somehow avoided, or at least reduced in magnitude. At present it is unclear how to construct such a model.

In the smooth models, the total dissipation rate ${Q}_{\mathrm{tot}}(r)$ is much lower than expected from a "phenomenological" turbulence model, Equation (27) with ${c}_{{\rm{d}}}=1$. The phenomenological model is based on the assumption that the cascade rates ${\epsilon }_{\pm }(r,{a}_{\perp })$ are approximately constant with wavenumber ${a}_{\perp }$ in the inertial range and that the cascade is everywhere direct. These assumptions are reasonable for a fully developed turbulent system in a closed box or periodic domain, but they are not appropriate for the extended corona where the cascade time for the dominant waves is comparable to the wave travel time and the dissipation is highly nonlocal. Indeed, we find that the cascade rate ${\epsilon }_{+}$ varies significantly with wavenumber in the inertial range. This fact prevents the straightforward application of the phenomenological formalism.

We predict that for models with a monotonically decreasing Alfvén speed ${v}_{{\rm{A}}}(r)$ a direct cascade rate occurs at all heights. To verify this prediction, we constructed a model (not shown) with a temperature ${T}_{0}(r)$ that decreases monotonically with height. The model was constructed by setting ${C}_{0}=0.3$, ${C}_{1}=-2$, and k = 4 in Equation (59) of paper I; the coronal base pressure was assumed to be 0.03 $\mathrm{dyne}\,{\mathrm{cm}}^{-2}$. In this model the Alfvén speed ${v}_{{\rm{A}}}(r)$ decreases monotonically with height. We find that indeed the cascade rate ${\epsilon }_{+}\gt 0$ at all heights for wavenumbers larger than those of the driver waves (${a}_{\perp }\gt 15$). However, such a model for ${v}_{{\rm{A}}}(r)$ is not realistic. In the chromospheric-corona transition region the temperature ${T}_{0}(r)$ must increase with height, and the density must rapidly decrease, leading to a local increase in Alfvén speed ${v}_{{\rm{A}}}(r)$ with height. In the corona we expect high Alfvén speeds, ${v}_{{\rm{A}}}\sim 1000\,\mathrm{km}\,{{\rm{s}}}^{-1}$. In the solar wind the Alfvén speed is again relatively low, ${v}_{{\rm{A}}}\sim 30\,\mathrm{km}\,{{\rm{s}}}^{-1}$ at 1 au. Therefore, for a realistic model of the fast solar wind, the Alfvén speed ${v}_{{\rm{A}}}(r)$ must have a maximum at some height in the corona. Hence, an inverse cascade of the kind found in the present models cannot easily be avoided.

The models presented here and in paper I differ significantly from other models in which the driver waves at the coronal base are assumed to have large perpendicular correlation lengths (10–30 Mm) and long correlation times (tens of minutes or longer; e.g., Verdini et al. 2009; Perez & Chandran 2013). Such large lengths and timescales are comparable to those of the supergranulation, which is a convective flow pattern observed in the solar photosphere. Long-period waves are more strongly reflected in the extended corona and produce stronger minority waves and higher wave dissipation rates than the present model. However, as argued in paper I, we do not believe that the supergranulation can play a significant role in producing the transverse waves that drive the solar wind. The reason is that supergranular flows have low velocity ($\sim 0.3\,\mathrm{km}\,{{\rm{s}}}^{-1}$), and the magnetic elements in the lower atmosphere respond quasi-statically to such weak, slowly varying flows. The buoyant magnetic flux elements in the photosphere are expected to be passively advected by these horizontal flows without much amplification of the motions with increasing height. Therefore, the supergranular flows are expected to produce velocities of only about 0.3 $\mathrm{km}\,{{\rm{s}}}^{-1}$ in the low corona, not the much larger velocities assumed by Verdini et al. (2009) and Perez & Chandran (2013).

In contrast, in the present model the waves are assumed to be produced by interactions of magnetic elements with the solar granulation, which is a different convective pattern with a typical length scale of 1 Mm and a timescale of a few minutes. The magnetic structures in the lower atmosphere respond dynamically to such short-period disturbances (e.g., van Ballegooijen et al. 2014), producing transverse waves that are amplified from about 1 km s⁻¹ in the photosphere to about 30 km s⁻¹ in the low corona. Therefore, in the present model we assume a relatively small correlation length (${\lambda }_{\perp ,\odot }=1$ Mm) and a short correlation time (${\tau }_{{\rm{c}}}=50$ s). This produces less wave reflection than in models with long-period waves, and we find that the minority waves are predominantly of the outward-propagating "anomalous" type, whereas Verdini et al. (2009) find that for long-period driving there is also a significant component of inward-propagating "classical" waves. Although we find an inverse cascade, its effects are limited to length scales smaller than those of the driver waves, and we do not find a strong tendency for the turbulence to cascade to larger scales. Distinguishing between the different turbulence models will require further observations of Alfvén waves in coronal holes, including their typical length and timescales and the relative amplitudes of outward- and inward-propagating waves. The observations by Morton et al. (2015) are an important step in that direction.

A number of authors have suggested that density fluctations and coupling between Alfvén and compressive MHD waves play an important role in the heating of the fast solar wind (e.g., Kudoh & Shibata 1999; Moriyasu et al. 2004; Suzuki & Inutsuka 2005, 2006; Matsumoto & Shibata 2010). Here we find that RMHD models with fixed density variations in the background atmosphere have inverse cascades that limit the rate at which wave energy can be dissipated. In such models, sufficient energy is available at large spatial scales, but the energy is not efficiently cascaded to small scales where it can be dissipated. To obtain a more efficient cascade process, it may be necessary to go beyond standard RMHD modeling and include the effects of perpendicular density variations, which give rise to phase mixing and resonant absorption of the waves (e.g., Heyvaerts & Priest 1983; De Groof & Goossens 2002; Goossens et al. 2011, 2012, 2013; Pascoe et al. 2012). Future modeling of Alfvén wave turbulence in the acceleration region of the fast wind should take such transverse density variations into account.

We thank the referee for useful comments that led us to explore alternative models for the background atmosphere. We thank Benjamin Chandran for comments that helped improve the presentation of the paper. We are most grateful to Alex Voss from the School of Computer Science at the University of St. Andrews for his support with the computational work. This project was supported under contract NNM07AB07C from NASA to the Smithsonian Astrophysical Observatory (SAO) and SP02H1701R from LMSAL to SAO.