3.1 Radial profiles of $\unicode[STIX]{x1D70C}$ , $B_{0}$ and $U$

We choose the radial profiles of $\unicode[STIX]{x1D70C}$ , $U$ and  $B_{0}$ to approximate the conditions found in coronal holes and the fast solar wind. Above the transition region, at $r>r_{\text{tr}}$ , we set

(3.4) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}=(10^{9}s^{-15.6}+2.51\times 10^{6}s^{-3.76}+1.85\times 10^{5}s^{-2})m_{\text{p}}~\text{cm}^{-3}, & \displaystyle\end{eqnarray}$$

(3.5) $$\begin{eqnarray}\displaystyle & \displaystyle B_{0}=1.5[s^{-6}(f_{\text{max}}-1)+s^{-2}]\,\text{G}, & \displaystyle\end{eqnarray}$$

and

(3.6) $$\begin{eqnarray}U=9.25\times 10^{12}\left(\frac{B_{0}}{1\,\text{G}}\right)\left(\frac{\unicode[STIX]{x1D70C}}{m_{\text{p}}~\text{cm}^{-3}}\right)^{-1}\text{cm}\,\text{s}^{-1},\end{eqnarray}$$

where

(3.7) $$\begin{eqnarray}\displaystyle & \displaystyle s=\frac{r}{R_{\odot }}, & \displaystyle\end{eqnarray}$$

(3.8) $$\begin{eqnarray}\displaystyle & \displaystyle f_{\text{max}}=9 & \displaystyle\end{eqnarray}$$

is the super-radial expansion factor, and $m_{\text{p}}$ is the proton mass. Equation (3.4) is adapted from the coronal-hole electron-density measurements of Feldman et al. (Reference Feldman, Habbal, Hoogeveen and Wang1997). We have modified those authors’ density profile by adding the $s^{-2}$ term in (3.4) so that the model extrapolates to a reasonable density at large  $r$ and by increasing the coefficient of the $s^{-15.6}$ term in order to match the low-corona density in the model of Cranmer & van Ballegooijen (Reference Cranmer and van Ballegooijen2005). Equation (3.5) is taken from Hollweg & Isenberg (Reference Hollweg and Isenberg2002). The general form of (3.6) follows from (2.17). The numerical coefficient on the right-hand side of (3.6) is chosen so that

(3.9a,b ) $$\begin{eqnarray}U(r_{\text{b}})=1.2~\text{km}~\text{s}^{-1}\quad U(1~\text{au})=750~\text{km}~\text{s}^{-1},\end{eqnarray}$$

where

(3.10) $$\begin{eqnarray}r_{\text{b}}=1.0027R_{\odot }\end{eqnarray}$$

is a heliocentric distance just larger than  $r_{\text{tr}}$ that we take to correspond to the base of the corona. Given the radial profiles in (3.4) through (3.6), the Alfvén critical point is at

(3.11) $$\begin{eqnarray}r_{\text{A}}=11.1R_{\odot },\end{eqnarray}$$

the Alfvén speed reaches its maximum value at

(3.12) $$\begin{eqnarray}r_{\text{m}}=1.71R_{\odot },\end{eqnarray}$$

and

(3.13a-c ) $$\begin{eqnarray}v_{\text{A}}(r_{\text{b}})=935~\text{km}~\text{s}^{-1}\quad v_{\text{A}}(r_{\text{m}})=2730~\text{km}~\text{s}^{-1}\quad v_{\text{A}}(r_{\text{A}})=U(r_{\text{A}})=627~\text{km}~\text{s}^{-1}.\end{eqnarray}$$

Below the transition region, we set

(3.14) $$\begin{eqnarray}\unicode[STIX]{x1D70C}=\unicode[STIX]{x1D70C}_{\text{ph}}\text{e}^{c(1-s)/s},\end{eqnarray}$$

where

(3.15) $$\begin{eqnarray}\unicode[STIX]{x1D70C}_{\text{ph}}=4.78\times 10^{16}m_{\text{p}}\text{cm}^{-3}\end{eqnarray}$$

is the photospheric density, $c=[s_{\text{tr}}/(1-s_{\text{tr}})]\ln (\unicode[STIX]{x1D70C}_{\text{tr},<}/\unicode[STIX]{x1D70C}_{\text{ph}})$ , $s_{\text{tr}}=r_{\text{tr}}/R_{\odot }$ and $\unicode[STIX]{x1D70C}_{\text{tr},<}$ is the density just below the transition region, which we take to be 100 times greater than the value of the density at $r=r_{\text{tr}}$ from (3.4). We then set (cf. van Ballegooijen et al. Reference van Ballegooijen, Asgari-Targhi, Cranmer and DeLuca2011)

(3.16) $$\begin{eqnarray}B=\left[\frac{(B_{\text{ph}}^{2}-B_{\text{tr}}^{2})(\unicode[STIX]{x1D70C}-\unicode[STIX]{x1D70C}_{\text{tr},<})}{\unicode[STIX]{x1D70C}_{\text{ph}}-\unicode[STIX]{x1D70C}_{\text{tr},<}}+B_{\text{tr}}^{2}\right]^{1/2},\end{eqnarray}$$

at $r<r_{\text{tr}}$ , where

(3.17) $$\begin{eqnarray}B_{\text{ph}}=1400~\text{G}\end{eqnarray}$$

is the assumed magnetic-field strength in the photospheric footpoint of the simulated flux tube, and $B_{\text{tr}}$ is the value of  $B$ at $r=r_{\text{tr}}$ from (3.5).

We plot the radial profiles of $\unicode[STIX]{x1D70C}$ , $B$ , $U$ and $v_{\text{A}}$ in figure 2. We also plot the $\boldsymbol{z}^{+}$ travel time between the photosphere and radius  $r$ ,

(3.18) $$\begin{eqnarray}T(r)=\int _{R_{\odot }}^{r}\frac{\text{d}r}{U+v_{\text{A}}}.\end{eqnarray}$$

Figure 2. The radial profiles of the solar-wind outflow velocity  $U$ , Alfvén speed  $v_{\text{A}}$ , plasma density  $\unicode[STIX]{x1D70C}$ divided by the proton mass  $m_{\text{p}}$ , background magnetic-field strength  $B_{0}$ and $\boldsymbol{z}^{+}$ travel time from the transition region  $T(r)$ in our direct numerical simulations. We use the same profiles when evaluating quantities in the analytic model that we present in § 4.

3.2 Boundary conditions

We take the $\boldsymbol{z}^{\pm }$ fluctuations to satisfy periodic boundary conditions in the $xy$ -plane. At the photosphere, we impose a time-dependent velocity field. We set the velocity Fourier components at the photosphere equal to zero when $k_{\bot }>3\times 2\unicode[STIX]{x03C0}/L_{\text{box}}(R_{\odot })$ , where

(3.19) $$\begin{eqnarray}k_{\bot }=\sqrt{k_{x}^{2}+k_{y}^{2}},\end{eqnarray}$$

and $k_{x}$ and $k_{y}$ are the $x$ and $y$ components of the wave vector  $\boldsymbol{k}$ . We set the amplitudes of the velocity Fourier components at $k_{\bot }\leqslant 3\times 2\unicode[STIX]{x03C0}/L_{\text{box}}(R_{\odot })$ equal to a constant, which we choose so that the r.m.s. amplitude of the fluctuating velocity at the photosphere is

(3.20) $$\begin{eqnarray}\unicode[STIX]{x1D6FF}v_{\text{ph},\text{rms}}=1.3~\text{km}~\text{s}^{-1},\end{eqnarray}$$

consistent with observational constraints on the velocities of solar granules (Richardson & Schwarzschild Reference Richardson and Schwarzschild1950). We then assign random values to the phases of these velocity Fourier components at the discrete set of times $t_{\text{n}}=n\,\unicode[STIX]{x1D70F}_{0}$ , where $\unicode[STIX]{x1D70F}_{0}=5~\text{min}$ in Run 1 and $\unicode[STIX]{x1D70F}_{0}=20~\text{min}$ in Runs 2 and 3. To determine the phases at times between successive $t_{\text{n}}$ , we use cubic interpolation in time. We define the correlation time of the photospheric velocity  $\unicode[STIX]{x1D70F}_{v}^{\text{(ph)}}$ to be the time lag over which the normalized velocity autocorrelation function decreases from 1 to 0.5. The resulting velocity correlation times are listed in table 1.

Our choices of $\unicode[STIX]{x1D70F}_{0}$ and $L_{\text{box}}(R_{\odot })$ determine (at least in part – see § 3.7) the correlation time  $\unicode[STIX]{x1D70F}_{\text{c}}$ and perpendicular correlation length  $L_{\bot }$ of the AWs launched by the Sun. (Since we only drive photospheric velocity modes with $k_{\bot }\leqslant 3\times 2\unicode[STIX]{x03C0}/L_{\text{box}}(R_{\odot })$ , $L_{\bot }$ is a few times smaller than  $L_{\text{box}}$ .) Estimates of $L_{\bot }(r_{\text{b}})$ range from  $\simeq 10^{3}~\text{km}$ (Cranmer et al. Reference Cranmer, van Ballegooijen and Edgar2007; Hollweg et al. Reference Hollweg, Cranmer and Chandran2010; van Ballegooijen & Asgari-Targhi Reference van Ballegooijen and Asgari-Targhi2016; van Ballegooijen & Asgari-Targhi Reference van Ballegooijen and Asgari-Targhi2017) to more than $10^{4}~\text{km}$ (Dmitruk et al. Reference Dmitruk, Matthaeus, Milano, Oughton, Zank and Mullan2002; Verdini & Velli Reference Verdini and Velli2007; Verdini et al. Reference Verdini, Grappin, Pinto and Velli2012), and estimates of $\unicode[STIX]{x1D70F}_{\text{c}}(r_{\text{b}})$ range from $\simeq 1$ –5 min (Cranmer & van Ballegooijen Reference Cranmer and van Ballegooijen2005; van Ballegooijen & Asgari-Targhi Reference van Ballegooijen and Asgari-Targhi2016; van Ballegooijen & Asgari-Targhi Reference van Ballegooijen and Asgari-Targhi2017) to one or more hours (Dmitruk & Matthaeus Reference Dmitruk and Matthaeus2003). Given the uncertainty in $L_{\bot }(r_{\text{b}})$ and $\unicode[STIX]{x1D70F}_{\text{c}}(r_{\text{b}})$ , we vary $L_{\text{box}}(R_{\odot })$ and $\unicode[STIX]{x1D70F}_{0}$ by factors of 4 and 5, respectively, in our different simulations in order to investigate how the values of $L_{\bot }(r_{\text{b}})$ and $\unicode[STIX]{x1D70F}_{\text{c}}(r_{\text{b}})$ influence the properties of the turbulence at larger  $r$ .

No information flows into the simulation domain through the outer boundary at $r=r_{\text{max}}$ , because $r_{\text{max}}>r_{\text{A}}$ . We thus do not impose an additional boundary condition at the outer boundary.

3.3 Hyper-dissipation

To dissipate the fluctuation energy that cascades to small wavelengths, we add a hyper-dissipation term of the form

(3.21) $$\begin{eqnarray}D_{g}=-\unicode[STIX]{x1D708}_{g}\left(\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}x^{2}}+\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}y^{2}}\right)^{4}\boldsymbol{g}\end{eqnarray}$$

to the right-hand side of (2.19), and a hyper-dissipation term of the form

(3.22) $$\begin{eqnarray}D_{f}=-\unicode[STIX]{x1D708}_{f}\left(\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}x^{2}}+\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}y^{2}}\right)^{4}\boldsymbol{f}\end{eqnarray}$$

to the right-hand side of (2.20). We choose the magnitude and radial dependence of the hyper-dissipation coefficients  $\unicode[STIX]{x1D708}_{g}$ and $\unicode[STIX]{x1D708}_{f}$ so that dissipation becomes important near the grid scale at all radii in each simulation. In particular, we take $\unicode[STIX]{x1D708}_{g}$ and $\unicode[STIX]{x1D708}_{f}$ to be proportional to  $[L_{\text{box}}(r)/L_{\text{box}}(R_{\odot })]^{8}$ .

3.4 Numerical algorithm

The REFLECT code solves (2.19) and (2.20) using a spectral element method based on a Chebyshev–Fourier basis (Canuto et al. Reference Canuto, Hussaini, Quarteroni and Zang1988). In each of our three simulations, we split the numerical domain into 1024 subdomains. Each subdomain covers the full flux-tube cross-section pictured in figure 1 using 256 grid points along both the $x$ and $y$ directions, but only part of the flux tube’s radial extent. Along the  $r$ axis, each subdomain contains 17 grid points, two of which are boundary grid points. The total number of radial grid points is 16 385. Except at $r_{\text{min}}$ and $r_{\max }$ , these boundary grid points are shared by neighbouring subdomains. Eight of the subdomains are in the chromosphere.

Figure 3. Panels (a,b,c) show the r.m.s. amplitudes of $\boldsymbol{z}^{\pm }$ in Runs 1 through 3 and in the analytic model described in § 4. The lower-right panel shows $\unicode[STIX]{x1D6FF}B_{\text{rms}}/B_{0}$ in Runs 1 through 3, where $\unicode[STIX]{x1D6FF}B_{\text{rms}}$ is the r.m.s. amplitude of the magnetic-field fluctuation.

A Chebyshev/Fourier transform of (2.19) and (2.20) leads to a system of ordinary differential equations for the Chebyshev–Fourier coefficients in each subdomain. These equations are coupled through matching conditions (continuity of $\unicode[STIX]{x1D6FF}\boldsymbol{v}$ and $\unicode[STIX]{x1D6FF}\boldsymbol{B}$ ) at the boundaries between neighbouring subdomains. The REFLECT code advances the solution forward in time using a third-order Runge–Kutta method, with an integrating factor to handle the hyper-dissipation terms. Within each subdomain, the REFLECT code discretizes the radial interval using a Gauss–Lobatto grid, which makes it possible to compute the Chebyshev transform using a fast cosine transform.

3.5 Duration of the simulations

We run each simulation from $t=0$ until $t=13.2~\text{h}$ . Between $t=0$ and $t=4~\text{h}$ , the magnetic and kinetic energies in the simulations fluctuate while trending upwards. For reference, it takes 1.3 h for an outward-propagating AW to travel from the photosphere to the Alfvén critical point at $r_{\text{A}}=11.1R_{\odot }$ , and 3 h for an outward-propagating AW to travel from the photosphere to  $r_{\text{max}}=21R_{\odot }$ (see figure 2). After $t\simeq 4~\text{h}$ , the magnetic and kinetic energies fluctuate around a steady value. We regard the turbulence as being in a statistical steady state at $t>6~\text{h}$ . All the numerical results that we present are calculated from time averages between $t=6~\text{h}$ and $t=13.2~\text{h}$ , except for the $z_{\text{HF},\text{rms}}^{+}$ and $z_{\text{LF},\text{rms}}^{+}$ profiles in Run 2; those profiles, because of technical difficulties, were only computed from averages between $t=12~\text{h}$ and $t=13~\text{h}$ .

3.6 Radial profiles of the fluctuation amplitudes

In figure 3, we plot the r.m.s. amplitudes of $\boldsymbol{z}^{\pm }$ , denoted $z_{\text{rms}}^{\pm }$ , as a function of  $r$ in Runs 1 through 3 and in the analytic model discussed in § 4. The lower-right panel of figure 3 shows the fractional variation in the magnetic-field strength as a function of  $r$ in our three numerical simulations. In all three simulations, $z_{\text{rms}}^{+}\simeq z_{\text{rms}}^{-}$ in the chromosphere, because of strong AW reflection at the transition region and photosphere. On the other hand, $z_{\text{rms}}^{+}\gg z_{\text{rms}}^{-}$ in the corona and solar wind because of the limited efficiency of reflection in these regions and because $\boldsymbol{z}^{-}$ fluctuations are rapidly cascaded to small scales by the large-amplitude $\boldsymbol{z}^{+}$ fluctuations.

The value of $z_{\text{rms}}^{+}$ increases between $r=R_{\odot }$ and $r=5R_{\odot }$ because of the radially decreasing density profile. Equation (2.19) implies that the r.m.s. amplitude of $\boldsymbol{g}$ ( $g_{\text{rms}}$ ) is independent of  $r$ when (i) the fluctuations are in a statistical steady state, (ii) $z_{\text{rms}}^{-}\ll z_{\text{rms}}^{+}$ and (iii) nonlinear interactions can be ignored. At $r<5R_{\odot }$ , $\unicode[STIX]{x1D70C}(r)\gg \unicode[STIX]{x1D70C}(r_{\text{A}})$ , and it follows from (2.15) that $z_{\text{rms}}^{+}\propto g_{\text{rms}}\unicode[STIX]{x1D70C}^{-1/4}$ . Equations (2.15) and (2.19) thus imply that the linear physics of AW propagation causes $z_{\text{rms}}^{+}$ to increase rapidly with increasing  $r$ at $r<5R_{\odot }$ , since $z_{\text{rms}}^{-}\ll z_{\text{rms}}^{+}$ in this region. When nonlinear interactions are taken into account, $g_{\text{rms}}$ becomes a decreasing function of  $r$ , but the linear physics ‘wins out’ at $r<5R_{\odot }$ , in the sense that $z_{\text{rms}}^{+}\propto g_{\text{rms}}\unicode[STIX]{x1D70C}^{-1/4}$ remains an increasing function of  $r$ . Since the rate of non-WKB reflection vanishes at $r=r_{\text{m}}=1.71R_{\odot }$ , the $z^{-}$ fluctuations seen at $r=r_{\text{m}}$ in all three simulations must be generated elsewhere. At $r<r_{\text{A}}$ , $z^{-}$ fluctuations propagate with a negative radial velocity once they are produced, and thus the $z^{-}$ fluctuations seen at $r=r_{\text{m}}$ in the simulations originate at $r>r_{\text{m}}$ .

3.7 Two components of outward-propagating fluctuations

In our simulations, the transition region, which acts like an AW antenna, is characterized by two time scales at the perpendicular outer scale of the turbulence, which we take to be

(3.23) $$\begin{eqnarray}L_{\bot }={\textstyle \frac{1}{3}}L_{\text{box}}.\end{eqnarray}$$

The first time scale is the correlation time of the photospheric velocity field, $\unicode[STIX]{x1D70F}_{v}^{\text{(ph)}}$ , which we define as the time increment required for the normalized velocity autocorrelation function at the photosphere to decrease from 1 to 0.5. This time increment is 3.3 min, 9.6 min and 9.3 min in Runs 1, 2 and 3, respectively, as displayed in table 1. The second time scale is the nonlinear time scale

(3.24) $$\begin{eqnarray}\unicode[STIX]{x1D70F}_{\text{nl}}=\frac{L_{\bot }}{z_{\text{rms}}^{\pm }}\end{eqnarray}$$

of the balanced turbulence (‘balanced’ meaning that $z_{\text{rms}}^{+}\simeq z_{\text{rms}}^{-}$ ) just below the transition region at $r=r_{\text{tr},<}=r_{\text{tr}}-\unicode[STIX]{x1D716}$ , where $\unicode[STIX]{x1D716}$ is an infinitesimal distance, and $z_{\text{rms}}^{\pm }(r_{\text{tr},<})\simeq 30~\text{km}~\text{s}^{-1}$ . (Section 3.10 discusses an effect that shortens this second time scale relative to the estimate in (3.24) in Runs 1 and 2.) Although the right-hand side of (3.24) contains a $\pm$ sign, we do not include a $\pm$ sign on the left-hand side, because we will only evaluate (3.24) at locations at which $z_{\text{rms}}^{+}\simeq z_{\text{rms}}^{-}$ . We define

(3.25) $$\begin{eqnarray}\unicode[STIX]{x1D70F}_{\text{nl}}^{\text{(tr)}}=\unicode[STIX]{x1D70F}_{\text{nl}}(r_{\text{tr},<}).\end{eqnarray}$$

Given the values of $L_{\text{box}}(r_{\text{tr}})$ listed in table 1, $\unicode[STIX]{x1D70F}_{\text{nl}}^{\text{(tr)}}$ is 0.8 min, 0.8 min and 3 min in Runs 1, 2 and 3, respectively, values that are several times smaller than  $\unicode[STIX]{x1D70F}_{v}^{\text{(ph)}}$ . This suggests that the transition region in our simulations launches two populations of $\boldsymbol{z}^{+}$ fluctuations characterized by different time scales and hence different radial correlation lengths.

Figure 4. Root mean square amplitudes of $\boldsymbol{z}_{\text{HF}}^{+}$ and $\boldsymbol{z}_{\text{LF}}^{+}$ (defined in (3.26) through (3.28) and (3.32)) in Runs 1 through 3 and in the analytic model described in § 4.

To investigate this possibility, we define

(3.26) $$\begin{eqnarray}\boldsymbol{g}_{\text{LF}}(\tilde{x},{\tilde{y}},r,t)=\frac{1}{2\unicode[STIX]{x1D6E5}}\int _{r_{\text{i}}}^{r_{\text{i}}+2\unicode[STIX]{x1D6E5}}\text{d}r^{\prime }\boldsymbol{g}(\tilde{x},{\tilde{y}},r^{\prime },t),\end{eqnarray}$$

(3.27) $$\begin{eqnarray}g_{\text{LF},\text{rms}}=\langle |\boldsymbol{g}_{\text{LF}}|^{2}\rangle ^{1/2},\end{eqnarray}$$

and

(3.28) $$\begin{eqnarray}g_{\text{HF},\text{rms}}=\sqrt{g_{\text{rms}}^{2}-g_{\text{LF},\text{rms}}^{2}},\end{eqnarray}$$

where $\tilde{x}=x/L_{\text{box}}$ , ${\tilde{y}}=y/L_{\text{box}}$ , and $\langle \cdots \rangle$ denotes an average over $x$ , $y$ and  $t$ . The quantity

(3.29) $$\begin{eqnarray}\unicode[STIX]{x1D6E5}=c_{\text{av}}\unicode[STIX]{x1D70F}_{\text{nl}}^{\text{(tr)}}v_{\text{A}}(r_{\text{b}})\end{eqnarray}$$

is the approximate radial correlation length in the low corona of a $\boldsymbol{z}^{+}$ fluctuation that is generated by a disturbance at the transition region whose correlation time is  $\unicode[STIX]{x1D70F}_{\text{nl}}^{\text{(tr)}}$ ,

(3.30) $$\begin{eqnarray}r_{\text{i}}=\left\{\begin{array}{@{}ll@{}}r_{\text{min}} & \text{if}~r<r_{\text{min}}+\unicode[STIX]{x1D6E5},\\ r-\unicode[STIX]{x1D6E5} & \text{if}~r_{\text{min}}+\unicode[STIX]{x1D6E5}\leqslant r\leqslant r_{\text{max}}-\unicode[STIX]{x1D6E5},\\ r_{\text{max}}-2\unicode[STIX]{x1D6E5} & \text{if}~r>r_{\text{max}}-\unicode[STIX]{x1D6E5}\end{array}\right.\end{eqnarray}$$

and $c_{\text{av}}$ is a dimensionless constant of order unity. We set

(3.31) $$\begin{eqnarray}c_{\text{av}}\simeq 0.6,\end{eqnarray}$$

which enables us to carry out the radial average in (3.26) in a computationally efficient way, using an integer number of subdomains. Given the above definitions, $\unicode[STIX]{x1D6E5}=0.08R_{\odot }$ in Runs 1 and 2, and $\unicode[STIX]{x1D6E5}=0.32R_{\odot }$ in Run 3. We define

(3.32a,b ) $$\begin{eqnarray}z_{\text{LF},\text{rms}}^{+}=\frac{\unicode[STIX]{x1D702}^{1/4}g_{\text{LF},\text{rms}}}{1+\unicode[STIX]{x1D702}^{1/2}}\quad z_{\text{HF},\text{rms}}^{+}=\frac{\unicode[STIX]{x1D702}^{1/4}g_{\text{HF},\text{rms}}}{1+\unicode[STIX]{x1D702}^{1/2}}.\end{eqnarray}$$

We emphasize that, although we use the subscripts ‘LF’ and ‘HF’ as shorthand for ‘low frequency’ and ‘high frequency’, the defining difference between $z_{\text{LF},\text{rms}}^{+}$ and $z_{\text{HF},\text{rms}}^{+}$ is the difference in their radial correlation lengths.

In figure 4 we plot the radial profiles of $z_{\text{LF},\text{rms}}^{+}$ and $z_{\text{HF},\text{rms}}^{+}$ in our numerical simulations and the analytic model of § 4. As this figure shows, all three simulations contain both $\boldsymbol{z}_{\text{LF}}^{+}$ and $\boldsymbol{z}_{\text{HF}}^{+}$ fluctuations, and these fluctuations evolve in different ways as they propagate away from the Sun. In all three runs, $z_{\text{HF},\text{rms}}^{+}\simeq z_{\text{LF},\text{rms}}^{+}$ in the low corona. As $r$ increases, $z_{\text{HF},\text{rms}}^{+}/z_{\text{LF},\text{rms}}^{+}$ decreases, particularly in Run 2, suggesting that the high-frequency component of $\boldsymbol{z}^{+}$ cascades and dissipates more rapidly than the low-frequency component.

3.8 Alignment

Figure 5 shows the characteristic value of the sine of the angle between $\boldsymbol{z}^{+}$ and $\boldsymbol{z}^{-}$ ,

(3.33) $$\begin{eqnarray}\sin \unicode[STIX]{x1D703}=\frac{\langle |\boldsymbol{z}^{+}\times \boldsymbol{z}^{-}|\rangle }{\langle |\boldsymbol{z}^{+}|\rangle \langle |\boldsymbol{z}^{-}|\rangle },\end{eqnarray}$$

in both our numerical simulations and the model we present in § 4. As $r$ increases, $\sin \unicode[STIX]{x1D703}$ decreases, particularly in Run 2, causing nonlinear interactions between $\boldsymbol{z}^{+}$ and $\boldsymbol{z}^{-}$ to weaken (see, e.g. Boldyrev Reference Boldyrev2005, Reference Boldyrev2006; Perez & Chandran Reference Perez and Chandran2013; Chandran, Schekochihin & Mallet Reference Chandran, Schekochihin and Mallet2015b ).Footnote ³

Figure 5. The characteristic value of the sine of the alignment angle  $\unicode[STIX]{x1D703}$ between $\boldsymbol{z}^{+}$ and $\boldsymbol{z}^{-}$ , defined in (3.33), in Runs 1 through 3 and in the analytic model of § 4 (using (4.8)).

3.9 Turbulent heating

In figure 6 we plot the rate  $Q$ at which energy is dissipated per unit mass by hyper-dissipation in our simulations (see Perez & Chandran Reference Perez and Chandran2013) as a function of  $r$ , as well as the turbulent heating rate in the analytic model described in § 4. The amplitudes of the turbulent fluctuations in our simulations are consistent with the results of several observational studies that were summarized in figure 9 of Cranmer & van Ballegooijen (Reference Cranmer and van Ballegooijen2005), including non-thermal line widths in coronal holes inferred from SUMER (Solar Ultraviolet Measurements of Emitted Radiation) and UVCS (Ultraviolet Coronagraph Spectrometer) measurements (Banerjee et al. Reference Banerjee, Teriaca, Doyle and Wilhelm1998; Esser et al. Reference Esser, Fineschi, Dobrzycka, Habbal, Edgar, Raymond, Kohl and Guhathakurta1999). For comparison, the r.m.s. amplitudes of the fluctuating velocity  $\unicode[STIX]{x1D6FF}v_{\text{rms}}$ at $r=r_{\text{tr}}$ in Runs 1, 2 and 3 are, respectively, $30.4~\text{km}~\text{s}^{-1}$ , $30.0~\text{km}~\text{s}^{-1}$ and $26.7~\text{km}~\text{s}^{-1}$ .Footnote ⁴ The values of $\unicode[STIX]{x1D6FF}v_{\text{rms}}$ at $r=2R_{\odot }$ in Runs 1, 2 and 3 are, respectively, $170~\text{km}~\text{s}^{-1}$ , $157~\text{km}~\text{s}^{-1}$ and $146~\text{km}~\text{s}^{-1}$ . Because the turbulence amplitudes in our simulations are consistent with the aforementioned observations, the turbulent heating rate in each of our simulations can be used to estimate the rate at which transverse, non-compressive MHD turbulence would heat the solar wind as a function of  $r$ if the correlation lengths and correlation time at $r=r_{\text{b}}$ in the simulation were realistic.Footnote ⁵

Figure 6. The turbulent-heating rate per unit mass  $Q$ in Runs 1 through 3 and in the analytic model of § 4. The dotted line labelled C11 is the turbulent-heating rate in the solar-wind model of Chandran et al. (Reference Chandran, Dennis, Quataert and Bale2011), which approximates the heating needed to power the fast solar wind.

To estimate the amount of turbulent heating that would be needed to power the solar wind, we also plot in figure 6 the turbulent-heating rate in the one-dimensional (flux-tube) solar-wind model of Chandran et al. (Reference Chandran, Dennis, Quataert and Bale2011). This model included Coulomb collisions, super-radial expansion of the magnetic field, separate energy equations for the protons and electrons, proton temperature anisotropy, a transition between Spitzer conductivity near the Sun and a Hollweg collisionless heat flux at larger  $r$ and enhanced pitch-angle scattering by temperature-anisotropy instabilities in regions in which the plasma is either mirror or firehose unstable. The model agreed with a number of remote observations of coronal holes and in situ measurements of fast-solar-wind streams.

The turbulent-heating rate in the Chandran et al. (Reference Chandran, Dennis, Quataert and Bale2011) model, which we denote  $Q_{\text{C11}}$ , is for the most part comparable to (i.e. within a factor of 3 of) the heating rate in our numerical simulations. The simulated heating rates in Runs 1 and 3 are in fact strikingly close to  $Q_{\text{C11}}$ at $r\gtrsim 4R_{\odot }$ . However, in all three runs, $Q>Q_{\text{C11}}$ at $r=2-3R_{\odot }$ . This latter discrepancy is largest in the case of Run 2, in which $Q\simeq 3Q_{\text{C11}}$ at $r=2R_{\odot }$ . Although Run 2 has the largest heating rate of all three simulations at $r=2R_{\odot }$ , the simulated heating rate in Run 2 is smaller than  $Q_{\text{C11}}$ at $r\gtrsim 5R_{\odot }$ by a factor of  ${\sim}2$ .

The only region in which the simulated heating rate differs from $Q_{\text{C11}}$ by a factor ${\gtrsim}4$ is at $r<1.3R_{\odot }$ in Run 3, where $Q/Q_{\text{C11}}$ falls below 0.1. Even in Runs 1 and 2, the simulated heating rate at $r<1.3R_{\odot }$ is smaller than  $Q_{\text{C11}}$ by a factor of  ${\sim}2$ . The finding that $Q\lesssim 0.5Q_{\text{C11}}$ at $r<1.3R_{\odot }$ in all three runs may indicate the presence of additional heating mechanisms in the actual low corona, such as compressive fluctuations, a possibility previously considered by Cranmer et al. (Reference Cranmer, van Ballegooijen and Edgar2007) and Verdini et al. (Reference Verdini, Velli, Matthaeus, Oughton and Dmitruk2010).

Recently, van Ballegooijen & Asgari-Targhi (Reference van Ballegooijen and Asgari-Targhi2016), van Ballegooijen & Asgari-Targhi (Reference van Ballegooijen and Asgari-Targhi2017) carried out a series of direct numerical simulations of reflection-driven MHD turbulence and concluded that such turbulence is unable to provide enough heating to power the solar wind. The reason we reach a different conclusion is likely that we use the two-fluid solar-wind model of Chandran et al. (Reference Chandran, Dennis, Quataert and Bale2011) to estimate the amount of heating required, whereas van Ballegooijen & Asgari-Targhi (Reference van Ballegooijen and Asgari-Targhi2016), van Ballegooijen & Asgari-Targhi (Reference van Ballegooijen and Asgari-Targhi2017) used a one-fluid solar-wind model (A. van Ballegooijen, private communication). In the Chandran et al. (Reference Chandran, Dennis, Quataert and Bale2011) two-fluid model, the electron temperature is lower than the proton temperature, and thus less heat is conducted back to the chromosphere than in a one-fluid solar-wind model.

3.10 Simulation results: Elsasser power spectra

We define the perpendicular Elsasser power spectra

(3.34) $$\begin{eqnarray}E^{\pm }(k_{\bot },r)=k_{\bot }\int _{0}^{2\unicode[STIX]{x03C0}}\text{d}\unicode[STIX]{x1D719}\overline{|\tilde{\boldsymbol{z}}^{\pm }(k_{\bot },\unicode[STIX]{x1D719},r,t)|^{2}},\end{eqnarray}$$

where $\tilde{\boldsymbol{z}}^{\pm }(k_{\bot },\unicode[STIX]{x1D719})$ is the Fourier transform of  $\boldsymbol{z}^{\pm }$ in  $x$ and  $y$ (see figure 1), $\unicode[STIX]{x1D719}$ is the polar angle in the $(k_{x},k_{y})$ plane and $\overline{(\ldots \,)}$ indicates a time average. As illustrated in figure 7(a), we find that $E^{\pm }(k_{\bot })$ exhibits an approximate power-law scaling of the form

(3.35) $$\begin{eqnarray}E^{\pm }(k_{\bot },r)\propto k_{\bot }^{-\unicode[STIX]{x1D6FC}^{\pm }(r)}\end{eqnarray}$$

from $k_{\bot }\simeq 3\times 2\unicode[STIX]{x03C0}/L_{\text{box}}$ to $k_{\bot }\simeq 15\times 2\unicode[STIX]{x03C0}/L_{\text{box}}$ at all  $r$ in all three of our simulations. We evaluate $\unicode[STIX]{x1D6FC}^{\pm }(r)$ by fitting $E^{\pm }(k_{\bot },r)$ to a power law within this range of wave numbers, and plot the resulting values of $\unicode[STIX]{x1D6FC}^{\pm }(r)$ in figure 7.

Figure 7. (a) The Elsasser power spectra  $E^{\pm }(k_{\bot },r)$ defined in  (3.34) as functions of perpendicular wavenumber  $k_{\bot }$ at $r=20R_{\odot }$ in Run 1. (b,c,d) The spectral indices $\unicode[STIX]{x1D6FC}^{+}(r)$ and  $\unicode[STIX]{x1D6FC}^{-}(r)$ defined in (3.35) in our three numerical simulations.

Although we drive only large-scale ( $k_{\bot }\leqslant 3\times 2\unicode[STIX]{x03C0}/L_{\text{box}}$ ) velocity fluctuations at the photosphere, figure 7 shows that there is broad-spectrum turbulence throughout the chromosphere. This is because of the strong reflection of $\boldsymbol{z}^{+}$ fluctuations at the transition region and the strong reflection of  $\boldsymbol{z}^{-}$ fluctuations (at all  $k_{\bot }$ ) at the photosphere (enforced by the fixed-velocity boundary condition at $r=R_{\odot }$ ), which together lead to comparatively ‘balanced’ turbulence (meaning that $z_{\text{rms}}^{+}\simeq z_{\text{rms}}^{-}$ ) in the chromosphere, as shown previously by van Ballegooijen et al. (Reference van Ballegooijen, Asgari-Targhi, Cranmer and DeLuca2011). In the low chromosphere, $\unicode[STIX]{x1D6FC}^{\pm }\simeq 1.3{-}1.5$ in all four simulations, which is similar to the value $\unicode[STIX]{x1D6FC}^{\pm }\simeq 3/2$ that arises in numerical simulations of homogeneous, balanced, RMHD turbulence (Mason, Cattaneo & Boldyrev Reference Mason, Cattaneo and Boldyrev2008; Beresnyak Reference Beresnyak2012; Perez et al. Reference Perez, Mason, Boldyrev and Cattaneo2012). On the other hand, $\unicode[STIX]{x1D6FC}^{+}$ decreases from $\simeq 1.5$ to $\simeq 0.8$ as $r$ increases from $1.001R_{\odot }$ to $r_{\text{tr}}=1.0026R_{\odot }$ in Runs 1 and 2. This spectral flattening arises because the Alfvén-speed gradient in the upper chromosphere acts as a high-pass filter on outward-propagating AWs in Runs 1 and 2, causing lower- $k_{\bot }$ (and hence lower-frequency – see Goldreich & Sridhar (Reference Goldreich and Sridhar1995)) $\boldsymbol{z}^{+}$ fluctuations to undergo non-WKB reflection, and allowing higher- $k_{\bot }$ (and hence higher-frequency) $\boldsymbol{z}^{+}$ fluctuations to propagate unhindered to the transition region (Velli Reference Velli1993; Réville, Tenerani & Velli Reference Réville, Tenerani and Velli2018). The difference in Run 3 is that $L_{\text{box}}$ is larger, and thus $\boldsymbol{z}^{+}$ fluctuations do not reach sufficiently large $k_{\bot }$ values that they can avoid non-WKB reflection in the upper chromosphere.Footnote ⁶ The idea that $\boldsymbol{z}^{+}$ fluctuations at the high- $k_{\bot }$ end of the inertial range propagate through the chromosphere more easily in Runs 1 and 2 than in Run 3 is consistent with the fact that $z_{\text{rms}}^{+}(r_{\text{b}})$ and $\unicode[STIX]{x1D6FF}B_{\text{rms}}$ are somewhat larger in Runs 1 and 2 than in Run 3 (see table 1 and figure 3).

As $r$ increases from $r_{\text{tr}}$ to $r_{\text{A}}$ and beyond, $\unicode[STIX]{x1D6FC}^{\pm }$ approaches $\simeq 3/2$ in all three runs. In Runs 1 and 2, the increase in $\unicode[STIX]{x1D6FC}^{+}$ as $r$ increases from $r_{\text{tr}}$ to $r_{\text{A}}$ is not steady. In Run 1, $\unicode[STIX]{x1D6FC}^{+}$ decreases as $r$ increases from $2.8R_{\odot }$ to $4.2R_{\odot }$ , and in Run 2, $\unicode[STIX]{x1D6FC}^{+}$ plateaus around a value of 1 between $r=2R_{\odot }$ and $r=3R_{\odot }$ . This behaviour suggests that, in these two simulations, the turbulent dynamics at $2R_{\odot }\lesssim r\lesssim 4R_{\odot }$ drives $\unicode[STIX]{x1D6FC}^{+}$ towards a value close to 1, and the tendency for $\unicode[STIX]{x1D6FC}^{+}$ to evolve towards 3/2 sets in at $r\gtrsim 4R_{\odot }$ . We discuss these trends further in § 6.

4.1 Amplitude of the inward-propagating fluctuations

To determine $z_{\text{rms}}^{-}$ , we assume that $\boldsymbol{z}^{-}$ cascades primarily via interactions with $\boldsymbol{z}_{\text{LF}}^{+}$ (or, equivalently, $\boldsymbol{g}_{\text{LF}}$ ). The outer-scale $\boldsymbol{z}^{-}$ cascade rate then depends upon the critical-balance parameter (Goldreich & Sridhar Reference Goldreich and Sridhar1995; Boldyrev Reference Boldyrev2006)

(4.9) $$\begin{eqnarray}\unicode[STIX]{x1D712}_{\text{LF}}^{-}=\frac{z_{\text{LF},\text{rms}}^{+}L_{r,\text{LF}}^{+}A}{L_{\bot }v_{\text{A}}},\end{eqnarray}$$

where $L_{r,\text{LF}}^{+}$ is the radial correlation length of the $\boldsymbol{g}_{\text{LF}}$ fluctuations. The critical-balance parameter $\unicode[STIX]{x1D712}_{\text{LF}}^{-}$ is an estimate of the fractional change in an outer-scale $\boldsymbol{z}^{-}$ fluctuation that results from a single ‘collision’ with an outer-scale $\boldsymbol{g}_{\text{LF}}$ fluctuation lasting a time $\unicode[STIX]{x0394}t\sim L_{r,\text{LF}}^{+}/v_{\text{A}}$ (Lithwick, Goldreich & Sridhar Reference Lithwick, Goldreich and Sridhar2007).

If $\unicode[STIX]{x1D712}_{\text{LF}}^{-}\ll 1$ , then each such collision causes only a small perturbation to the outer-scale $\boldsymbol{z}^{-}$ fluctuation, and the turbulence is weak. In this limit, the effects of successive collisions add like a random walk, and roughly

(4.10) $$\begin{eqnarray}N=(\unicode[STIX]{x1D712}_{\text{LF}}^{-})^{-2}\end{eqnarray}$$

collisions are needed for nonlinear interactions to cause an order-unity change in the outer-scale $\boldsymbol{z}^{-}$ fluctuation. The outer-scale $\boldsymbol{z}^{-}$ cascade time scale  $t_{\text{NL}}^{-}$ is then  ${\sim}NL_{r,\text{LF}}^{+}/v_{\text{A}}$ . The generation of outer-scale $\boldsymbol{z}^{-}$ (or $\boldsymbol{f}$ ) fluctuations by non-WKB reflection in this weak-turbulence regime can also be viewed as a random-walk-like process. Equation (2.20) implies that, in a reference frame  $S^{-}$ that propagates with radial velocity  $U-v_{\text{A}}$ , the increment to $\boldsymbol{f}$ from non-WKB reflection during a time $\unicode[STIX]{x0394}t=L_{r,\text{LF}}^{+}/v_{\text{A}}$ is of order

(4.11) $$\begin{eqnarray}\unicode[STIX]{x0394}f\sim \left(\frac{U-v_{\text{A}}}{2v_{\text{A}}}\right)\left|\frac{\text{d}v_{\text{A}}}{\text{d}r}\right|g_{\text{LF},\text{rms}}\unicode[STIX]{x0394}t.\end{eqnarray}$$

It follows from (2.15) that the corresponding increment to $\boldsymbol{z}^{-}$ is of order

(4.12) $$\begin{eqnarray}\unicode[STIX]{x0394}z^{-}\sim \left(\frac{U+v_{\text{A}}}{2v_{\text{A}}}\right)\left|\frac{\text{d}v_{\text{A}}}{\text{d}r}\right|z_{\text{LF},\text{rms}}^{+}\unicode[STIX]{x0394}t.\end{eqnarray}$$

The r.m.s. value of $\boldsymbol{z}^{-}$ is approximately the ‘amount’ of $\boldsymbol{z}^{-}$ that ‘builds up’ in frame  $S^{-}$ by non-WKB reflection during the cascade/damping time scale ${\sim}N\unicode[STIX]{x0394}t$ . The resulting value of  $z_{\text{rms}}^{-}$ is ${\sim}N^{1/2}\unicode[STIX]{x0394}z^{-}$ , or, equivalently,

(4.13) $$\begin{eqnarray}z_{\text{rms}}^{-}\sim \frac{L_{\bot }}{A}\left(\frac{U+v_{\text{A}}}{2v_{\text{A}}}\right)\left|\frac{\text{d}v_{\text{A}}}{\text{d}r}\right|.\end{eqnarray}$$

If $\unicode[STIX]{x1D712}_{\text{LF}}^{-}\gtrsim 1$ , then the outer-scale $\boldsymbol{z}^{-}$ fluctuations are sheared coherently throughout their lifetimes, the turbulence is strong and $t_{\text{NL}}^{-}\sim L_{\bot }/(z_{\text{LF},\text{rms}}^{+}A)$ . In this case, $z_{\text{rms}}^{-}$ is approximately the rate at which $\boldsymbol{z}^{-}$ fluctuations are produced by non-WKB reflection multiplied by $t_{\text{NL}}^{-}$ , which again leads to (4.13). This estimate, with $A\rightarrow 1$ , is the same as that obtained by Chandran & Hollweg (Reference Chandran and Hollweg2009) for the strong-turbulence limit. In the limits $U\rightarrow 0$ and $A\rightarrow 1$ , equation (4.13) is also the same as the estimate by Dmitruk et al. (Reference Dmitruk, Matthaeus, Milano, Oughton, Zank and Mullan2002) for the strong-turbulence limit.

Since (4.13) holds in both the weak and strong-turbulence regimes, we set

(4.14) $$\begin{eqnarray}z_{\text{rms}}^{-}=\frac{c^{-}L_{\bot }}{A}\left(\frac{U+v_{\text{A}}}{2v_{\text{A}}}\right)\left|\frac{\text{d}v_{\text{A}}}{\text{d}r}\right|\end{eqnarray}$$

regardless of the value of  $\unicode[STIX]{x1D712}_{\text{LF}}^{-}$ , where $c^{-}$ is a dimensionless piecewise constant function that has one value at $r>r_{\text{m}}$ and a smaller value at $r\leqslant r_{\text{m}}$ , where $r_{\text{m}}$ is defined in (3.12). Before discussing $c^{-}$ further, we note an immediate consequence of (4.14), that $z_{\text{rms}}^{-}$ increases as  $A$ (or equivalently  $\sin \unicode[STIX]{x1D703}$ ) decreases. This is because reducing  $\sin \unicode[STIX]{x1D703}$ decreases the rate at which outer-scale $\boldsymbol{z}^{-}$ fluctuations cascade without decreasing the rate at which they are produced by non-WKB reflection.

4.2 Suppression of inward-propagating fluctuations at $r<r_{\text{m}}$

The reason we take $c^{-}$ to have a smaller value at $r<r_{\text{m}}$ than at $r>r_{\text{m}}$ is that the non-WKB-reflection source term for $\boldsymbol{z}^{-}$ fluctuations reverses direction at $r=r_{\text{m}}$ , since $\text{d}v_{\text{A}}/\text{d}r$ changes sign. Since $\boldsymbol{g}_{\text{LF}}$ has a large radial correlation length, when $\boldsymbol{z}^{-}$ fluctuations produced via non-WKB reflection at $r>r_{\text{m}}$ propagate to $r<r_{\text{m}}$ , they tend to cancel out the $\boldsymbol{z}^{-}$ fluctuations that are produced via non-WKB reflection at $r<r_{\text{m}}$ , reducing  $z_{\text{rms}}^{-}$ . If the $\boldsymbol{z}^{-}$ fluctuations at $r=r_{\text{m}}$ can propagate a radial distance ${\sim}(r_{\text{m}}-R_{\odot })$ before cascading and dissipating, then this cancellation effect is large. On the other hand, if the $\boldsymbol{z}^{-}$ fluctuations at $r=r_{\text{m}}$ can only propagate a radial distance  $\ll (r_{\text{m}}-R_{\odot })$ before cascading and dissipating, then little cancellation occurs. To account for this phenomenology, we set

(4.15) $$\begin{eqnarray}c^{-}=\left\{\begin{array}{@{}ll@{}}c_{\text{I}}^{-}\quad & \text{if}~r\leqslant r_{\text{m}},\\ c_{\text{O}}^{-}\quad & \text{if}~r>r_{\text{m}},\end{array}\right.\end{eqnarray}$$

where $c_{\text{O}}^{-}$ is a dimensionless free parameter,

(4.16) $$\begin{eqnarray}\displaystyle & \displaystyle c_{\text{I}}^{-}=\frac{c_{\text{O}}^{-}}{1+M}, & \displaystyle\end{eqnarray}$$

(4.17) $$\begin{eqnarray}\displaystyle & \displaystyle M=\frac{v_{\text{A}\text{m}}L_{\bot m}}{z_{\text{LFm}}^{+}L_{\unicode[STIX]{x1D735}}}, & \displaystyle\end{eqnarray}$$

$L_{\unicode[STIX]{x1D735}}$ is a free parameter with dimensions of length, and $v_{\text{Am}}$ , $L_{\bot m}$ , and $z_{\text{LF}\text{m}}^{+}$ are the values of $v_{\text{A}}$ , $L_{\bot }$ and $z_{\text{LF},\text{rms}}^{+}$ at $r=r_{\text{m}}$ . As we argue below (see (4.20)), $A$ is of order unity at $r<r_{\text{m}}$ , which means that $M$ is the approximate radial distance an outer-scale $\boldsymbol{z}^{-}$ fluctuation at $r=r_{\text{m}}$ propagates before cascading to smaller scales, divided by  $L_{\unicode[STIX]{x1D735}}$ . We can rewrite  $M$ in terms of quantities evaluated at $r=r_{\text{b}}$ by making the approximations that $v_{\text{A}}\gg U$ at $r\leqslant r_{\text{m}}$ and $g_{\text{LF},\text{rms}}(r_{\text{m}})\simeq g_{\text{LF},\text{rms}}(r_{\text{b}})$ and by using (2.17). This yields

(4.18) $$\begin{eqnarray}M=\frac{v_{\text{Ab}}L_{\bot \text{b}}}{z_{\text{LF}\text{b}}^{+}L_{\unicode[STIX]{x1D735}}}\left(\frac{v_{\text{A}\text{m}}}{v_{\text{Ab}}}\right)^{1/2},\end{eqnarray}$$

where $v_{\text{Ab}}$ , $L_{\bot b}$ and $z_{\text{LF}\text{b}}^{+}$ are the values of $v_{\text{A}}$ , $L_{\bot }$ and $z_{\text{LF},\text{rms}}^{+}$ at $r=r_{\text{b}}$ .

4.3 Alignment factor and critical-balance parameter

To estimate the alignment factor  $A$ introduced in (4.7), we first note that nonlinear interactions between counter-propagating AWs produce negative residual energy, with $\boldsymbol{z}^{-}$ anti-parallel to  $\boldsymbol{z}^{+}$ (i.e. an excess of magnetic energy over kinetic energy) (Müller & Grappin Reference Müller and Grappin2005; Boldyrev et al. Reference Boldyrev, Perez, Borovsky and Podesta2011). At $r>r_{\text{m}}$ , $\text{d}v_{\text{A}}/\text{d}r<0$ , and it follows from (2.15) and (2.20) that non-WKB reflection also acts to produce negative residual energy. On the other hand, at $r<r_{\text{m}}$ , $\text{d}v_{\text{A}}/\text{d}r>0$ , and non-WKB reflection acts to produce positive residual energy. In other words, at $r<r_{\text{m}}$ , linear processes (non-WKB reflection) and nonlinear processes have competing effects on the alignment of  $\boldsymbol{z}^{-}$ . Based on these arguments, we conjecture that at $r<r_{\text{m}}$ the outer-scale fluctuations do not develop significant alignment, and that at $r>r_{\text{m}}$ the outer-scale $\boldsymbol{z}_{\text{LF}}^{+}$ and $\boldsymbol{z}^{-}$ fluctuations become increasingly aligned as the $\boldsymbol{z}_{\text{LF}}^{+}$ fluctuations ‘decay’ via nonlinear interactions. We also conjecture that $A$ is a decreasing function of $\unicode[STIX]{x1D70F}_{v}^{\text{(ph)}}$ , because a larger $\unicode[STIX]{x1D70F}_{v}^{\text{(ph)}}$ increases the efficiency of non-WKB reflection, which produces $\boldsymbol{z}^{-}$ fluctuations that are aligned with  $\boldsymbol{z}_{\text{LF}}^{+}$ . In addition, we conjecture that $A$ is a decreasing function of

(4.19) $$\begin{eqnarray}\unicode[STIX]{x1D6E4}=\frac{z_{\text{LF},\text{rms}}^{+}L_{r,\text{LF}}^{+}}{L_{\bot }v_{\text{A}}},\end{eqnarray}$$

which is the critical-balance parameter  $\unicode[STIX]{x1D712}_{\text{LF}}^{-}$ in (4.9) without the factor of $A$ . There are two reasons for taking $A$ to decrease with increasing  $\unicode[STIX]{x1D6E4}$ . The first is that when $\unicode[STIX]{x1D6E4}\ll 1$ , outer-scale $\boldsymbol{z}^{-}$ fluctuations can propagate through many different outer-scale $\boldsymbol{z}_{\text{LF}}^{+}$ fluctuations before cascading to smaller scales. The $\boldsymbol{z}^{-}$ fluctuations that are co-located with a particular outer-scale $\boldsymbol{z}_{\text{LF}}^{+}$ ‘eddy’ of radial extent  ${\sim}L_{r,\text{LF}}^{+}$ are thus a mixture of the $\boldsymbol{z}^{-}$ fluctuations produced by the non-WKB reflection of that $\boldsymbol{z}_{\text{LF}}^{+}$ eddy and $\boldsymbol{z}^{-}$ fluctuations that were initially produced by the non-WKB reflection of $\boldsymbol{z}_{\text{LF}}^{+}$ eddies located farther from the Sun. The greater the number of distinct outer-scale $\boldsymbol{z}_{\text{LF}}^{+}$ eddies whose reflections contribute to the value of $\boldsymbol{z}^{-}$ at any single point, the less aligned the $\boldsymbol{z}^{-}$ field will be with any individual $\boldsymbol{z}_{\text{LF}}^{+}$ eddy. Moreover, when $\unicode[STIX]{x1D6E4}>1$ , shearing of the $\boldsymbol{z}^{-}$ fluctuations by $\boldsymbol{z}_{\text{LF}}^{+}$ rotates the $\boldsymbol{z}^{-}$ fluctuations into alignment with $\boldsymbol{z}_{\text{LF}}^{+}$ , and the resulting value of  $A$ is a decreasing function of  $\unicode[STIX]{x1D6E4}$ (Chandran et al. Reference Chandran, Schekochihin and Mallet2015b ). We quantify the foregoing conjectures by setting

(4.20) $$\begin{eqnarray}A=\left\{\begin{array}{@{}ll@{}}A_{0}\quad & \text{if }r<r_{\text{m}},\\ \displaystyle A_{0}\left[1+\frac{\unicode[STIX]{x1D70F}_{v}^{\text{(ph)}}\unicode[STIX]{x1D6E4}}{\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D703}}}\ln \left(\frac{g_{\text{LFm}}^{2}}{g_{\text{LF},\text{rms}}^{2}}\right)\right]^{-1}\quad & \text{if}~r>r_{\text{m}},\end{array}\right.\end{eqnarray}$$

where the dimensionless constant  $A_{0}$ and the time constant $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D703}}$ are free parameters.

In the linear, short-wavelength, AW propagation problem, if an AW is launched into a coronal hole by a boundary condition imposed at the transition region and photosphere, and if the AW period is  $P$ , then the radial wavelength of the AW at radius  $r$ is $(U+v_{\text{A}})P$ . That is, the wave period remains constant as the wave propagates away from the Sun, and the radial wavelength varies in proportion to the wave phase velocity. We take nonlinear, non-WKB $\boldsymbol{z}^{+}$ fluctuations to behave in the same way, setting

(4.21) $$\begin{eqnarray}\frac{L_{r,\text{LF}}^{+}}{L_{r,\text{LFb}}^{+}}=\frac{U+v_{\text{A}}}{U_{\text{b}}+v_{\text{Ab}}},\end{eqnarray}$$

where $L_{r,\text{LFb}}^{+}$ , $U_{\text{b}}$ and $v_{\text{Ab}}$ are the values, respectively, of $L_{r,\text{LF}}^{+}$ , $U$ and  $v_{\text{A}}$ evaluated at $r=r_{\text{b}}$ , and likewise for $L_{r,\text{HF}}^{+}$ . It then follows from (2.17), (3.3), (3.23) and (4.21) that

(4.22) $$\begin{eqnarray}\unicode[STIX]{x1D6E4}=\unicode[STIX]{x1D6E4}_{\text{b}}\frac{g_{\text{LF},\text{rms}}}{g_{\text{LF}\text{b}}}\sqrt{\frac{v_{\text{A}}}{v_{\text{Ab}}}},\end{eqnarray}$$

where $\unicode[STIX]{x1D6E4}_{\text{b}}$ is the value of  $\unicode[STIX]{x1D6E4}$ at $r=r_{\text{b}}$ .

4.4 Solving for the fluctuation-amplitude profiles

Upon combining (4.6), (4.7), (4.14) and (4.15), we obtain

(4.23) $$\begin{eqnarray}\frac{\text{d}}{\text{d}r}\ln g_{\text{LF},\text{rms}}^{2}=\left\{\begin{array}{@{}ll@{}}\displaystyle -c_{\text{I}}\frac{\text{d}}{\text{d}r}\ln v_{\text{A}}\quad & \text{if}~r\leqslant r_{\text{m}},\\ \displaystyle c_{\text{O}}\frac{\text{d}}{\text{d}r}\ln v_{\text{A}}\quad & \text{if}~r>r_{\text{m}},\end{array}\right.\end{eqnarray}$$

where

(4.24a,b ) $$\begin{eqnarray}c_{\text{I}}\equiv c_{\text{diss}}c_{\text{I}}^{-}\quad c_{\text{O}}\equiv c_{\text{diss}}c_{\text{O}}^{-}.\end{eqnarray}$$

After integrating (4.23), we find that

(4.25) $$\begin{eqnarray}\frac{g_{\text{LF},\text{rms}}^{2}}{g_{\text{LFb}}^{2}}=\left\{\begin{array}{@{}ll@{}}\displaystyle \left(\frac{v_{\text{Ab}}}{v_{\text{A}}}\right)^{c_{\text{I}}} & \text{if}~r_{\text{b}}<r<r_{\text{m}},\\ \displaystyle \left(\frac{v_{\text{Ab}}}{v_{\text{A}\text{m}}}\right)^{c_{\text{I}}}\left(\frac{v_{\text{A}}}{v_{\text{A}\text{m}}}\right)^{c_{\text{O}}} & \text{if}~r>r_{\text{m}},\end{array}\right.,\end{eqnarray}$$

where $g_{\text{LF}\text{b}}$ is the value of $g_{\text{LF},\text{rms}}$ at $r=r_{\text{b}}$ . Upon combining (4.4), (4.5), (4.14) and (4.15), we obtain

(4.26) $$\begin{eqnarray}\frac{\text{d}}{\text{d}r}\ln g_{\text{HF},\text{rms}}^{2}=\left\{\begin{array}{@{}ll@{}}\displaystyle -\frac{c_{\text{I}}}{A}\frac{\text{d}}{\text{d}r}\ln v_{\text{A}} & \text{if}~r\leqslant r_{\text{m}}\\ \displaystyle \frac{c_{\text{O}}}{A}\frac{\text{d}}{\text{d}r}\ln v_{\text{A}} & \text{if}~r>r_{\text{m}}.\end{array}\right.\end{eqnarray}$$

With the aid of (4.20) and (4.22), we integrate (4.26) to obtain

(4.27) $$\begin{eqnarray}\frac{g_{\text{HF},\text{rms}}^{2}}{g_{\text{HFb}}^{2}}=\left\{\begin{array}{@{}ll@{}}\displaystyle \left(\frac{v_{\text{Ab}}}{v_{\text{A}}}\right)^{c_{\text{I}}/A_{0}} & \text{if}~r<r_{\text{m}},\\ \displaystyle \left(\frac{v_{\text{Ab}}}{v_{\text{Am}}}\right)^{c_{\text{I}}/A_{0}}w^{c_{\text{O}}/A_{0}}\text{e}^{-H} & \text{if}~r>r_{\text{m}},\end{array}\right.\end{eqnarray}$$

where

(4.28) $$\begin{eqnarray}H=\frac{c_{\unicode[STIX]{x1D703}}c_{\text{O}}^{2}\unicode[STIX]{x1D6E4}_{\text{b}}}{\unicode[STIX]{x1D70E}^{2}A_{0}}\left(\frac{v_{\text{Am}}}{v_{\text{Ab}}}\right)^{(1-c_{\text{I}})/2}(\unicode[STIX]{x1D70E}w^{\unicode[STIX]{x1D70E}}\ln w-w^{\unicode[STIX]{x1D70E}}+1),\end{eqnarray}$$

(4.29a,b ) $$\begin{eqnarray}w=\frac{v_{\text{A}}}{v_{\text{A}\text{m}}}\quad \unicode[STIX]{x1D70E}=\frac{1+c_{\text{O}}}{2},\end{eqnarray}$$

$g_{\text{HFb}}$ is the value of $g_{\text{HF},\text{rms}}$ at $r=r_{\text{b}}$ and $c_{\unicode[STIX]{x1D703}}=\unicode[STIX]{x1D70F}_{v}^{\text{(ph)}}/\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D703}}$ .

4.5 Turbulent-heating rate

The turbulent-heating rate in our model is

(4.30) $$\begin{eqnarray}Q=\frac{\unicode[STIX]{x1D70C}}{2}\left[\unicode[STIX]{x1D6FE}_{\text{HF}}^{+}(z_{\text{HF},\text{rms}}^{+})^{2}+\unicode[STIX]{x1D6FE}_{\text{LF}}^{+}(z_{\text{LF},\text{rms}}^{+})^{2}+\unicode[STIX]{x1D6FE}^{-}(z_{\text{rms}}^{-})^{2}\right],\end{eqnarray}$$

where

(4.31) $$\begin{eqnarray}\unicode[STIX]{x1D6FE}^{-}=\unicode[STIX]{x1D6FE}_{\text{LF}}^{-}+\unicode[STIX]{x1D6FE}_{\text{HF}}^{-}\end{eqnarray}$$

is the cascade rate of the outer-scale $\boldsymbol{z}^{-}$ fluctuations, and $\unicode[STIX]{x1D6FE}_{\text{LF}}^{-}$ ( $\unicode[STIX]{x1D6FE}_{\text{HF}}^{-}$ ) is the contribution to $\unicode[STIX]{x1D6FE}^{-}$ from interactions between $\boldsymbol{z}^{-}$ fluctuations and LF (HF) $\boldsymbol{z}^{+}$ fluctuations. To allow for either weakly turbulent ( $\unicode[STIX]{x1D712}_{\text{LF}}^{-}<1$ ) or strongly turbulent ( $\unicode[STIX]{x1D712}_{\text{LF}}^{-}\geqslant 1$ ) shearing of $\boldsymbol{z}^{-}$ fluctuations by $\boldsymbol{z}_{\text{LF}}^{+}$ fluctuations, we set

(4.32) $$\begin{eqnarray}\unicode[STIX]{x1D6FE}_{\text{LF}}^{-}=\frac{c_{\text{diss}}z_{\text{LF},\text{rms}}^{+}A}{L_{\bot }}\times \left\{\begin{array}{@{}ll@{}}\unicode[STIX]{x1D712}_{\text{LF}}^{-} & \text{if}~\unicode[STIX]{x1D712}_{\text{LF}}^{-}\leqslant 1,\\ 1 & \text{if}~\unicode[STIX]{x1D712}_{\text{LF}}^{-}>1.\end{array}\right.\end{eqnarray}$$

In analogy to (4.9), we define the critical-balance parameter for the shearing of $\boldsymbol{z}^{-}$ fluctuations by $\boldsymbol{z}_{\text{HF}}^{+}$ fluctuations to be

(4.33) $$\begin{eqnarray}\unicode[STIX]{x1D712}_{\text{HF}}^{-}=\frac{z_{\text{HF},\text{rms}}^{+}L_{r,\text{HF}}^{+}}{L_{\bot }v_{\text{A}}},\end{eqnarray}$$

where we have omitted the factor of  $A$ , because we take $\boldsymbol{z}^{-}$ to be aligned with $\boldsymbol{z}_{\text{LF}}^{+}$ but not with $\boldsymbol{z}_{\text{HF}}^{+}$ . We then set

(4.34) $$\begin{eqnarray}\unicode[STIX]{x1D6FE}_{\text{HF}}^{-}=\frac{c_{\text{diss}}z_{\text{HF},\text{rms}}^{+}}{L_{\bot }}\times \left\{\begin{array}{@{}ll@{}}\unicode[STIX]{x1D712}_{\text{HF}}^{-} & \text{if}~\unicode[STIX]{x1D712}_{\text{HF}}^{-}\leqslant 1,\\ 1 & \text{if}~\unicode[STIX]{x1D712}_{\text{HF}}^{-}>1,\end{array}\right.\end{eqnarray}$$

4.6 Comparison with simulation results

To compare our model with one of our numerical simulations, we treat $L_{\bot }(r_{\text{b}})=L_{\text{box}}(r_{\text{b}})/3$ , $L_{r,\text{LF}}^{+}(r_{\text{b}})$ , $L_{r,\text{HF}}^{+}(r_{\text{b}})$ and $z_{\text{rms}}^{+}(r_{\text{b}})$ as boundary conditions in our model, which we determine using the measured values of these quantities in that particular simulation. Also, motivated by figure 4, we set

(4.35) $$\begin{eqnarray}\frac{g_{\text{HFb}}^{2}}{g_{\text{LFb}}^{2}}=1\end{eqnarray}$$

in all our model solutions. We take $L_{r,\text{HF}}^{+}(r_{\text{b}})$ in our simulations to be the radial separation  $\unicode[STIX]{x0394}r$ at which $C(r_{\text{b}},\unicode[STIX]{x0394}r)=1/2$ , where

(4.36) $$\begin{eqnarray}C(r,\unicode[STIX]{x0394}r)=\frac{\langle \boldsymbol{g}(\tilde{x},{\tilde{y}},r,t)\boldsymbol{\cdot }\boldsymbol{g}(\tilde{x},{\tilde{y}},r+\unicode[STIX]{x0394}r,t)\rangle }{\langle |\boldsymbol{g}(\tilde{x},{\tilde{y}},r,t)|^{2}\rangle }\end{eqnarray}$$

is the radial autocorrelation function of the  $\boldsymbol{g}$ fluctuations and $\tilde{x}$ and ${\tilde{y}}$ are defined following (3.28). On the other hand, because figure 4 shows that the LF fluctuations are energetically dominant at $r=r_{\text{max}}$ , we define $L_{r,\text{LF}}^{+}(r_{\text{max}})$ to be the value of  $\unicode[STIX]{x0394}r$ at which $C(r_{\text{max}},-\unicode[STIX]{x0394}r)=1/2$ . Applying (4.21), we then set $L_{r,\text{LF}}^{+}(r_{\text{b}})=L_{r,\text{LF}}^{+}(r_{\text{max}})(U_{\text{b}}+v_{\text{Ab}})/[U(r_{\text{max}})+v_{\text{A}}(r_{\text{max}})]=0.886L_{\text{LF}}^{+}(r_{\text{max}})$ . The values of $L_{\text{box}}(r_{\text{b}})$ , $L_{r,\text{LF}}^{+}(r_{\text{b}})$ , $L_{r,\text{HF}}^{+}(r_{\text{b}})$ and $z_{\text{rms}}^{+}(r_{\text{b}})$ in our three simulations are listed in table 3.

Table 3. Boundary conditions in our analytic model for matching Runs 1 through 3.

Note: $z_{\text{rms}}^{+}$ is the r.m.s. amplitude of the outward-propagating Elsasser variable, $L_{r,\text{LF}}^{+}$ is the radial correlation length of the low-frequency outward-propagating Heinemann–Olbert variable $\boldsymbol{g}_{\text{LF}}$ , $L_{r,\text{HF}}^{+}$ is the radial correlation length of the high-frequency outward-propagating Heinemann–Olbert variable $\boldsymbol{g}_{\text{HF}}$ , $L_{\bot }$ is the perpendicular outer scale (see (3.23)) and $r_{\text{b}}$ is the radius of the coronal base defined in (3.10).

Table 4. Best-fit free parameters in our analytic model.

Note: The quantity $c_{\text{diss}}$ is a coefficient appearing in the cascade/damping rates  $\unicode[STIX]{x1D6FE}_{\text{HF}}^{+}$ and $\unicode[STIX]{x1D6FE}_{\text{LF}}^{+}$ ((4.5) and (4.7)),  $c_{\text{O}}^{-}$ is a coefficient in our estimate of  $z_{\text{rms}}^{-}$ (see (4.14) through (4.16)), $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D703}}$ and $A_{0}$ are constants appearing in our estimate of the alignment angle (4.20) and $L_{\unicode[STIX]{x1D735}}$ is a length scale that affects the degree to which $\boldsymbol{z}^{-}$ fluctuations produced by non-WKB reflection at $r>r_{\text{m}}\simeq 1.7R_{\odot }$ cancel out the $\boldsymbol{z}^{-}$ fluctuations produced by non-WKB reflection at $r<r_{\text{m}}$ (see (4.14) through (4.17)).

We take the free parameters $c_{\text{diss}}$ , $c_{\text{O}}^{-}$ , $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D703}}$ , $A_{0}$ and  $L_{\unicode[STIX]{x1D735}}$ to be the same regardless of the simulation with which we are comparing our model. We then vary these free parameters to optimize the agreement between our model and all three simulations. We list the resulting parameters in table 4.

Figures 3 through 6 show the radial profiles of $z_{\text{rms}}^{+}$ , $z_{\text{rms}}^{-}$ , $z_{\text{HF},\text{rms}}^{+}$ , $z_{\text{LF},\text{rms}}^{+}$ , $\sin \unicode[STIX]{x1D703}$ and  $Q$ that result from our model using the best-fit parameters in table 4 and the boundary conditions in table 3. As these figures show, our model reproduces a number of trends seen in the simulations. For example, in both the model and simulations, $z_{\text{HF},\text{rms}}^{+}/z_{\text{LF},\text{rms}}^{+}$ and $\sin \unicode[STIX]{x1D703}$ decrease with increasing  $r$ , particularly in Run 2. The radial decrease in $z_{\text{HF}}^{+}/z_{\text{LF}}^{+}$ is consistent with our expectation that high-frequency $\boldsymbol{z}^{+}$ fluctuations cascade and dissipate more rapidly than low-frequency $\boldsymbol{z}^{+}$ fluctuations, because high-frequency $\boldsymbol{z}^{+}$ fluctuations are not aligned with $\boldsymbol{z}^{-}$ . In our model, the radial decrease in $\sin \unicode[STIX]{x1D703}$ is related both to the comparatively rapid cascade of the unaligned high-frequency $\boldsymbol{z}^{+}$ fluctuations and the fact that the low-frequency $\boldsymbol{z}^{+}$ fluctuations become increasingly aligned with  $\boldsymbol{z}^{-}$ as they interact nonlinearly with  $\boldsymbol{z}^{-}$ . We note that the decrease in $\sin \unicode[STIX]{x1D703}$ coincides with an increase in $z_{\text{rms}}^{-}$ for the reasons described following (4.14). The model reproduces the $z_{\text{rms}}^{\pm }$ profiles in the simulations fairly accurately. The turbulent-heating rates in the model and simulations also agree quite well, but the heating rate in the model is somewhat smaller than in Run 3 at $r>3R_{\odot }$ . The most notable failing of the model is that $z_{\text{rms}}^{-}=Q=0$ at $r=r_{\text{m}}$ , because our estimate of $z_{\text{rms}}^{-}$ is proportional to the local value of $\text{d}v_{\text{A}}/\text{d}r$ , which vanishes at $r=r_{\text{m}}$ . A more realistic model would account for the fact that $\boldsymbol{z}^{-}$ fluctuations propagate a finite distance before cascading and dissipating, which would smooth out the profiles of $z_{\text{rms}}^{-}$ and  $Q$ in the vicinity of $r=r_{\text{m}}$ . Importantly, despite the aforementioned differences between the model and our numerical results, varying the boundary conditions in the model to match the measured conditions in the simulations largely accounts for the differences between the $z_{\text{rms}}^{\pm }$ and $Q$ profiles in Runs 1, 2, and 3 without any modification to the free parameters in table 4. This suggests that the model provides a reasonably accurate representation of the dominant physical processes that control these radial profiles.

5.1 Balanced, homogeneous RMHD turbulence

In ‘balanced turbulence’, the statistical properties of $\boldsymbol{z}^{+}$ and $\boldsymbol{z}^{-}$ fluctuations are identical. In particular,

(5.1a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D6FF}z_{\unicode[STIX]{x1D706}}^{+}=\unicode[STIX]{x1D6FF}z_{\unicode[STIX]{x1D706}}^{-}\quad l_{\unicode[STIX]{x1D706}}^{+}=l_{\unicode[STIX]{x1D706}}^{-},\end{eqnarray}$$

and the cross-helicity (the difference between the energies per unit mass of $\boldsymbol{z}^{+}$ and $\boldsymbol{z}^{-}$ fluctuations) is zero. In homogeneous RMHD turbulence, the strongest nonlinear interactions are local in scale, meaning that $\unicode[STIX]{x1D6FF}z_{\unicode[STIX]{x1D706}}^{\pm }$ fluctuations are cascaded primarily by $\boldsymbol{z}^{\mp }$ fluctuations at perpendicular scales comparable to  $\unicode[STIX]{x1D706}$ . To understand how a $\unicode[STIX]{x1D6FF}z_{\unicode[STIX]{x1D706}}^{\pm }$ wave packet cascades, it is helpful to consider a propagating ‘slice’ of the wave packet – i.e. a single cross-section of the wave packet in the plane perpendicular to the background magnetic field (see, e.g. Lithwick et al. Reference Lithwick, Goldreich and Sridhar2007). This slice ‘collides’ with a series of counter-propagating $\unicode[STIX]{x1D6FF}z_{\unicode[STIX]{x1D706}}^{\mp }$ wave packets. Each collision has a duration of

(5.2) $$\begin{eqnarray}t_{\unicode[STIX]{x1D706}}^{\pm }\sim \frac{l_{\unicode[STIX]{x1D706}}^{\mp }}{v_{\text{A}}}.\end{eqnarray}$$

The instantaneous rate at which $\unicode[STIX]{x1D6FF}z_{\unicode[STIX]{x1D706}}^{\mp }$ wave packets shear $\unicode[STIX]{x1D6FF}z_{\unicode[STIX]{x1D706}}^{\pm }$ wave packets is  ${\sim}\unicode[STIX]{x1D6FF}z_{\unicode[STIX]{x1D706}}^{\mp }/\unicode[STIX]{x1D706}$ . During a single collision, the aforementioned ‘slice’ of the $\unicode[STIX]{x1D6FF}z_{\unicode[STIX]{x1D706}}^{\pm }$ fluctuation undergoes a fractional distortion of order (Goldreich & Sridhar Reference Goldreich and Sridhar1995; Goldreich & Sridhar Reference Goldreich and Sridhar1997; Lithwick et al. Reference Lithwick, Goldreich and Sridhar2007)

(5.3) $$\begin{eqnarray}\unicode[STIX]{x1D712}_{\unicode[STIX]{x1D706}}^{\pm }=\frac{\unicode[STIX]{x1D6FF}z_{\unicode[STIX]{x1D706}}^{\mp }l_{\unicode[STIX]{x1D706}}^{\mp }}{\unicode[STIX]{x1D706}v_{\text{A}}}.\end{eqnarray}$$

5.1.1 Weak balanced turbulence

If $\unicode[STIX]{x1D712}_{\unicode[STIX]{x1D706}}^{\pm }\ll 1$ , a $\unicode[STIX]{x1D6FF}z_{\unicode[STIX]{x1D706}}^{\pm }$ wave packet undergoes only a small fractional change during each collision, and the turbulence is weak. Ng & Bhattacharjee (Reference Ng and Bhattacharjee1996, Reference Ng and Bhattacharjee1997) and Goldreich & Sridhar (Reference Goldreich and Sridhar1997) advanced a phenomenological model of weak, incompressible, MHD turbulence in which the effects of consecutive collisions are uncorrelated and add like a random walk. After $N$ collisions, the r.m.s. fractional change in a $\unicode[STIX]{x1D6FF}z_{\unicode[STIX]{x1D706}}^{\pm }$ wave packet is ${\sim}N^{1/2}\unicode[STIX]{x1D712}_{\unicode[STIX]{x1D706}}^{\pm }$ . After $N\sim (\unicode[STIX]{x1D712}_{\unicode[STIX]{x1D706}}^{\pm })^{-2}$ collisions, the r.m.s. fractional distortion of the wave packet grows to a value of order unity, and the energy contained within the wave packet cascades to smaller scales. The cascade time scale is thus

(5.4) $$\begin{eqnarray}\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D706}}^{\pm }\sim (\unicode[STIX]{x1D712}_{\unicode[STIX]{x1D706}}^{\pm })^{-2}t_{\unicode[STIX]{x1D706}}^{\pm }\sim \frac{\unicode[STIX]{x1D706}^{2}v_{\text{A}}}{(\unicode[STIX]{x1D6FF}z_{\unicode[STIX]{x1D706}}^{\mp })^{2}l_{\unicode[STIX]{x1D706}}^{\mp }}.\end{eqnarray}$$

Because neither the $\unicode[STIX]{x1D6FF}z_{\unicode[STIX]{x1D706}}^{+}$ nor $\unicode[STIX]{x1D6FF}z_{\unicode[STIX]{x1D706}}^{-}$ wave packet is altered significantly during any single collision, the leading and trailing edges of a $\unicode[STIX]{x1D6FF}z_{\unicode[STIX]{x1D706}}^{\pm }$ wave packet are sheared in virtually the same way during each collision, and the parallel length scale of the wave packets does not change as the fluctuation energy cascades to smaller  $\unicode[STIX]{x1D706}$ (Shebalin, Matthaeus & Montgomery Reference Shebalin, Matthaeus and Montgomery1983); i.e.

(5.5) $$\begin{eqnarray}l_{\unicode[STIX]{x1D706}}^{\pm }\propto \unicode[STIX]{x1D706}^{0}.\end{eqnarray}$$

In the inertial range, the $\boldsymbol{z}^{\pm }$ energy-cascade rate (per unit mass), $\unicode[STIX]{x1D716}^{\pm }$ , is independent of scale:

(5.6) $$\begin{eqnarray}\unicode[STIX]{x1D716}^{\pm }\sim \frac{(\unicode[STIX]{x1D6FF}z_{\unicode[STIX]{x1D706}}^{\pm })^{2}}{\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D706}}^{\pm }}\propto \unicode[STIX]{x1D706}^{0}.\end{eqnarray}$$

Equations (5.1), (5.4), (5.5) and (5.6) imply that

(5.7) $$\begin{eqnarray}\unicode[STIX]{x1D6FF}z_{\unicode[STIX]{x1D706}}^{\pm }\propto \unicode[STIX]{x1D706}^{1/2}.\end{eqnarray}$$

The scaling of the one-dimensional power spectrum of the $\boldsymbol{z}^{\pm }$ fluctuations, denoted $E^{\pm }(k_{\bot })$ , follows from the relation

(5.8) $$\begin{eqnarray}k_{\bot }E^{\pm }(k_{\bot })\sim (\unicode[STIX]{x1D6FF}z_{\unicode[STIX]{x1D706}}^{\pm })^{2}|_{\unicode[STIX]{x1D706}=k_{\bot }^{-1}},\end{eqnarray}$$

where $k_{\bot }$ is the component of the wave vector perpendicular to the background magnetic field. Equations (5.7) and (5.8) imply that

(5.9) $$\begin{eqnarray}E^{\pm }(k_{\bot })\propto k_{\bot }^{-2}.\end{eqnarray}$$

The scaling in (5.9) has been found in direct numerical simulations (Perez & Boldyrev Reference Perez and Boldyrev2008) as well as in exact solutions to the weak-turbulence wave kinetic equations for incompressible MHD turbulence (Galtier et al. Reference Galtier, Nazarenko, Newell and Pouquet2000). It is worth noting, however, that in weak-turbulence theory all AWs are cascaded by $k_{\Vert }=0$ modes, where $k_{\Vert }$ is the wave-vector component along  $\boldsymbol{B}_{0}$ , and these zero-frequency modes violate the assumptions of weak-turbulence theory. Several studies have addressed this issue, as well as its consequences for imbalanced turbulence (Boldyrev & Perez Reference Boldyrev and Perez2009; Schekochihin, Nazarenko & Yousef Reference Schekochihin, Nazarenko and Yousef2012; Meyrand, Kiyani & Galtier Reference Meyrand, Kiyani and Galtier2015), as discussed further in § 5.2.1.

5.1.2 Strong balanced turbulence

If $\unicode[STIX]{x1D712}_{\unicode[STIX]{x1D706}}^{\pm }\gtrsim 1$ , then each slice of a $\unicode[STIX]{x1D6FF}z_{\unicode[STIX]{x1D706}}^{\pm }$ wave packet is strongly distorted during a single collision, the turbulence is strong and $\boldsymbol{z}^{\pm }$ energy at scale  $\unicode[STIX]{x1D706}$ cascades to smaller scales on the time scale

(5.10) $$\begin{eqnarray}\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D706}}^{\pm }\sim \frac{\unicode[STIX]{x1D706}}{\unicode[STIX]{x1D6FF}z_{\unicode[STIX]{x1D706}}^{\mp }},\end{eqnarray}$$

leading to a scale-independent energy-cascade rate

(5.11) $$\begin{eqnarray}\unicode[STIX]{x1D716}^{\pm }\sim \frac{(\unicode[STIX]{x1D6FF}z_{\unicode[STIX]{x1D706}}^{\pm })^{2}\unicode[STIX]{x1D6FF}z_{\unicode[STIX]{x1D706}}^{\mp }}{\unicode[STIX]{x1D706}}\propto \unicode[STIX]{x1D706}^{0}.\end{eqnarray}$$

Equations (5.1) and (5.11) imply that

(5.12) $$\begin{eqnarray}\unicode[STIX]{x1D6FF}z_{\unicode[STIX]{x1D706}}^{\pm }\propto \unicode[STIX]{x1D706}^{1/3},\end{eqnarray}$$

which implies via (5.8) that

(5.13) $$\begin{eqnarray}E^{\pm }(k_{\bot })\propto k_{\bot }^{-5/3}.\end{eqnarray}$$

Goldreich & Sridhar (Reference Goldreich and Sridhar1995) conjectured that in strong, balanced, RMHD turbulence (and also in anisotropic, incompressible, MHD turbulence), the linear and nonlinear time scales of each wave packet are comparable, i.e.

(5.14) $$\begin{eqnarray}\unicode[STIX]{x1D712}_{\unicode[STIX]{x1D706}}^{\pm }\sim 1.\end{eqnarray}$$

Numerical simulations confirm that this ‘critical-balance’ conjecture describes strong RMHD turbulence not only on average (Cho & Vishniac Reference Cho and Vishniac2000), but structure by structure (Mallet, Schekochihin & Chandran Reference Mallet, Schekochihin and Chandran2015). Together, equations (5.12) and (5.14) imply that

(5.15) $$\begin{eqnarray}l_{\unicode[STIX]{x1D706}}^{\pm }\propto \unicode[STIX]{x1D706}^{2/3}.\end{eqnarray}$$

Several studies have argued, on the basis of numerical simulations and theoretical arguments, that the inertial-range power spectrum in strong, balanced, RMHD turbulence is flatter than in the Goldreich–Sridhar model and closer to $k_{\bot }^{-3/2}$ , because of scale-dependent dynamic alignment (Boldyrev Reference Boldyrev2005, Reference Boldyrev2006; Mason et al. Reference Mason, Cattaneo and Boldyrev2008; Perez et al. Reference Perez, Mason, Boldyrev and Cattaneo2012) and/or intermittency (Maron & Goldreich Reference Maron and Goldreich2001; Chandran et al. Reference Chandran, Schekochihin and Mallet2015b ; Mallet & Schekochihin Reference Mallet and Schekochihin2017). On the other hand, Beresnyak (Reference Beresnyak2012, Reference Beresnyak2014) argued for a scaling closer to  $k_{\bot }^{-5/3}$ based on the Reynolds-number scaling of the amplitude of dissipation-scale structures. A possible resolution of the disagreement between these two sets of studies was provided by Loureiro & Boldyrev (Reference Loureiro and Boldyrev2017a ), Loureiro & Boldyrev (Reference Loureiro and Boldyrev2017b ) and Mallet, Schekochihin & Chandran (Reference Mallet, Schekochihin and Chandran2017a ), Mallet, Schekochihin & Chandran (Reference Mallet, Schekochihin and Chandran2017b ), who investigated the disruption of sheet-like structures in RMHD turbulence by the tearing instability and magnetic reconnection (see also Pucci & Velli Reference Pucci and Velli2014; Pucci et al. Reference Pucci, Velli, Tenerani and Del Sarto2018; Vech et al. Reference Vech, Mallet, Klein and Kasper2018).

5.2 Imbalanced RMHD turbulence in homogeneous plasmas

In ‘imbalanced turbulence’, one of the Elsasser variables, say  $\boldsymbol{z}^{+}$ , has a substantially higher r.m.s. amplitude than the other,

(5.16) $$\begin{eqnarray}z_{\text{rms}}^{+}>z_{\text{rms}}^{-}.\end{eqnarray}$$

Equation (5.16) includes the highly imbalanced case, in which $z_{\text{rms}}^{+}\gg z_{\text{rms}}^{-}$ , as well as moderately imbalanced turbulence, in which, e.g. $z_{\text{rms}}^{+}\simeq 2z_{\text{rms}}^{-}$ .

5.2.1 Weak imbalanced turbulence

When (5.16) is satisfied and

(5.17a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D712}_{\unicode[STIX]{x1D706}}^{+}\ll 1\quad \unicode[STIX]{x1D712}_{\unicode[STIX]{x1D706}}^{-}\ll 1,\end{eqnarray}$$

the turbulence is both imbalanced and weak. Galtier et al. (Reference Galtier, Nazarenko, Newell and Pouquet2000) showed that in the weak-turbulence theory of imbalanced incompressible MHD turbulence,

(5.18) $$\begin{eqnarray}\unicode[STIX]{x1D6FC}^{+}+\unicode[STIX]{x1D6FC}^{-}=4,\end{eqnarray}$$

where

(5.19) $$\begin{eqnarray}E^{\pm }(k_{\bot })\propto k_{\bot }^{-\unicode[STIX]{x1D6FC}^{\pm }},\end{eqnarray}$$

the homogeneous-turbulence version of (3.35). Lithwick & Goldreich (Reference Lithwick and Goldreich2003) argued that in weak incompressible MHD turbulence, the spectra are ‘pinned’ at the dissipation wavenumber  $k_{\bot d}$ , with $E^{+}(k_{\bot d})=E^{-}(k_{\bot d})$ , and that the more energetic Elsasser variable has the steeper inertial-range power spectrum. Boldyrev & Perez (Reference Boldyrev and Perez2009) espoused a different picture, in which a ‘condensate’ of magnetic fluctuations at $k_{\Vert }=0$ dominates the energy cascade, leading to a state in which  $\unicode[STIX]{x1D6FC}^{+}=\unicode[STIX]{x1D6FC}^{-}=2$ . Schekochihin et al. (Reference Schekochihin, Nazarenko and Yousef2012) developed a theory accounting for both weakly turbulent AWs with non-zero  $k_{\Vert }$ and two-dimensional modes with $k_{\Vert }=0$ , and found that $\unicode[STIX]{x1D6FC}^{+}=\unicode[STIX]{x1D6FC}^{-}=2$ for the weakly turbulent modes and  $\unicode[STIX]{x1D6FC}^{+}=\unicode[STIX]{x1D6FC}^{-}=1$ for the two-dimensional modes in the imbalanced case.

5.2.2 Strong imbalanced turbulence

When (5.16) is satisfied and $\unicode[STIX]{x1D712}_{\unicode[STIX]{x1D706}}^{+}$ or $\unicode[STIX]{x1D712}_{\unicode[STIX]{x1D706}}^{-}$ is ${\gtrsim}1$ , the turbulence is considered strong. A number of authors have developed models of strong imbalanced MHD turbulence (e.g. Beresnyak & Lazarian Reference Beresnyak and Lazarian2008; Chandran Reference Chandran2008a ; Beresnyak & Lazarian Reference Beresnyak and Lazarian2009; Perez & Boldyrev Reference Perez and Boldyrev2009; Perez & Boldyrev Reference Perez and Boldyrev2010; Podesta & Bhattacharjee Reference Podesta and Bhattacharjee2010). Here we focus on the study by Lithwick et al. (Reference Lithwick, Goldreich and Sridhar2007) (hereafter LGS), who explored an assumption about the forcing of outer-scale $\boldsymbol{z}^{-}$ fluctuations that turns out to be particularly relevant to inhomogeneous reflection-driven MHD turbulence in the solar wind.

LGS assumed, in addition to (5.16), that

(5.20) $$\begin{eqnarray}\unicode[STIX]{x1D712}_{\unicode[STIX]{x1D706}}^{-}\gtrsim 1.\end{eqnarray}$$

Equation (5.20) implies that $\unicode[STIX]{x1D6FF}z_{\unicode[STIX]{x1D706}}^{-}$  fluctuations are sheared on a time scale  $\unicode[STIX]{x1D706}/\unicode[STIX]{x1D6FF}z_{\unicode[STIX]{x1D706}}^{+}$ that is comparable to or less than the time  $l_{\unicode[STIX]{x1D706}}^{+}/v_{\text{A}}$ for a slice of a $\unicode[STIX]{x1D6FF}z_{\unicode[STIX]{x1D706}}^{-}$ wave packet to pass through a counter-propagating $\unicode[STIX]{x1D6FF}z_{\unicode[STIX]{x1D706}}^{+}$ wave packet. The cascade time scale for $\unicode[STIX]{x1D6FF}z_{\unicode[STIX]{x1D706}}^{-}$ wave packets is therefore

(5.21) $$\begin{eqnarray}\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D706}}^{-}\sim \frac{\unicode[STIX]{x1D706}}{\unicode[STIX]{x1D6FF}z_{\unicode[STIX]{x1D706}}^{+}}.\end{eqnarray}$$

LGS argued that, since a $\unicode[STIX]{x1D6FF}z_{\unicode[STIX]{x1D706}}^{-}$ wave packet cascades after it has propagated along the background magnetic field for a distance  ${\sim}v_{\text{A}}\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D706}}^{-}$ , the parallel correlation length of the $\unicode[STIX]{x1D6FF}z_{\unicode[STIX]{x1D706}}^{-}$ wave packet is

(5.22) $$\begin{eqnarray}l_{\unicode[STIX]{x1D706}}^{-}\sim v_{\text{A}}\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D706}}^{-}\sim \frac{v_{\text{A}}\unicode[STIX]{x1D706}}{\unicode[STIX]{x1D6FF}z_{\unicode[STIX]{x1D706}}^{+}}.\end{eqnarray}$$

LGS further argued that, since $\unicode[STIX]{x1D6FF}z_{\unicode[STIX]{x1D706}}^{+}$ wave packets separated by a distance $l_{\unicode[STIX]{x1D706}}^{-}$ along the magnetic field are sheared by uncorrelated $\unicode[STIX]{x1D6FF}z_{\unicode[STIX]{x1D706}}^{-}$ fluctuations,

(5.23) $$\begin{eqnarray}l_{\unicode[STIX]{x1D706}}^{+}\sim l_{\unicode[STIX]{x1D706}}^{-}.\end{eqnarray}$$

It follows from (5.3), (5.22) and (5.23) that

(5.24a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D712}_{\unicode[STIX]{x1D706}}^{-}\sim 1\quad \unicode[STIX]{x1D712}_{\unicode[STIX]{x1D706}}^{+}\sim \frac{\unicode[STIX]{x1D6FF}z_{\unicode[STIX]{x1D706}}^{-}}{\unicode[STIX]{x1D6FF}z_{\unicode[STIX]{x1D706}}^{+}}<1.\end{eqnarray}$$

The apparent implication of the second half of (5.24a,b ), particularly when $\unicode[STIX]{x1D6FF}z_{\unicode[STIX]{x1D706}}^{-}/\unicode[STIX]{x1D6FF}z_{\unicode[STIX]{x1D706}}^{+}\ll 1$ , is that $\unicode[STIX]{x1D6FF}z_{\unicode[STIX]{x1D706}}^{+}$ wave packets cascade in a weakly turbulent manner, through multiple, uncorrelated collisions with $\unicode[STIX]{x1D6FF}z_{\unicode[STIX]{x1D706}}^{-}$ wave packets, each of which leads to a small fractional change in the $\unicode[STIX]{x1D6FF}z_{\unicode[STIX]{x1D706}}^{+}$ wave packet of order $\unicode[STIX]{x1D712}_{\unicode[STIX]{x1D706}}^{+}$ (see § 5.1.1). LGS argued, however, that each $\unicode[STIX]{x1D6FF}z_{\unicode[STIX]{x1D706}}^{+}$ wave packet is in fact sheared coherently throughout its lifetime, even when $\unicode[STIX]{x1D712}_{\unicode[STIX]{x1D706}}^{+}\sim \unicode[STIX]{x1D6FF}z_{\unicode[STIX]{x1D706}}^{-}/\unicode[STIX]{x1D6FF}z_{\unicode[STIX]{x1D706}}^{+}\ll 1$ . To establish this conclusion, LGS considered the ‘ $\boldsymbol{z}^{+}$ frame’, which moves with $\boldsymbol{z}^{+}$  fluctuations at speed  $v_{\text{A}}$ along  $\boldsymbol{B}_{0}$ relative to the background plasma. They then proposed a thought experiment in which the amplitude of $\boldsymbol{z}^{-}$ is infinitesimal, so that $\boldsymbol{z}^{-}$ has negligible effect upon  $\boldsymbol{z}^{+}$ . The $\boldsymbol{z}^{+}$ vector field is then time-independent in the $\boldsymbol{z}^{+}$ frame. If the $\boldsymbol{z}^{+}$ fluctuations are initialized with a power-law spectrum spanning the entire inertial range, and if $\boldsymbol{z}^{-}$ fluctuations are continuously injected at the outer scale with an arbitrarily long coherence time  $T$ in the $\boldsymbol{z}^{+}$ frame, then the $\boldsymbol{z}^{-}$ fluctuations will cascade to small scales and set up not just a statistical steady state, but an actual steady state in the $\boldsymbol{z}^{+}$ frame in which the $\boldsymbol{z}^{-}$ vector field is independent of time. This latter conclusion follows because $\boldsymbol{z}^{-}$ is nonlinearly distorted by $\boldsymbol{z}^{+}$ , which is constant in time in the $\boldsymbol{z}^{+}$ frame. The $\unicode[STIX]{x1D6FF}z_{\unicode[STIX]{x1D706}}^{-}$ fluctuations encountered by a $\unicode[STIX]{x1D6FF}z_{\unicode[STIX]{x1D706}}^{+}$ wave packet are therefore coherent for an arbitrarily long time, and in particular for a time much longer than the crossing time

(5.25) $$\begin{eqnarray}t_{\text{cross},\unicode[STIX]{x1D706}}^{+}\sim \frac{l_{\unicode[STIX]{x1D706}}^{-}}{v_{\text{A}}}\end{eqnarray}$$

required for a slice of the $\unicode[STIX]{x1D6FF}z_{\unicode[STIX]{x1D706}}^{+}$ wave packet to propagate through a $\unicode[STIX]{x1D6FF}z_{\unicode[STIX]{x1D706}}^{-}$ wave packet.

Building upon this thought experiment, LGS proceeded to consider the more realistic case in which $\unicode[STIX]{x1D6FF}z_{\unicode[STIX]{x1D706}}^{-}$ is finite, but still small compared to $\unicode[STIX]{x1D6FF}z_{\unicode[STIX]{x1D706}}^{+}$ at all  $\unicode[STIX]{x1D706}$ . They made a key assumption, which we call the ‘coherence assumption’, that the coherence time  $T$ (at a fixed position in the $\boldsymbol{z}^{+}$ frame) of the forcing of outer-scale $\boldsymbol{z}^{-}$ fluctuations is at least as long as the $\boldsymbol{z}^{+}$ cascade time at the outer scale, as was the case in the thought experiment above. When the coherence assumption holds, the dominant mechanism for decorrelating the $\unicode[STIX]{x1D6FF}z_{\unicode[STIX]{x1D706}}^{-}$ fluctuations encountered by a $\unicode[STIX]{x1D6FF}z_{\unicode[STIX]{x1D706}}^{+}$ wave packet is the variation of the $\boldsymbol{z}^{+}$ vector field, not the crossing of counter-propagating wave packets, and the $\unicode[STIX]{x1D6FF}z_{\unicode[STIX]{x1D706}}^{+}$ wave packet is sheared coherently throughout its lifetime. The  $\boldsymbol{z}^{+}$ cascade time scale at scale  $\unicode[STIX]{x1D706}$ then becomes

(5.26) $$\begin{eqnarray}\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D706}}^{+}\sim \frac{\unicode[STIX]{x1D706}}{\unicode[STIX]{x1D6FF}z_{\unicode[STIX]{x1D706}}^{-}},\end{eqnarray}$$

and

(5.27) $$\begin{eqnarray}\unicode[STIX]{x1D716}^{+}\sim \frac{(\unicode[STIX]{x1D6FF}z_{\unicode[STIX]{x1D706}}^{+})^{2}}{\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D706}}^{+}}\sim \frac{(\unicode[STIX]{x1D6FF}z_{\unicode[STIX]{x1D706}}^{+})^{2}\unicode[STIX]{x1D6FF}z_{\unicode[STIX]{x1D706}}^{-}}{\unicode[STIX]{x1D706}}.\end{eqnarray}$$

Because of (5.24a,b ), $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D706}}^{-}\sim \unicode[STIX]{x1D706}/\unicode[STIX]{x1D6FF}z_{\unicode[STIX]{x1D706}}^{+}$ and

(5.28) $$\begin{eqnarray}\unicode[STIX]{x1D716}^{-}\sim \frac{(\unicode[STIX]{x1D6FF}z_{\unicode[STIX]{x1D706}}^{-})^{2}}{\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D706}}^{-}}\sim \frac{(\unicode[STIX]{x1D6FF}z_{\unicode[STIX]{x1D706}}^{-})^{2}\unicode[STIX]{x1D6FF}z_{\unicode[STIX]{x1D706}}^{+}}{\unicode[STIX]{x1D706}}.\end{eqnarray}$$

Setting $\unicode[STIX]{x1D716}^{\pm }\propto \unicode[STIX]{x1D706}^{0}$ , LGS combined equations (5.27) and (5.28) to obtain

(5.29) $$\begin{eqnarray}\unicode[STIX]{x1D6FF}z_{\unicode[STIX]{x1D706}}^{+}\propto \unicode[STIX]{x1D6FF}z_{\unicode[STIX]{x1D706}}^{-}\propto \unicode[STIX]{x1D706}^{1/3},\end{eqnarray}$$

which, via (5.8), implies that

(5.30) $$\begin{eqnarray}E^{\pm }(k_{\bot })\propto k_{\bot }^{-5/3}.\end{eqnarray}$$

5.3 Anomalous coherence in reflection-driven MHD turbulence

Velli et al. (Reference Velli, Grappin and Mangeney1989) (hereafter VGM) proposed a model of reflection-driven MHD turbulence in which the Elsasser power spectra were isotropic functions of the wavenumber  $k$ , denoted  $E^{\pm }(k)$ . They divided the $\boldsymbol{z}^{\pm }$ fluctuations into ‘primary’ and ‘secondary’ components, where the primary components of $\boldsymbol{z}^{\pm }$ had the usual phase velocities of $\boldsymbol{U}\pm \boldsymbol{v}_{\text{A}}$ . The secondary components of $\boldsymbol{z}^{\pm }$ were driven modes produced by the reflection of $\boldsymbol{z}^{\mp }$ fluctuations and as a consequence had phase velocities of $\boldsymbol{U}\mp \boldsymbol{v}_{\text{A}}$ . VGM considered the super-Alfvénic solar wind at $r>r_{\text{A}}$ and took $\boldsymbol{z}^{-}$ to be dominated by secondary fluctuations. VGM estimated the r.m.s. amplitude of the secondary component of $\boldsymbol{z}^{-}$ at scale  $k^{-1}$ , which we denote $z_{k,\text{s}}^{-}$ , to be

(5.31) $$\begin{eqnarray}z_{k,\text{s}}^{-}\sim \frac{z_{k,\text{p}}^{+}}{kv_{\text{A}}\unicode[STIX]{x1D70F}_{\text{r}}},\end{eqnarray}$$

where $z_{k,\text{p}}^{+}$ is the r.m.s. amplitude of the primary component of $\boldsymbol{z}^{+}$ at scale  $1/k$ , $\unicode[STIX]{x1D70F}_{\text{r}}$ is the reflection time scale (which depends only on the radial profile of the background flow), $z_{k,p}^{+}/\unicode[STIX]{x1D70F}_{\text{r}}$ is the rate at which $z_{k,\text{s}}^{-}$ fluctuations are produced by the reflection of $z_{k,\text{p}}^{+}$ fluctuations and $1/(kv_{\text{A}})$ is the time it takes for the secondary $\boldsymbol{z}^{-}$ fluctuations at scale  $1/k$ to propagate out of the primary $\boldsymbol{z}^{+}$ fluctuations that produced them. VGM argued that the secondary $\boldsymbol{z}^{-}$ fluctuations shear the $\boldsymbol{z}^{+}$ fluctuations coherently in time, since both fluctuation types have phase velocities of $\boldsymbol{U}+\boldsymbol{v}_{\text{A}}$ . They then set the $\boldsymbol{z}^{+}$ cascade power to be

(5.32) $$\begin{eqnarray}\unicode[STIX]{x1D716}^{+}\sim kz_{k,\text{s}}^{-}(z_{k,\text{p}}^{+})^{2}\sim \frac{(z_{k,\text{p}}^{+})^{3}}{v_{\text{A}}\unicode[STIX]{x1D70F}_{\text{r}}}\end{eqnarray}$$

and took $\unicode[STIX]{x1D716}^{+}$ to be independent of  $k$ , obtaining $z_{k,\text{p}}^{+}\propto k^{0}$ . Equations (5.8) and (5.31) then yield

(5.33a,b ) $$\begin{eqnarray}E^{+}(k)\propto k^{-1}\quad E^{-}(k)\propto k^{-3}.\end{eqnarray}$$

It is useful to compare the VGM model with the LGS model discussed in § 5.2.2. In both models, the $\boldsymbol{z}^{-}$ fluctuations are anomalously coherent in the reference frame of the $\boldsymbol{z}^{+}$ fluctuations. In the LGS model, this coherence is introduced via the ‘coherence assumption’ discussed in § 5.2.2. VGM argued that this coherence arises because of the physics of AW reflection. A key difference between the models is that VGM neglected the ‘tertiary’ small-scale $\boldsymbol{z}^{-}$ fluctuations that are produced as secondary $\boldsymbol{z}^{-}$ fluctuations cascade to small scales. In the LGS model, these tertiary $\boldsymbol{z}^{-}$ fluctuations are anomalously coherent in the $\boldsymbol{z}^{+}$ reference frame and drive the Elsasser spectra towards a $k_{\bot }^{-5/3}$ scaling rather than a $k^{-1}$  scaling.

5.4 Inverse cascade in reflection-driven MHD turbulence

van Ballegooijen & Asgari-Targhi (Reference van Ballegooijen and Asgari-Targhi2017) carried out direct numerical simulations of reflection-driven MHD turbulence in the solar corona and solar wind using a methodology similar to the one we have employed. Using their simulation data, they computed the rate  $\unicode[STIX]{x1D716}^{\pm }(k_{\bot },r,t)$ at which nonlinear interactions transfer $\boldsymbol{z}^{\pm }$ energy from perpendicular wavenumbers less than  $k_{\bot }$ to perpendicular wavenumbers greater than  $k_{\bot }$ (their equation (17) divided by  $\unicode[STIX]{x1D70C}$ , with $R\rightarrow L_{\text{box}}/2$ , $f_{\pm ,k}\rightarrow \unicode[STIX]{x1D719}_{\pm ,k}$ and $a\rightarrow k_{\bot }L_{\text{box}}/2$ ),

(5.34) $$\begin{eqnarray}\unicode[STIX]{x1D716}^{\pm }(k_{\bot },r,t)=\frac{1}{[L_{\text{box}}(r)]^{2}}\mathop{\sum }_{k_{\bot l}>k_{\bot }}\mathop{\sum }_{k_{\bot j}<k_{\bot }}\mathop{\sum }_{k_{\bot i}}M_{lji}(k_{\bot i}^{2}-k_{\bot j}^{2}-k_{\bot l}^{2})\unicode[STIX]{x1D719}_{\pm ,l}\unicode[STIX]{x1D719}_{\pm ,j}\unicode[STIX]{x1D719}_{\mp ,i},\end{eqnarray}$$

where $\unicode[STIX]{x1D719}_{\pm k}$ is the Fourier transform (in $x$ and  $y$ ) of the Elsasser streamfunction  $\unicode[STIX]{x1D719}_{\pm }$ (defined such that $\boldsymbol{z}^{\pm }=\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{\pm }\times \boldsymbol{B}_{0}/B_{0}$ ), and $M_{lji}$ is a dimensionless mode-coupling coefficient that depends upon $k_{\bot l}$ , $k_{\bot j}$ and $k_{\bot i}$ , but not upon the mode amplitudes. They found that $\unicode[STIX]{x1D716}^{+}(k_{\bot },r,t)$ became negative across a broad range of  $k_{\bot }$ within a modest range of radii just larger than the radius (or radii) at which $\text{d}v_{\text{A}}/\text{d}r$ changes signs – e.g. just beyond the Alfvén-speed maximum at $r=r_{\text{m}}$ in the subset of their simulations in which the background density was smooth.

To explain their findings, they considered two locations, one just inside the $r=r_{\text{m}}$ surface at $r=r_{1}$ and one just outside the $r=r_{\text{m}}$ surface at $r=r_{2}$ , such that $|\text{d}v_{\text{A}}/\text{d}r|$ was the same at the two radii. They noted that, because $z_{\text{rms}}^{+}\gg z_{\text{rms}}^{-}$ , $\boldsymbol{z}^{-}$ fluctuations cascade much more rapidly than $\boldsymbol{z}^{+}$ fluctuations. There thus exists a range of values of $r_{2}-r_{1}$ for which the time $\unicode[STIX]{x0394}t_{12}$ required for $\boldsymbol{z}^{+}$ fluctuations to propagate from $r_{1}$ to $r_{2}$ is small compared to the outer-scale $\boldsymbol{z}^{+}$ cascade time scale but large compared to the outer-scale $\boldsymbol{z}^{-}$ cascade time scale. For values of $r_{2}-r_{1}$ in this range,

(5.35) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D719}_{+,k}(r_{2},t) & \simeq & \displaystyle \unicode[STIX]{x1D719}_{+,k}(r_{1},t-\unicode[STIX]{x0394}t_{12})\end{eqnarray}$$

(5.36) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D719}_{-,k}(r_{2},t) & \simeq & \displaystyle -\unicode[STIX]{x1D719}_{-,k}(r_{1},t-\unicode[STIX]{x0394}t_{12}),\end{eqnarray}$$

where $\unicode[STIX]{x0394}t_{12}$ is the $\boldsymbol{z}^{+}$ propagation time between $r_{1}$ and  $r_{2}$ . Equation (5.35) holds because nonlinear interactions do not have enough time to substantially alter the $\boldsymbol{z}^{+}$ fluctuations during their transit from  $r_{1}$ to $r_{2}$ . Equation (5.36) follows from the change in sign of $\text{d}v_{\text{A}}/\text{d}r$ at $r=r_{\text{m}}$ .Footnote ⁸ Because $\unicode[STIX]{x1D719}_{-,k}$ changes sign and $\unicode[STIX]{x1D719}_{+,k}$ remains almost unchanged, $\unicode[STIX]{x1D716}^{+}$ in (5.34) changes sign between $r_{1}$ and $r_{2}$ . Between the coronal base and the Alfvén-speed maximum, nonlinear interactions set up the usual direct cascade of energy from large scales to small scales, causing $\unicode[STIX]{x1D716}^{+}$ to be positive at $r_{1}$ . At $r_{2}$ , $\unicode[STIX]{x1D716}^{+}$ thus becomes negative, indicating an inverse cascade.

van Ballegooijen & Asgari-Targhi (Reference van Ballegooijen and Asgari-Targhi2017) found that as the $\boldsymbol{z}^{+}$ fluctuations propagate farther beyond $r=r_{\text{m}}$ , they gradually adjust to the new value of $\boldsymbol{z}^{-}$ , and the direct cascade of energy from large scales to small scales resumes. This transition back to a direct cascade occurs first at large  $k_{\bot }$ (at which the nonlinear time is short) and later at small  $k_{\bot }$ . In one of their simulations, there is an inverse cascade of $\boldsymbol{z}^{+}$ energy throughout the region between the Alfvén-speed maximum at $1.4R_{\odot }$ and an outer radius of  $r=2.5R_{\odot }$ . In a second simulation, there is an inverse cascade between the Alfvén-speed maximum at $1.6R_{\odot }$ and an outer radius of  $4R_{\odot }$ . Since the outer-scale $\boldsymbol{z}^{+}$ cascade time is comparable to the time required for $v_{\text{A}}$ to change by a factor of 2 in the $\boldsymbol{z}^{+}$ reference frame (Dmitruk et al. Reference Dmitruk, Matthaeus, Milano, Oughton, Zank and Mullan2002; Chandran & Hollweg Reference Chandran and Hollweg2009), the above results indicate that the inverse cascade persists (in the $\boldsymbol{z}^{+}$ reference frame) for a time comparable to the outer-scale $\boldsymbol{z}^{+}$ energy-cascade time scale.

Although $\unicode[STIX]{x1D716}^{+}$ became negative between $r\simeq r_{\text{m}}$ and $r\simeq 2r_{\text{m}}$ in the numerical simulations of van Ballegooijen & Asgari-Targhi (Reference van Ballegooijen and Asgari-Targhi2017), the energy-dissipation rate (computed from the dissipation terms added to the governing equations) decreased by only a factor of  $\simeq 2$ within the inverse-cascade region. The reason for this is that the direct-cascade region at $r<r_{\text{m}}$ had already ‘done the work’ of transporting $\boldsymbol{z}^{+}$ energy to large  $k_{\bot }$ , and the inverse cascade between $r\simeq r_{\text{m}}$ and $r\simeq 2r_{\text{m}}$ was unable to completely evacuate the high- $k_{\bot }$ part of the spectrum.