## Reflection-driven magnetohydrodynamic turbulence in the solar atmosphere and solar wind

## Benjamin D. G. Chandran 1 † and Jean C. Perez 2

1 Department of Physics and Astronomy, University of New Hampshire, Durham, New Hampshire 03824, USA

2 Department of Aerospace, Physics and Space Sciences, Florida Institute of Technology, Melbourne, Florida, 32901, USA

(Received xx; revised xx; accepted xx)

We present three-dimensional direct numerical simulations and an analytic model of reflection-driven magnetohydrodynamic (MHD) turbulence in the solar wind. Our simulations describe transverse, non-compressive MHD fluctuations within a narrow magnetic flux tube that extends from the photosphere, through the chromosphere and corona, and out to a heliocentric distance r of 21 solar radii ( R /circledot ) . We launch outward-propagating ' z + fluctuations' into the simulation domain by imposing a randomly evolving photospheric velocity field. As these fluctuations propagate away from the Sun, they undergo partial reflection, producing inward-propagating ' z -fluctuations.' Counter-propagating fluctuations subsequently interact, causing fluctuation energy to cascade to small scales and dissipate. Our analytic model incorporates dynamic alignment, allows for strongly or weakly turbulent nonlinear interactions, and divides the z + fluctuations into two populations with different characteristic radial correlation lengths. The inertial-range power spectra of z + and z -fluctuations in our simulations evolve toward a k -3 / 2 ⊥ scaling at r &gt; 10 R /circledot , where k ⊥ is the wave-vector component perpendicular to the background magnetic field. In two of our simulations, the z + power spectra are much flatter between the coronal base and r /similarequal 4 R /circledot . We argue that these spectral scalings are caused by: (1) high-pass filtering in the upper chromosphere; (2) the anomalous coherence of inertial-range z -fluctuations in a reference frame propagating outwards with the z + fluctuations; and (3) the change in the sign of the radial derivative of the Alfvén speed at r = r m /similarequal 1 7 . R /circledot , which disrupts this anomalous coherence between r = r m and r /similarequal 2 r m . At r &gt; 1 3 . R /circledot , the turbulent-heating rate in our simulations is comparable to the turbulent-heating rate in a previously developed solar-wind model that agreed with a number of observational constraints, consistent with the hypothesis that MHD turbulence accounts for much of the heating of the fast solar wind.

Key words: solar wind - Sun: corona - turbulence - waves

## 1. Introduction

One model for the origin of the solar wind relies upon Alfvén waves (AWs) with wavelengths much larger than the proton gyroradius and frequencies much smaller than the proton cyclotron frequency. In this model, photospheric motions and/or magnetic reconnection in the solar atmosphere launch AWs into the corona and solar wind, where

the AWs undergo partial non-WKB (Wentzel-Kramers-Brillouin) reflection (Velli et al. 1989; Zhou &amp; Matthaeus 1989). Subsequent interactions between counter-propagating AW packets transfer fluctuation energy from large scales to small scales. At sufficiently small scales, the fluctuation energy dissipates. Large-scale AWs also exert an outward force on the plasma. Several studies have found that this dissipation and momentum deposition can account for much of the heating and acceleration of the solar wind (e.g., Cranmer et al. 2007; Verdini et al. 2010; Chandran et al. 2011; van der Holst et al. 2014).

A number of authors have investigated different aspects of reflection-driven magnetohydrodynamic (MHD) turbulence. For example, Heinemann &amp; Olbert (1980), Velli (1993), and Hollweg &amp; Isenberg (2007) investigated the linear AW propagation problem, accounting for radial variations in the density, outflow velocity, and magnetic-field strength. Dmitruk et al. (2002), Cranmer &amp; van Ballegooijen (2005), Verdini &amp; Velli (2007), Chandran &amp; Hollweg (2009), and Zank et al. (2018) investigated the radial evolution of MHD turbulence in the solar atmosphere and solar wind accounting for reflection and nonlinear interactions. Cranmer et al. (2007), Verdini et al. (2010), Chandran et al. (2011), van der Holst et al. (2014), and Usmanov et al. (2014) incorporated reflection-driven MHD turbulence into one-dimensional (1-D) and 3-D solar-wind models. Verdini et al. (2009) and Verdini et al. (2012) carried out numerical simulations of reflection-driven MHD turbulence, in which they approximated the nonlinear terms in the governing equations using a shell model. Dmitruk &amp; Matthaeus (2003) carried out direct numerical simulations of reflection-driven MHD turbulence (i.e., without approximating the nonlinear terms) in the corona in the absence of a background flow. van Ballegooijen et al. (2011) carried out direct numerical simulations of reflection-driven MHD turbulence in the chromosphere and corona without a background flow. Perez &amp; Chandran (2013), van Ballegooijen &amp; Asgari-Targhi (2016), and van Ballegooijen &amp; Asgari-Targhi (2017) carried out direct numerical simulations of reflection-driven MHD turbulence from the low corona to the Alfvén critical point (at a heliocentric distance r of r A ∼ 10 R /circledot ) and beyond, taking into account the solar-wind outflow velocity.

In Section 3 of this paper, we present three new direct numerical simulations of reflection-driven MHD turbulence extending from the photosphere, through the chromosphere, through a coronal hole, and out to r = 21 R /circledot . These simulations go beyond previous simulations extending to r /greaterorsimilar r A by incorporating the chromosphere. This enables us to account, at least in an approximate way, for the strong turbulence that develops in the chromosphere, which launches a broad spectrum of fluctuations into the corona (van Ballegooijen et al. 2011). Our simulations also reach larger r than the simulations of Perez &amp; Chandran (2013) and contain 16 times as many grid points in the field-perpendicular plane as the simulations of van Ballegooijen &amp; Asgari-Targhi (2017).

To offer some insight into the physical processes at work in our simulations, we present an analytic model of reflection-driven MHD turbulence in Section 4. This model accounts for the generation of inward-propagating AWs by non-WKB reflection, nonlinear interactions between counter-propagating AW packets, and the development of alignment between outward-propagating and inward-propagating fluctuations. For reasons that we describe in Sections 3 and 4, we divide the outward-propagating fluctuations into two populations with different characteristic radial correlation lengths. Our model reproduces our numerical results reasonably well.

The power-law scalings of the inertial-range power spectra in our simulations vary with radius. We discuss the causes of these variations in Section 6, after reviewing several

relevant studies in Section 5. We briefly discuss other wave-launching parameter regimes in Section 7 and phase mixing in Section 8, and we present our conclusions in Section 9.

## 2. Transverse, non-compressive fluctuations in a radially stratified corona and solar wind

We focus exclusively on non-compressive fluctuations, which are observed to dominate the energy density of solar-wind turbulence (Tu &amp; Marsch 1995), and which carry an energy flux in the low corona that is sufficient to power the solar wind (De Pontieu et al. 2007). A disadvantage of our approach is that we neglect nonlinear couplings between compressive and non-compressive fluctuations (see, e.g. Cho &amp; Lazarian 2003; Chandran 2005; Luo &amp; Melrose 2006; Chandran 2008; Yoon &amp; Fang 2009; Shoda et al. 2019), which are likely important in the solar atmosphere and solar wind. For example, the plasma density varies by a factor of ∼ 6 over a distance of a few thousand km perpendicular to the background magnetic field B 0 in the low corona (Raymond et al. 2014), which suggests that phase mixing (Heyvaerts &amp; Priest 1983) is an efficient mechanism for cascading AW energy to small scales measured perpendicular to B 0 near the Sun. † We also neglect the parametric decay of AWs into slow magnetosonic waves and counter-propagating AWs (e.g., Galeev &amp; Oraevskii 1963; Sagdeev &amp; Galeev 1969; Cohen &amp; Dewar 1974; Tenerani et al. 2017), which may cause outward-propagating AWs in the fast solar wind to acquire a k -1 ‖ spectrum by the time these fluctuations reach r = 0 3 . au (Chandran 2018), where k ‖ is the wave-vector component parallel to the background magnetic field, and 1 au is the mean Earth-Sun distance. Nevertheless, the simulations that we report in Section 3 describe an important subset of the full turbulent dynamics.

Our analysis begins with the continuity, momentum, and induction equations of ideal MHD,

$$\frac { \partial \rho } { \partial t } + \nabla \cdot ( \rho v ) = 0,$$

and

$$\rho \left ( \frac { \partial v } { \partial t } + v \cdot \nabla v \right ) = - \nabla p _ { t o t } + \frac { B \cdot \nabla B } { 4 \pi } - \rho \nabla \Phi,$$

$$\frac { \partial B } { \partial t } = \bigtriangledown \times ( v \times B ), \quad \ \ \ \ ( 2. 3 ) \\ \dots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots$$

where ρ , v , and B are the mass density, velocity, and magnetic field, Φ is the gravitational potential, p tot = p + B / π 2 8 is the total pressure, and p is the plasma pressure. We set

$$v = U + \delta v \quad B = B _ { 0 } + \delta B$$

and take the background flow velocity U to be aligned with B 0 . We neglect density fluctuations, setting

$$\delta \rho & = 0. & ( 2. 5 )$$

We assume that the fluctuations are transverse and non-compressive, i.e.,

$$\delta v \cdot B _ { 0 } = 0 \quad \delta B \cdot B _ { 0 } = 0 \quad \nabla \cdot \delta v = 0, \quad \quad \quad ( 2. 6 ) \\ \quad \_ \quad \cdot \_ \quad \cdot \_ \quad \cdot \_ \cdot \_ \cdot \_ \cdot \_ \cdot$$

and we take ρ , U , and B 0 to be steady-state solutions of equations (2.1) through (2.3) (as well as the MHD energy equation). The Alfvén velocity and Elsasser variables are

† In contrast, Helios radio occultation observations show that the fractional density variations drop to 0 1 . -0 2 . at r ∈ (5 R , /circledot 20 R /circledot ) (Hollweg et al. 2010).

given by

$$\sum _ { \Lambda } v _ { \Lambda } = \frac { B _ { 0 } } { \sqrt { 4 \pi \rho } } \quad z ^ { \pm } = \delta v \mp \delta b, \quad \sum _ { \Lambda } \quad \dots$$

where δ b = δ B / √ 4 πρ . Rewriting equations (2.2) and (2.3) in terms of z ± , we obtain (Velli et al. 1989; Zhou &amp; Matthaeus 1990)

$$( \text{$\text{$even$ all}_{\infty$, $\infty$}),$} \text{$\infty$} ) \\ \frac { \partial z ^ { \pm } } { \partial t } + ( U \pm v _ { A } ) \cdot \nabla z ^ { \pm } + z ^ { \mp } \cdot \nabla ( U \mp v _ { A } ) + \frac { 1 } { 2 } \left ( z ^ { - } - z ^ { + } \right ) \left ( \nabla \cdot v _ { A } \mp \frac { 1 } { 2 } \nabla \cdot U \right ) \\ = & - \left ( z ^ { \mp } \cdot \nabla z ^ { \pm } + \frac { \nabla p _ { \text{tot} } } { \rho } \right ). \\ \text{$as in homogeneous MHD turbulence, the $\rho^{-1} \nabla p _ { \text{tot} }$ term in $(2.8)$} \text{cancels the compressive}$$

$$( 2. 8 )$$

As in homogeneous MHD turbulence, the ρ -1 ∇ p tot term in (2.8) cancels the compressive part of the z ∓ · ∇ z ± term to maintain the condition ∇· z ± = 0 .

We assume that the background magnetic field B 0 possesses a field line that is purely radial. Working, temporarily, in spherical coordinates ( r, θ, φ ) , with θ = 0 coinciding with this radial field line, we restrict our analysis to

$$\theta \ll 1.$$

We further assume that and

$$\frac { 1 } { B _ { 0 } } \frac { \partial B _ { 0 }$$

Since z ∓ · B 0 = 0 , these assumptions imply that to leading order in θ (Chandran et al. 2015 a )

ˆ

∂

0

∂r

$$z ^ { \mp } \cdot \nabla ( U \mp v _ { \L$$

b

· ∇

=

,

(2.12)

and where

$$\begin{array} {$$

We take B 0 to be directed away from the Sun, so that z + ( z -) corresponds to outwardpropagating (inward-propagating) fluctuations (when viewed in the local plasma frame), and we define vector versions of the variables introduced by Heinemann &amp; Olbert (1980),

$$g = \frac { ( 1 + \eta ^ { 1 / 2 } ) z ^ { + }$$

where

$$\eta = \rho / \rho _ { a },$$

and ρ a is the value of ρ at the Alfvén critical point, at which U = v A . Mass conservation and flux conservation imply that

$$\frac { \rho U } { B _ { 0 } } = \text{ constant},$$

which in turn implies that

$$v _ { A \phi } = U _ { \phi } = \partial U / \partial \phi = \partial v _ { \mathrm A } / \partial \phi = 0$$

$$v _ { \text{A} } & = \eta ^ { 1 / 2 } U. & ( 2. 1 8 )$$

With the use of (2.15) and (2.18), we rewrite z ± in (2.8) in terms of g and f , obtaining the nonlinear Heinemann-Olbert equations (Heinemann &amp; Olbert 1980; Chandran &amp; Hollweg 2009),

$$\frac { \partial g } { \partial t } + ( U + v _ { A } ) \frac { \partial g } { \partial r } - \left ( \frac { U + v _ { A } } { 2 v _ { A } } \right ) \frac { d v _ { A } } { d r } f = - z ^ { - } \cdot \nabla g - \left ( \frac { 1 + \eta ^ { 1 / 2 } } { \eta ^ { 1 / 4 } } \right ) \frac { \nabla p _ { t o t } } { \rho } \quad ( 2. 1 9 )$$

$$\frac { \partial f } { \partial t } + ( U - v _ { A } ) \frac { \partial f } { \partial r } - \left ( \frac { U - v _ { A } } { 2 v _ { A } } \right ) \frac { \text{d} v _ { A } } { \text{d} r } g = - z ^ { + } \cdot \nabla f - \left ( \frac { 1 - \eta ^ { 1 / 2 } } { \eta ^ { 1 / 4 } } \right ) \frac { \nabla p _ { \text{tot} } } { \rho }. \quad ( 2. 2 0 ) \\ \text{Fatalations} \, ( 2. 1 9 ) \text{ and } ( 2. 2 0 ) \text{ are equivalent to the equations solved by Perez } \& \text{ Chandran}$$

Because (2.6) is also satisfied by non-compressive fluctuations in reduced MHD (RMHD), (2.19) and (2.20) could be viewed as an inhomogeneous version of RMHD. However, the way in which we have arrived at (2.19) and (2.20) - in particular, starting with (2.5) and (2.6) as assumptions - differs from the usual derivation of the RMHD equations (see, e.g., Schekochihin et al. 2009), which begins by assuming that δB /lessmuch B 0 and λ /lessmuch l , where λ ( ) l is the characteristic length scale of the fluctuations measured perpendicular (parallel) to B 0 . We conjecture that (2.19) and (2.20) may provide a reasonable description of transverse, non-compressive fluctuations and their mutual interactions even when the assumptions δB /lessmuch B 0 and λ /lessmuch l fail. For example, if collisionless damping (Barnes 1966) or passive-scalar mixing (Schekochihin et al. 2016; Meyrand et al. 2019) removes compressive and longitudinal fluctuations, then (2.5) and (2.6) may be reasonable approximations even if δB ∼ B 0 and λ ∼ l . We note that neither our derivation of (2.19) and (2.20), nor the derivation of RMHD as a limit of the Vlasov equation (Schekochihin et al. 2009), requires that β = 8 πp/B 2 be ordered as either large or small.

Equations (2.19) and (2.20) are equivalent to the equations solved by Perez &amp; Chandran (2013) and van Ballegooijen &amp; Asgari-Targhi (2016, 2017). †

## 3. Direct Numerical Simulations

We have carried out three direct numerical simulations of (2.19) and (2.20) using the pseudo-spectral/Chebyshev REFLECT code (Perez &amp; Chandran 2013). In each simulation, the numerical domain is a narrow magnetic flux tube with a square cross section, as illustrated in Fig. 1. This flux tube extends from the photosphere at r = r min = 1 R /circledot , through the chromosphere, the 'transition region' (the narrow layer at the top of the chromosphere), and a coronal hole, and then out to a heliocentric distance of

$$r _ { \max } & = 2 1 R _ { \odot }.$$

We model the transition region in our simulations as a discontinuity in the density at

$$r _ { \text{tr} } = 1. 0 0 2 6 R _ { \odot },$$

a distance of roughly 1800 km above the photosphere. (We have collected in table 2 several heliocentric distances that we refer to repeatedly in the discussion to follow.) The walls of the simulation domain are parallel to the background magnetic field B 0 . As r

† Equations (2.1), (2.2), (2.3), (2.19), (2.20), and the plasma internal-energy equation possess two conservation laws involving f and g . The first is total-energy conservation, and the second is some times referred to as 'non-WKB wave-action conservation' (Heinemann &amp; Olbert 1980; Cranmer &amp; van Ballegooijen 2005; Verdini &amp; Velli 2007; Chandran et al. 2015 a ). This second conservation relation can be derived from the equation of cross-helicity conservation (Chandran et al. 2015 a ).

Table 1. Simulation Parameters

| Quantity                                                                                          | Run 1                                                              | Run 2                                                              | Run 3                                                              |
|---------------------------------------------------------------------------------------------------|--------------------------------------------------------------------|--------------------------------------------------------------------|--------------------------------------------------------------------|
| δv ph , rms τ (ph) v L box (1 R /circledot ) L box (1 . 0026 R /circledot ) Number of grid points | 1 . 3 km/s 3 . 3 min 4 . 1 × 10 2 km 4 . 1 × 10 3 km 256 2 × 16385 | 1 . 3 km/s 9 . 6 min 4 . 1 × 10 2 km 4 . 1 × 10 3 km 256 2 × 16385 | 1 . 3 km/s 9 . 3 min 1 . 6 × 10 3 km 1 . 6 × 10 4 km 256 2 × 16385 |

δv ph rms , is the r.m.s. amplitude of the velocity fluctuation at the photosphere, τ (ph) v is the correlation time of the photospheric velocity, and L box is the perpendicular dimension (along either the x or y directions) of the numerical domain.

Table 2. Glossary of Heliocentric Distances

| Symbol                 | Numerical Value                                                                                     | Corresponding location                                                                                   |
|------------------------|-----------------------------------------------------------------------------------------------------|----------------------------------------------------------------------------------------------------------|
| r tr r b r m r A r max | 1 . 0026 R /circledot 1 . 0027 R /circledot 1 . 71 R /circledot 11 . 1 R /circledot 21 R /circledot | transition region coronal base Alfvén-speed maximum Alfvén critical point maximum r in simulation domain |

increases and B 0 ( r ) decreases, the width L box of the simulation domain perpendicular to B 0 grows according to the relation

$$L _ { \text{box} } ( r ) = L _ { \text{box} } ( 1 R _ { \odot } ) \left [ \frac { B _ { 0 } ( 1 R _ { \odot } ) } { B _ { 0 } ( r ) } \right ] ^ { 1 / 2 }. \quad & ( 3. 3 ) \\ \text{one channel, but down to the button, but the function, too, } \text{$\infty$} \, ( 2 \, 1 \, 1 \, \cdot \, )$$

Because B 0 drops sharply between the photosphere and the transition region (see (3.16) below), L box ( r tr ) /similarequal 10 L box (1 R /circledot ) . The values of L box (1 R /circledot ) and L box ( r tr ) in our three simulations are listed in table 1. We discuss why we choose these values for L box (1 R /circledot ) in Section 3.2.

At r &gt; r tr , the field lines of B 0 are nearly radial, even though we allow for super-radial expansion of the magnetic field. This is because the flux-tube width is much smaller than the characteristic radial distance over which B 0 varies by a factor of order unity. Because the flux tube is narrow and B 0 is nearly radial, we can ignore the curvature of the field-perpendicular surfaces to a good approximation at r &gt; r tr . We thus use Cartesian coordinates, x and y , to denote position in the plane perpendicular to the radial line that runs down the centre of the simulation domain.

At r &lt; r tr , our assumption in Section 2 that B 0 is nearly radial breaks down, because the flux tube expands so rapidly with height above the photosphere. Because of this, and because we neglect compressive fluctuations, our simulations provide only a crude approximation of chromospheric turbulence. Nevertheless, we retain the chromosphere in our simulations, because turbulence in the actual chromosphere launches a broad spec-

Figure 1. Numerical domain of the REFLECT Code.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00880v2-figure-4cf28d5a-bed9-4296-b75e-5d38b300bad8.jpg)

trum of AWs into the corona (van Ballegooijen et al. 2011), and our model chromosphere gives us a way of approximating this turbulent wave-launching process.

$$3. 1. \ R a d i a l \ P r o f i l e s \ o f \, \rho, \, B _ { 0 }, \, a n d \ U$$

We choose the radial profiles of ρ , U , and B 0 to approximate the conditions found in coronal holes and the fast solar wind. Above the transition region, at r &gt; r tr , we set

$$\rho = \left ( 1 0 ^ { 9 } s ^ { - 1 5. 6 } + 2. 5 1 \times 1 0 ^ { 6 } s ^ { - 3. 7 6 } + 1. 8 5 \times 1 0 ^ { 5 } s ^ { - 2 } \right ) m _ { p } \text{ cm} ^ { - 3 }, \quad \quad ( 3. 4 ) \\ B _ { 0 } = 1. 5 \left [ s ^ { - 6 } ( f _ { \max } - 1 ) + s ^ { - 2 } \right ] \text{ G},$$

and

$$B _ { 0 } = 1. 5 \left [ s ^ { - 6 } ( f _ { \max } - 1 ) + s ^ { - 2 } \right ] G, \quad \quad \quad$$

where

$$U = 9. 2 5 \times 1 0 ^ { 1 2 } \left ( \frac { B _ { 0 } } { 1 \, G } \right ) \left ( \frac { \rho } { m$$

$$s = \frac { r } { R _ { \odot } },$$

$$f _ { \max } = 9$$

is the super-radial expansion factor, and m p is the proton mass. Equation (3.4) is adapted from the coronal-hole electron-density measurements of Feldman et al. (1997). 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 &amp; van Ballegooijen (2005). Equation (3.5) is taken from Hollweg &amp; Isenberg (2002). 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

$$U ( r _ { b } ) = 1. 2 \, \text{km} / s \quad U ( 1 \, \text{au} ) = 7 5 0 \, \text{km} / s,$$

where

$$r _ { b } & = 1. 0 0 2 7 R _ { \odot } & ( 3. 1 0 )$$

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

$$r _ { \mathrm A } & = 1 1. 1 R _ { \odot }, & ( 3. 1 1 )$$

the Alfvén speed reaches its maximum value at

$$r _ { m } & = 1. 7 1 R _ { \odot }, & ( 3. 1 2 )$$

and

$$v _ { A } ( r _ { b } ) = 9 3 5 \, \text{km} / s \quad v _ { A } ( r _ { m } ) = 2 7 3 0 \, \text{km} / s \quad v _ { A } ( r _ { A } ) = U ( r _ { A } ) = 6 2 7 \, \text{km} / s. \ \ ( 3. 1 3 )$$

Below the transition region, we set

$$\rho = \rho _ { \text{ph} } e ^ { c ( 1 - s ) / s },$$

where

$$\rho _ { \text{ph} } = 4. 7 8 \times 1 0 ^ { 1 6 } m _ { \text{p} } \, \text{cm} ^ { - 3 } \quad \text{(3.15)} \\ \text{longit} \nu \, \varepsilon = \int _ { \ell } ( \nu ) \, \int _ { \ell } ( \nu ) \, \cdot \, \int _ { \ell } \, - \, \nu \, \int _ { \ell } \, \text{and} \, \alpha \, \cdot \, \text{is the} }$$

is the photospheric density, c = [ s tr / (1 -s tr )] ln( ρ tr ,&lt; /ρ ph ) , s tr = r tr /R /circledot , and ρ tr ,&lt; 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 tr from (3.4). We then set (cf. van Ballegooijen et al. 2011)

$$B = \left [ \frac { ( B _ { \text{ph} } ^ { 2 } - B _ { \text{tr} } ^ { 2 } ) ( \rho - \rho _ { \text{tr}, < } ) } { \rho _ { \text{ph} } - \rho _ { \text{tr}, < } } + B _ { \text{tr} } ^ { 2 } \right ] ^ { 1 / 2 },$$

at r &lt; r tr , where

$$B _ { \text{ph} } = 1 4 0 0 \, G$$

is the assumed magnetic-field strength in the photospheric footpoint of the simulated flux tube, and B tr is the value of B at r = r tr from (3.5).

We plot the radial profiles of ρ , B U , , and v A in figure 2. We also plot the z + travel time between the photosphere and radius r ,

$$\cdot \\ T ( r ) = \int _ { R _ { \odot } } ^ { r } \frac { \mathrm d r } { U + v _ { \mathrm A } }.$$

## 3.2. Boundary Conditions

We take the z ± 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 ⊥ &gt; 3 × 2 π/L box ( R /circledot ) , where

$$k _ { \perp } = \sqrt { k _ { x } ^ { 2 } + k _ { y } ^ { 2 } }, & & ( 3. 1 9 ) \\ \text{ad $y$ components of the wave vector $k$. We set the amplitudes}$$

and k x and k y are the x and y components of the wave vector k . We set the amplitudes of the velocity Fourier components at k ⊥ /lessorequalslant 3 × 2 π/L box ( R /circledot ) equal to a constant, which we choose so that the root-mean-square (r.m.s.) amplitude of the fluctuating velocity at the photosphere is

$$\delta v _ { \text{ph,rms} } = 1. 3 \, \text{km/s},$$

consistent with observational constraints on the velocities of solar granules (Richardson &amp; Schwarzsc 1950). We then assign random values to the phases of these velocity Fourier components at the discrete set of times t n = nτ , 0 where τ 0 = 5 min in Run 1 and τ 0 = 20 min in Runs 2 and 3. To determine the phases at times between successive t n , we use cubic interpolation in time. We define the correlation time of the photospheric velocity τ (ph) v 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 τ 0 and L box ( R /circledot ) determine (at least in part - see Section 3.7) the correlation time τ c and perpendicular correlation length L ⊥ of the AWs launched by the Sun. (Since we only drive photospheric velocity modes with k ⊥ /lessorequalslant 3 × 2 π/L box ( R /circledot ) , L ⊥ is a few times smaller than L box .) Estimates of L ⊥ ( r b ) range from /similarequal 10 3 km (Cranmer et al. 2007; Hollweg et al. 2010; van Ballegooijen &amp; Asgari-Targhi 2016, 2017) to more than 10 4 km (Dmitruk et al. 2002; Verdini &amp; Velli 2007; Verdini et al. 2012), and estimates of τ c ( r b ) range from /similarequal 1 -5 minutes (Cranmer &amp; van Ballegooijen 2005; van Ballegooijen &amp; Asgari-Targhi 2016, 2017) to one or more hours (Dmitruk &amp; Matthaeus

Figure 2. The radial profiles of the solar-wind outflow velocity U , Alfvén speed v A , plasma density ρ divided by the proton mass m p , background magnetic-field strength B 0 , and 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 Section 4.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00880v2-figure-960ec607-134c-4b1e-a561-f109fefdc76a.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00880v2-figure-5309fb6e-f27c-4263-84a4-764d1d689a5c.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00880v2-figure-58c3e4df-6489-43de-8fed-bf20508ea573.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00880v2-figure-c4e9da0c-89b6-489f-b5ce-2b5ee0a074e8.jpg)

2003). Given the uncertainty in L ⊥ ( r b ) and τ c ( r b ) , we vary L box ( R /circledot ) and τ 0 by factors of 4 and 5, respectively, in our different simulations in order to investigate how the values of L ⊥ ( r b ) and τ c ( r b ) influence the properties of the turbulence at larger r .

No information flows into the simulation domain through the outer boundary at r = r max , because r max &gt; r 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 hyperdissipation term of the form

$$D _ { g } = - \nu _ { g } \left ( \frac { \partial ^ { 2 } } { \partial x ^ { 2 } } + \frac { \partial ^ { 2 } } { \partial y ^ { 2 } } \right ) ^ { 4 } g \quad \quad \quad$$

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

$$D _ { f } & = - \nu _ { f } \left ( \frac { \partial ^ { 2 } } { \partial x ^ { 2 } } + \frac { \partial ^ { 2 } } { \partial y ^ { 2 } } \right ) ^ { 4 } f \\ \vdots \lambda \$$

to the right-hand side of (2.20). We choose the magnitude and radial dependence of the hyper-dissipation coefficients ν g and ν f so that dissipation becomes important near the grid scale at all radii in each simulation. In particular, we take ν g and ν f to be proportional to [ L box ( r /L ) box ( R /circledot )] 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. 1988). 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 16385. Except at r min and r max , these boundary grid points are shared by neighboring subdomains. Eight of the subdomains are in the chromosphere.

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 δ v and δ B ) at the boundaries between neighboring 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 . hr. Between t = 0 and t = 4 hr, the magnetic and kinetic energies in the simulations fluctuate while trending upwards. For reference, it takes 1.3 hours for an outward-propagating AW to travel from the photosphere to the Alfvén critical point at r A = 11 1 . R /circledot , and 3 hours for an outwardpropagating AW to travel from the photosphere to r max = 21 R /circledot (see figure 2). After t /similarequal 4 hr, the magnetic and kinetic energies fluctuate around a steady value. We regard the turbulence as being in a statistical steady state at t &gt; 6 hr. All the numerical results that we present are calculated from time averages between t = 6 hr and t = 13 2 . hr, except for the z + HF rms , and z + LF rms , profiles in Run 2; those profiles, because of technical difficulties, were only computed from averages between t = 12 hr and t = 13 hr.

## 3.6. Radial profiles of the fluctuation amplitudes

In figure 3, we plot the r.m.s. amplitudes of z ± , denoted z ± rms , as a function of r in Runs 1 through 3 and in the analytic model discussed in Section 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 + rms /similarequal z -rms in the chromosphere, because of strong AW reflection at the transition region and photosphere. On the other hand, z + rms /greatermuch z -rms in the corona and solar wind because of the limited efficiency of reflection in these regions and because z -fluctuations are rapidly cascaded to small scales by the large-amplitude z + fluctuations.

The value of z + rms increases between r = R /circledot and r = 5 R /circledot because of the radially decreasing density profile. Equation (2.19) implies that the r.m.s. amplitude of g ( g rms ) is independent of r when (i) the fluctuations are in a statistical steady state, (ii) z -rms /lessmuch z + rms , and (iii) nonlinear interactions can be ignored. At r &lt; 5 R /circledot , ρ r ( ) /greatermuch ρ r ( A ) , and it follows from (2.15) that z + rms ∝ g rms ρ -1 / 4 . Equations (2.15) and (2.19) thus imply that the linear physics of AW propagation causes z + rms to increase rapidly with increasing r at r &lt; 5 R /circledot , since z -rms /lessmuch z + rms in this region. When nonlinear interactions are taken into account, g rms becomes a decreasing function of r , but the linear physics 'wins out' at r &lt; 5 R /circledot , in the sense that z + rms ∝ g rms ρ -1 / 4 remains an increasing function of r . Since the rate of non-WKB reflection vanishes at r = r m = 1 71 . R /circledot , the z -fluctuations

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00880v2-figure-e399111a-9a41-4187-8021-ae405f539465.jpg)

Figure 3. The top panels and bottom-left panel show the r.m.s. amplitudes of z ± in Runs 1 through 3 and in the analytic model described in Section 4. The lower-right panel shows δB rms /B 0 in Runs 1 through 3, where δB rms is the r.m.s. amplitude of the magnetic-field fluctuation.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00880v2-figure-a53a00e4-50f4-4d97-8180-fd3db8143e16.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00880v2-figure-72b06f23-644c-4f7b-b196-69e27742e416.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00880v2-figure-5620b69c-27ac-47bc-8be3-5c07441c7ad6.jpg)

seen at r = r m in all three simulations must be generated elsewhere. At r &lt; r A , z -fluctuations propagate with a negative radial velocity once they are produced, and thus the z -fluctuations seen at r = r m in the simulations originate at r &gt; r 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

$$L _ { \perp } & = \frac { 1 } { 3 } L _ { \text{box} }. & \quad \text{(3.23)}$$

The first time scale is the correlation time of the photospheric velocity field, τ (ph) v , 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

$$\tau _ { \text{nl} } & = \frac { L _ { \perp } } { z _ { \text{rms} } ^ { \pm } }$$

of the balanced turbulence ('balanced' meaning that z + rms /similarequal z -rms ) just below the transition region at r = r tr ,&lt; = r tr -/epsilon1 , where /epsilon1 is an infinitesimal distance, and z ± rms ( r tr ,&lt; ) /similarequal 30 km/s. (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

Figure 4. R.m.s. amplitudes of z + HF and z + LF (defined in (3.26) through (3.28) and (3.32)) in Runs 1 through 3 and in the analytic model described in Section 4.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00880v2-figure-b3caaea7-c955-4899-b2cc-cb46265174f2.jpg)

of (3.24) contains a ± sign, we do not include a ± sign on the left-hand side, because we will only evaluate (3.24) at locations at which z + rms /similarequal z -rms . We define

$$\tau _ { n l } ^ { ( t r ) } = \tau _ { n l } ( r _ { t r, < } ).$$

Given the values of L box ( r tr ) listed in table 1, τ (tr) nl is 0.8 min, 0.8 min, and 3 min in Runs 1, 2, and 3, respectively, values that are several times smaller than τ (ph) v . This suggests that the transition region in our simulations launches two populations of z + fluctuations characterized by different time scales and hence different radial correlation lengths.

To investigate this possibility, we define

$$g _ { \text{LF} } ( \tilde { x }, \tilde { y }, r, t ) = \frac { 1 } { 2 \Delta } \int _ { r _ { i } } ^ { r _ { i } + 2 \Delta } \text{d} r ^ { \prime } g \left ( \tilde { x }, \tilde { y }, r ^ { \prime }, t \right ),$$

and

$$g _ { \text{HF,rms} } = \sqrt { g _ { \text{rms} } ^ { 2 } - g _ { \text{LF,rms} } ^ { 2 } }, & & ( 3. 2 8 ) \\ \tilde { y } \, = \, y / L _ { \text{box} }, & \, \text{and} \, \langle \dots \rangle \, \text{ denotes an average over $x$, $y$, and $t$. The}$$

$$g _ { \text{LF,rms} } & = \left \langle | g _ { \text{LF} } | ^ { 2 } \right \rangle ^ { 1 / 2 }, \quad & & ( 3. 2 7 ) \\ & \sqrt { \circ } \quad & \circ \circ \circ$$

where ˜ = x x/L box , ˜ = y y/L box , and 〈 . . . 〉 denotes an average over x , y , and t . The quantity

$$\Delta = c _ { \tt a v } \tau _ { \tt n l } ^ { ( \tt t r ) } v _ { \tt A } ( r _ { \tt b } )$$

is the approximate radial correlation length in the low corona of a z + fluctuation that is

generated by a disturbance at the transition region whose correlation time is τ (tr) nl ,

$$r _ { i } = \begin{cases} \ r _ { \min } & \text{ if } r < r _ { \min } + \Delta \\ r - \Delta & \text{ if } r _ { \min } + \Delta \leqslant r \leqslant r _ { \max } - \Delta \\ r _ { \max } - 2 \Delta & \text{ if } r > r _ { \max } - \Delta \end{cases} \quad \text{ (3.30)}$$

 and c av is a dimensionless constant of order unity. We set

$$c _ { \text{av} } \simeq 0. 6, \quad \ \ & \quad ( 3. 3 1 ) \\.. \quad.. \quad.. \quad.. \quad.. \quad.. \quad.. \quad.. \quad.$$

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, ∆ = 0 08 . R /circledot in Runs 1 and 2, and ∆ = 0 32 . R /circledot in Run 3. We define

$$z _ { \text{LF,rms} } ^ { + } = \frac { \eta ^ { 1 / 4 } g _ { \text{LF,rms} } } { 1 + \eta ^ { 1 / 2 } } \quad z _ { \text{HF,rms} } ^ { + } = \frac { \eta ^ { 1 / 4 } g _ { \text{HF,rms} } } { 1 + \eta ^ { 1 / 2 } }.$$

We emphasize that, although we use the subscripts 'LF' and 'HF' as shorthand for 'lowfrequency' and 'high-frequency,' the defining difference between z + LF rms , and z + HF rms , is the difference in their radial correlation lengths.

In figure 4 we plot the radial profiles of z + LF rms , and z + HF rms , in our numerical simulations and the analytic model of Section 4. As this figure shows, all three simulations contain both z + LF and z + HF fluctuations, and these fluctuations evolve in different ways as they propagate away from the Sun. In all three runs, z + HF rms , /similarequal z + LF rms , in the low corona. As r increases, z + HF rms , /z + LF rms , decreases, particularly in Run 2, suggesting that the highfrequency component of 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 z + and z -,

$$\sin \theta & = \frac { \langle | z ^ { + } \times z ^ { - } | \rangle } { \langle | z ^ { + } | \rangle \langle | z ^ { - } | \rangle }, & ( 3. 3 3 ) \\ \text{mulations and the model we resent in Section 4. As $r$ increases.}$$

in both our numerical simulations and the model we present in Section 4. As r increases, sin θ decreases, particularly in Run 2, causing nonlinear interactions between z + and z -to weaken (see, e.g., Boldyrev 2005, 2006; Perez &amp; Chandran 2013; Chandran et al. 2015 b ). †

## 3.9. Turbulent heating

In figure 6 we plot the rate Q at which energy is dissipated per unit mass by hyperdissipation in our simulations (see Perez &amp; Chandran 2013) as a function of r , as well as the turbulent heating rate in the analytic model described in Section 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 &amp; van Ballegooijen (2005), 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. 1998; Esser et al. 1999). For comparison, the r.m.s. amplitudes of the fluctuating velocity δv rms at r = r tr in Runs 1, 2, and 3 are, respectively, 30 4 . km/s, 30 0 . km/s, and 26 7 . km/s. ‡ The values of δv rms at r =

- † We note that a different alignment angle, between δ v and δ b fluctuations, was the basis for Boldyrev's (2006) theory of scale-dependent dynamic alignment.
- ‡ In contrast to z ± , δ v is continuous across the transition region, and it makes no difference whether we evaluate δv rms at, just above, or just below the transition region.

Figure 5. The characteristic value of the sine of the alignment angle θ between z + and z -, defined in (3.33), in Runs 1 through 3 and in the analytic model of Section 4 (using (4.8)).

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00880v2-figure-7c3aa7a7-3716-4a62-8fd7-5374aa94d324.jpg)

2 R /circledot in Runs 1, 2, and 3 are, respectively, 170 km/s, 157 km/s, and 146 km/s. 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 b in the simulation were realistic. ¶

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. (2011). 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. (2011) model, which we denote Q 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 C11 at r /greaterorsimilar 4 R /circledot . However, in all three runs, Q &gt; Q C11 at r = 2 -3 R /circledot . This latter discrepancy is largest in the case of Run 2, in which Q /similarequal 3 Q C11 at r = 2 R /circledot .

¶ An important caveat to this statement is that we have neglected the interaction between non-compressive fluctuations and compressive fluctuations, including phase mixing, which we discuss in Section 8.

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

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00880v2-figure-2ffb9fe5-e870-4980-b349-1ee2953d122f.jpg)

Although Run 2 has the largest heating rate of all three simulations at r = 2 R /circledot , the simulated heating rate in Run 2 is smaller than Q C11 at r /greaterorsimilar 5 R /circledot by a factor of ∼ 2 .

The only region in which the simulated heating rate differs from Q C11 by a factor /greaterorsimilar 4 is at r &lt; 1 3 . R /circledot in Run 3, where Q/Q C11 falls below 0.1. Even in Runs 1 and 2, the simulated heating rate at r &lt; 1 3 . R /circledot is smaller than Q C11 by a factor of ∼ 2 . The finding that Q /lessorsimilar 0 5 . Q C11 at r &lt; 1 3 . R /circledot 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. (2007) and Verdini et al. (2010).

Recently, van Ballegooijen &amp; Asgari-Targhi (2016, 2017) 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. (2011) to estimate the amount of heating required, whereas van Ballegooijen &amp; Asgari-Targhi (2016, 2017) used a one-fluid solar-wind model (A. van Ballegooijen, private communication). In the Chandran et al. (2011) 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.

Figure 7. Top-left panel: the Elsasser power spectra E ± ( k ⊥ , r ) defined in (3.34) as functions of perpendicular wavenumber k ⊥ at r = 20 R /circledot in Run 1. Top-right and bottom panels: the spectral indices α + ( r ) and α -( r ) defined in (3.35) in our three numerical simulations.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00880v2-figure-ae525850-4c24-4fac-b181-823304cf049e.jpg)

## 3.10. Simulation results: Elsasser power spectra

We define the perpendicular Elsasser power spectra

$$E ^ { \pm } ( k _ { \perp }, r ) = k _ { \perp } \int _ { 0 } ^ { 2 \pi } \text{d} \phi & \overline { \left | \tilde { z } ^ { \pm } ( k _ { \perp }, \phi, r, t ) \right | ^ { 2 } }, \\ \wedge \colon = \text{$k_{k}$} \, \Xi \, \text{$\Xi$} \, \text{$\Xi$} \, \text{$\Xi$} \, \text{$\Xi$} \, \text{$\Xi$} \, \text{$\Xi$} \, \text{$\Xi$} \, \text{$\Xi$} \, \text{$\Xi$} \, \text{$\Xi$} \, \text{$\Xi$} \, \text{$\Xi$} \, \text{$\Xi$} \, \text{$\Xi$} \, \text{$\Xi$} \, \text{$\Xi$} \, \text{$\Xi$} \, \text{$\Xi$} \, \text{$\Xi$} \, \text{$\Xi$} \, \text{$\Xi$}}$$

where ˜ z ± ( k ⊥ , φ ) is the Fourier transform of z ± in x and y (see figure 1), φ is the polar angle in the ( k , k x y ) plane, and ( . . . ) indicates a time average. As illustrated in the topleft panel of figure 7, we find that E ± ( k ⊥ ) exhibits an approximate power-law scaling of the form

$$E ^ { \pm } ( k _ { \perp }, r ) \in k _ { \perp } ^ { - \alpha ^ { \pm } ( r ) } & & ( 3. 3 5 ) \\ \intertext { t o l. } \sim 1 5 \vee 2 \pi / L. \quad \text{at all $r$ in all three of our simultaneous $\mathfrak{We}$}$$

from k ⊥ /similarequal 3 × 2 π/L box to k ⊥ /similarequal 15 × 2 π/L box at all r in all three of our simulations. We evaluate α ± ( r ) by fitting E ± ( k ⊥ , r ) to a power law within this range of wave numbers, and plot the resulting values of α ± ( r ) in figure 7.

Although we drive only large-scale ( k ⊥ /lessorequalslant 3 × 2 π/L 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 z + fluctuations at the transition region and the strong reflection of z -fluctuations (at all k ⊥ ) at the photosphere (enforced by the fixed-velocity boundary condition at r = R /circledot ), which together lead to comparatively 'balanced' turbulence (meaning that z + rms /similarequal z -rms ) in the chromosphere,

as shown previously by van Ballegooijen et al. (2011). In the low chromosphere, α ± /similarequal 1 3 . -1 5 . in all four simulations, which is similar to the value α ± /similarequal 3 / 2 that arises in numerical simulations of homogeneous, balanced, RMHD turbulence (Mason et al. 2008; Perez et al. 2012; Beresnyak 2012). On the other hand, α + decreases from /similarequal 1 5 . to /similarequal 0 8 . as r increases from 1 001 . R /circledot to r tr = 1 0026 . R /circledot 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 lowerk ⊥ (and hence lower-frequency - see Goldreich &amp; Sridhar (1995)) z + fluctuations to undergo nonWKB reflection, and allowing higherk ⊥ (and hence higher-frequency) z + fluctuations to propagate unhindered to the transition region (Velli 1993; Réville et al. 2018). The difference in Run 3 is that L box is larger, and thus z + fluctuations do not reach sufficiently large k ⊥ values that they can avoid non-WKB reflection in the upper chromosphere. † The idea that z + fluctuations at the highk ⊥ 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 + rms ( r b ) and δB rms are somewhat larger in Runs 1 and 2 than in Run 3 (see table 1 and figure 3).

As r increases from r tr to r A and beyond, α ± approaches /similarequal 3 / 2 in all three runs. In Runs 1 and 2, the increase in α + as r increases from r tr to r A is not steady. In Run 1, α + decreases as r increases from 2 8 . R /circledot to 4 2 . R /circledot , and in Run 2, α + plateaus around a value of 1 between r = 2 R /circledot and r = 3 R /circledot . This behavior suggests that, in these two simulations, the turbulent dynamics at 2 R /circledot /lessorsimilar r /lessorsimilar 4 R /circledot drives α + towards a value close to 1, and the tendency for α + to evolve towards 3/2 sets in at r /greaterorsimilar 4 R /circledot . We discuss these trends further in Section 6.

## 4. Two-component analytic model of reflection-driven MHD turbulence

Chandran &amp; Hollweg (2009) (hereafter CH09) developed an analytic model of reflection-driven MHD turbulence in the solar corona and solar wind. This model can reproduce the radial profile of z + rms in our numerical simulations fairly accurately, provided the constant χ introduced in Equation (34) of CH09 is treated as an adjustable free parameter that is allowed to take on different values in different simulations. With the best-fit values of χ , the CH09 model also reproduces the turbulent-heating profiles in Runs 1 and 2 reasonably well. However, the model is significantly less accurate at reproducing Q r ( ) in Run 3 and deviates markedly from the z -rms profiles in all three runs. Moreover, the CH09 model does not explain the differences between the best-fit values of χ for Runs 1, 2, and 3 (which are, respectively, 0.55, 0.72, and 0.36), or explain how these values can be determined from the perpendicular correlation length and correlation time of the AWs launched by the Sun. These shortcomings indicate that there are important physical processes operating in our numerical simulations that were not accounted for by CH09.

In order to elucidate these processes, we develop a new analytic model of reflectiondriven MHD turbulence at

$$r \geqslant r _ { b },$$

where r b is the radius of the coronal base defined in (3.10). The reader who is not

† The transition region in our simulations is a density discontinuity, which reflects z + fluctuations with an efficiency that is independent of k ⊥ (van Ballegooijen et al. 2011); reflection at the transition region helps explain why α -is comparatively small in the upper chromosphere in Runs 1 and 2, but it does not act like a high-pass filter in our simulations.

interested in the technical details may wish to skip to Section 4.6, which summarizes the free parameters and boundary conditions of the model and compares the model with our simulation results.

We begin by dividing the z + fluctuations into two components as described in Section 3.7:

$$z ^ { + } = z _ { \text{HF} } ^ { + } + z _ { \text{LF} } ^ { + } \quad g = g _ { \text{HF} } + g _ { \text{LF} }.$$

The quantities z + HF and z + LF have different radial correlation lengths (see (3.26)), but we take them to have the same perpendicular outer scale L ⊥ . † We make the simplifying approximation that the HF and LF fluctuations are uncorrelated; i.e., g HF · g LF = 0 , where ( . . . ) indicates a time average. Non-WKB reflection is more efficient for lowfrequency AWs than for high-frequency AWs (Heinemann &amp; Olbert 1980; Velli 1993). We thus take f to be a 'low-frequency quantity' that is correlated with g LF but not g HF . Upon taking the dot product of (2.19) with 2 g HF , averaging, and assuming a statistical steady state, we obtain

$$( U + v _ { \text{A} } ) \frac { \text{d} } { \text{d} r } g _ { \text{HF}, \text{rms} } ^ { 2 } = 2 \overline { R \cdot g _ { \text{HF} } },$$

where g HF rms , = | g HF 2 1 | / 2 (with analogous definitions for g LF rms , , z + HF rms , , z + LF rms , , f rms , and z -rms ), and R represents the right-hand side of (2.19). The nonlinear term on the right-hand side of (2.19) acts to cascade g HF fluctuations to small scales at which the fluctuations dissipate. We set R g · HF = -γ + HF HF rms g 2 , , where γ + HF is the cascade rate of the outer-scale g HF fluctuations. Equation (4.3) then becomes

$$( U + v _ { \mathrm A } ) \frac { d } { d r } g _ { \text{HF}, rms } ^ { 2 } = - 2 \gamma _ { \text{HF} } ^ { + } g _ { \text{HF}, rms } ^ { 2 }.$$

We follow Velli et al. (1989) and Verdini et al. (2009) in taking the outer-scale z -fluctuations to be anomalously coherent in a reference frame that propagates outward with the z + fluctuations, because the z -fluctuations are produced by sources that propagate outward at speed U + v A . We thus estimate γ + HF using a strong-turbulence scaling regardless of the value of z -HF rms , , setting

$$\gamma _ { \text{HF} } ^ { + } = \frac { c _ { \text{diss} } z _ { \text{rms} } ^ { - } } { L _ { \perp } },$$

where c diss is a dimensionless free parameter.

Using a similar procedure, but this time for g LF , we find that

$$( U + v _ { A } ) \frac { d } { d r } g _ { L F, r m s } ^ { 2 } = - 2 \gamma _ { L F } ^ { + } g _ { L F, r m s } ^ { 2 },$$

where γ + LF is the rate at which outer-scale g LF fluctuations cascade to small scales and dissipate. In writing (4.6), we have dropped a term containing f · g LF on the assumption that f /lessmuch g LF . We set

$$\gamma _ { \text{LF} } ^ { + } = \frac { c _ { \text{diss} } z _ { \text{rms} } ^ { - } A } { L _ { \perp } }, \quad \begin{array} { c } ( 4. 7 ) \\ \end{array} \quad \begin{array} { c } ( 4. 7 ) \\ \end{array} \quad \begin{array} { c } ( 4. 7 ) \\ \end{array} \quad \begin{array} { c } ( 4. 7 ) \\ \end{array} \quad \begin{array} { c } ( 4. 7 ) \\ \end{array} \quad \begin{array} { c } ( 4. 7 ) \\ \end{array}$$

where the dimensionless coefficient A models the weakening of nonlinear interactions between g LF and f as these two fluctuation types become increasingly aligned with each

† This is an over-simplification for Runs 1 and 2, because α + ( r b ) /similarequal 0 8 . in these runs, indicating that much of the z + HF energy is concentrated at high k ⊥ . We neglect this spectral flattening in our analytic model, however, because there is minimal flattening in Run 3 and because we wish to keep the model as simple as possible.

other. We discuss how we determine A in Section 4.3 below. In order to compare our model with our simulation results, we take A to be related to sin θ in (3.33) via the equation

$$\sin \theta = \frac { 0. 5 5 ( A g _ { L F, r m s } ^ { 2 } + g _ { H F, r m s } ^ { 2 } ) } { g _ { L F, r m s } ^ { 2 } + g _ { H F, r m s } ^ { 2 } }, \quad \begin{array} { c } & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & \\. \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \end{array}$$

which expresses the idea that only low-frequency g fluctuations align with f fluctuations, while both low-frequency and high-frequency g fluctuations contribute to the average that is used to compute sin θ in (3.33). The factor of 0.55 in (4.8) is included because this is the typical value of the right-hand side of (3.33) for outer-scale fluctuations in homogeneous RMHD turbulence (Chandran et al. 2015 ). b

## 4.1. Amplitude of the inward-propagating fluctuations

To determine z -rms , we assume that z -cascades primarily via interactions with z + LF (or, equivalently, g LF ). The outer-scale z -cascade rate then depends upon the criticalbalance parameter (Goldreich &amp; Sridhar 1995; Boldyrev 2006)

$$\chi _ { L F } ^ { - } = \frac { z _ { L F, r m s } ^ { + } L _ { r, L F } ^ { + } A } { L _ { \perp } v _ { A } },$$

where L + r, LF is the radial correlation length of the g LF fluctuations. The critical-balance parameter χ -LF is an estimate of the fractional change in an outer-scale z -fluctuation that results from a single 'collision' with an outer-scale g LF fluctuation lasting a time ∆t ∼ L + r, LF /v A (Lithwick et al. 2007).

If χ -LF /lessmuch 1 , then each such collision causes only a small perturbation to the outerscale z -fluctuation, and the turbulence is weak. In this limit, the effects of successive collisions add like a random walk, and roughly

$$N = ( \chi _ { L F } ^ { - } ) ^ { - 2 }$$

collisions are needed for nonlinear interactions to cause an order-unity change in the outer-scale z -fluctuation. The outer-scale z -cascade time scale t -NL is then ∼ NL + r, LF /v A . The generation of outer-scale z -(or f ) fluctuations by nonWKB 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 A , the increment to f from non-WKB reflection during a time ∆t = L + r, LF /v A is of order

$$\Delta f & \sim \left ( \frac { U - v _ { \Lambda } } { 2 v _ { \Lambda } } \right ) \left | \frac { d v _ { \Lambda } } { d r } \right | g _ { \text{LF,rms} } \Delta t. & & ( 4. 1 1 ) \\ ( 2. 1 5 ) & \text{ that the corresponding increment to $z^{-}$ is of order}$$

∣ ∣ It follows from (2.15) that the corresponding increment to z -is of order

$$\Delta z ^ { - } & \sim \left ( \frac { U + v _ { A } } { 2 v _ { A } } \right ) \left | \frac { d v _ { A } } { d r } \right | z _ { L F, r m s } ^ { + } \Delta t. & ( 4. 1 2 ) \\ \intertext { e of $z^{-}$ is approximately the ``amount`` of $z^{-}$ that ``builds up`` in frame $S^{-}$} \text{of} \text{tion} \text{ during the cascade} \text{-} \text{$m} \text{$time}$ scale} \text{-} N \text{At} \text{-} The resulting}$$

$$z _ { r m s } ^ { - } & \sim \frac { L _ { \perp } } { A } \left ( \frac { U + v _ { A } } { 2 v _ { A } } \right ) \left | \frac { d v _ { A } } { d r } \right |. \\ \intertext { n the outer-scale $z^{-}$ fluctuations are sheared coherently throughout}$$

∣ ∣ The r.m.s. value of z -is approximately the 'amount' of z -that 'builds up' in frame S -by non-WKB reflection during the cascade/damping time scale ∼ N∆t . The resulting value of z -rms is ∼ N 1 / 2 ∆z -, or, equivalently,

∣ ∣ If χ -LF /greaterorsimilar 1 , then the outer-scale z -fluctuations are sheared coherently throughout

their lifetimes, the turbulence is strong, and t -NL ∼ L ⊥ / z ( + LF rms , A ) . In this case, z -rms is approximately the rate at which z -fluctuations are produced by non-WKB reflection multiplied by t -NL , which again leads to (4.13). This estimate, with A → 1 , is the same as that obtained by Chandran &amp; Hollweg (2009) for the strong-turbulence limit. In the limits U → 0 and A → 1 , (4.13) is also the same as the estimate by Dmitruk et al. (2002) for the strong-turbulence limit.

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

$$z _ { r m s } ^ { - } = \frac { c ^ { - } L _ { \perp } } { A } \left ( \frac { U + v _ { A } } { 2 v _ { A } } \right ) \left | \frac { d v _ { A } } { d r } \right | \\ \text{value of } \chi _ { L F } ^ { - }, \, \text{where } c ^ { - } \text{ is a dimensionless piecewise constant function} \\ \quad \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \, \sim \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \colon \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \dodelta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \ delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \Dod \, \delta \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dot \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \delta \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \ dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \ dod \, \dot \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \ dod \, \dodo \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \ddod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \,\dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \, \dod \,$$

∣ ∣ regardless of the value of χ -LF , where c -is a dimensionless piecewise constant function that has one value at r &gt; r m and a smaller value at r /lessorequalslant r m , where r m is defined in (3.12). Before discussing c -further, we note an immediate consequence of (4.14), that z -rms increases as A (or equivalently sin θ ) decreases. This is because reducing sin θ decreases the rate at which outer-scale 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 &lt; r m

The reason we take c -to have a smaller value at r &lt; r m than at r &gt; r m is that the non-WKB-reflection source term for z -fluctuations reverses direction at r = r m , since d v A / d r changes sign. Since g LF has a large radial correlation length, when z -fluctuations produced via non-WKB reflection at r &gt; r m propagate to r &lt; r m , they tend to cancel out the z -fluctuations that are produced via non-WKB reflection at r &lt; r m , reducing z -rms . If the z -fluctuations at r = r m can propagate a radial distance ∼ ( r m -R /circledot ) before cascading and dissipating, then this cancellation effect is large. On the other hand, if the z -fluctuations at r = r m can only propagate a radial distance /lessmuch ( r m -R /circledot ) before cascading and dissipating, then little cancellation occurs. To account for this phenomenology, we set

$$c ^ { - } = \begin{cases} \ c _ { 1 } ^ { - } & \quad \text{if $r\leqslant r_{m}$} \\ \ c _ { 0 } ^ { - } & \quad \text$$

where c -O is a dimensionless free parameter,

$$c _ { 1 } ^ { - } = \frac { c _ { O } ^ { - } } { 1 + M },$$

$$M = \frac { v _ { \text{Am} } L _ { \perp m } } { z _ { \text{LFm} } ^ { + } L _ { \nabla } }, \quad \quad \quad$$

L ∇ is a free parameter with dimensions of length, and v Am , L ⊥ m , and z + LFm are the values of v A , L ⊥ , and z + LF rms , at r = r m . As we argue below (see (4.20)), A is of order unity at r &lt; r m , which means that M is the approximate radial distance an outer-scale z -fluctuation at r = r m propagates before cascading to smaller scales, divided by L ∇ . We can rewrite M in terms of quantities evaluated at r = r b by making the approximations that v A /greatermuch U at r /lessorequalslant r m and g LF rms , ( r m ) /similarequal g LF rms , ( r b ) and by using (2.17). This yields

$$M & = \frac { v _ { \text{Ab} } L _ { \perp b } } { z _ { \text{LFb} } ^ { + } L _ { \nabla } } \left ( \frac { v _$$

where v Ab , L ⊥ b , and z + LFb are the values of v A , L ⊥ , and z + LF rms , at r = r 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 z -anti-parallel to z + (i.e., an excess of magnetic energy over kinetic energy) (Müller &amp; Grappin 2005; Boldyrev et al. 2011). At r &gt; r m , d v A / d r &lt; 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 &lt; r m , d v A / d r &gt; 0 , and non-WKB reflection acts to produce positive residual energy. In other words, at r &lt; r m , linear processes (non-WKB reflection) and non-linear processes have competing effects on the alignment of z -. Based on these arguments, we conjecture that at r &lt; r m the outer-scale fluctuations do not develop significant alignment, and that at r &gt; r m the outer-scale z + LF and z -fluctuations become increasingly aligned as the z + LF fluctuations 'decay' via nonlinear interactions. We also conjecture that A is a decreasing function of τ (ph) v , because a larger τ (ph) v increases the efficiency of non-WKB reflection, which produces z -fluctuations that are aligned with z + LF . In addition, we conjecture that A is a decreasing function of

$$\Gamma = \frac { z _ { \text{LF,rms} } ^ { + } L _ { r, \text{LF} } ^ { + } } { L _ { \perp } v _ { \text{A}$$

which is the critical-balance parameter χ -LF in (4.9) without the factor of A . There are two reasons for taking A to decrease with increasing Γ . The first is that when Γ /lessmuch 1 , outer-scale z -fluctuations can propagate through many different outer-scale z + LF fluctuations before cascading to smaller scales. The z -fluctuations that are co-located with a particular outer-scale z + LF 'eddy' of radial extent ∼ L + r, LF are thus a mixture of the z -fluctuations produced by the non-WKB reflection of that z + LF eddy and z -fluctuations that were initially produced by the non-WKB reflection of z + LF eddies located farther from the Sun. The greater the number of distinct outer-scale z + LF eddies whose reflections contribute to the value of z -at any single point, the less aligned the z -field will be with any individual z + LF eddy. Moreover, when Γ &gt; 1 , shearing of the z -fluctuations by z + LF rotates the z -fluctuations into alignment with z + LF , and the resulting value of A is a decreasing function of Γ (Chandran et al. 2015 ). We quantify b the foregoing conjectures by setting

$$A = \begin{cases} \ A _ { 0 } & \text{ if } r < r _ { m } \\ A _ { 0 } \left [ 1 + \frac { \tau _ { v } ( \text{ph} ) } { \tau _ { \theta } } \ln \left ( \frac { g _ { L F m } ^ { 2 } } { g _ { L F, r m s } ^ { 2 } } \right ) \right ] ^ { - 1 } & \text{ if } r > r _ { m } \end{cases}, \quad ( 4. 2 0 )$$

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 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 z + fluctuations to behave in the same way, setting

$$\frac { L _ { r, L F } ^ { + } } { L _ { r, L F b } ^ { + } } = \frac { U + v _ { \text{A} } } { U _ { \text{b} } + v _ { \text{Ab} } },$$

where L + r, LFb , U b , and v Ab are the values, respectively, of L + r, LF , U , and v A evaluated at

r = r b , and likewise for L + r, HF . It then follows from (2.17), (3.3), (3.23), and (4.21) that

$$\Gamma = \Gamma _ { \text{b} } \frac { g _ { \text{LF}, \text{rms} } } { g _ { \text{LFb} } } \sqrt { \frac { v _ { \text{A} } } { v _ { \text{Ab} } } }, \quad \text{(4.22)} \\ \epsilon \Gamma \sim \Delta \sim - \Delta$$

where Γ b is the value of Γ at r = r b .

## 4.4. Solving for the fluctuation-amplitude profiles

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

$$c _ { I } \equiv c _ { \text{diss} } c _ { I } ^ { - } \quad c _ { O } \equiv c _ { \text{diss} } c _ { O } ^ { - }.$$

$$\text{comming} \, \left \{ \mathcal { U } \right \}, \, \left \{ \mathcal { U } \right \}, \, \left \{ \mathcal { U } \right \}, \, \left \{ \mathcal { U } \right \}, \, \left \{ \mathcal { U } \right \}, \, \right \}, \, \text{and} \, \left \{ \mathcal { U } \right \}, \, \text{we youam} \\ \frac { d } { d r } \ln g _ { L F, r m s } ^ { 2 } = \begin{cases} \begin{array} { c c } - c _ { I } \frac { d } { d r } \ln v _ { A } & \text{if $r\leqslant r_{m}$} \\ \frac { d } { d r } \frac { d } { d r } \ln v _ { A } & \text{if $r>r_{m}$} \end{array} \end{cases}, \quad \text{(4.23)} \\ \quad \text{$c_{I} \equiv c_{diss} c_{I}^{-}$} \quad c _ { O } \equiv c _ {diss} c_{O}^{-}$}.$$

After integrating (4.23), we find that

$$& \text{er integrating } ( 4. 2 3 ), \, \text{we nnd that} \\ & \frac { g _ { \text{LF,rms} } ^ { 2 } } { g _ { \text{LFb} } ^ { 2 } } = \begin{cases} \left ( \frac { v _ { \text{Ab} } } { v _ { \text{A} } } \right ) ^ { c _ { \text{I} } } & \text{if $r_{b} < r<r_{m}$} \\ \left ( \frac { v _ { \text{Ab} } } { v _ { \text{Am} } } \right ) ^ { c _ { \text{I} } } \left ( \frac { v _ { \text{A} } } { v _ { \text{Am} } } \right ) ^ { c _ { \text{Co} } } & \text{if $r > r_{m}$} \end{cases}, \quad & ( 4. 2 5 ) \\ & \text{ere $g_{LFb}$ is the value of $g_{LF,rms} at $r=r_{b}$. Upon combining } ( 4. 4 ), \, ( 4. 4 ), \, \text{and} \\ & \text{$15$} ), \, \text{we obtain}$$

 where g LFb is the value of g LF rms , at r = r b . Upon combining (4.4), (4.5), (4.14), and (4.15), we obtain

$$\frac { d } { d r } \ln g _ { \text{HF,rms} } ^ { 2 } = \begin{cases} \begin{array} { c c c } - \frac { c _ { \text{I} } } { A } \frac { d } { d r } \ln v _ { \text{A} } & \text{ if } r \leqslant r _ { \text{m} } \\ \frac { c _ { \text{O} } } { A } \frac { d } { d r } \ln v _ { \text{A} } & \text{ if } r > r _ { \text{m} } \end{array} \end{cases} \quad \begin{pmatrix} 4. 2 6 \\ \end{pmatrix} \\ e \text{ aid of } ( 4. 2 0 ) \text{ and } ( 4. 2 2 ), \, \text{ we integrate } ( 4. 2 6 ) \text{ to obtain } \begin{array} { c c c } \end{array} \quad \begin{array} { c c c } \end{cases}$$

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

$$\text{$+$ the and of} \, ( 4. 2 V ) \, \text{and} \, ( 4. 2 Z ), \, \text{we integrate} \, ( 4. 2 0 ) \, \text{to obtain} \\ \frac { g _ { \text{HF} } ^ { 2 }, _ { \text{rms} } } { g _ { \text{HF} b } ^ { 2 } } = \begin{cases} \begin{array} { c c } \left ( \frac { v _ { \text{Ab} } } { v _ { \text{A} } } \right ) ^ { c _ { 1 } / A _ { 0 } } & \text{if $r<r_{m}$} \\ \left ( \frac { v _ { \text{Ab} } } { v _ { \text{Am} } } \right ) ^ { c _ { 1 } / A _ { 0 } } w ^ { c o / A _ { 0 } } e ^ { - H } & \text{if $r>r_{m}$} \end{cases} \\ \text{e} \\ \text{$1$} \, - \, c _ { \text{O} } ^ { 2 } \text{$L_{b}$} \left ( v _ { \text{Am} } \right ) ^ { ( 1 - c _ { 1 } ) / 2 } / \underset { 1 \sim \dots \pi \ 1 \cdots \dots \dots \pi + 1 \vee \infty } \, \end{cases}$$

$$H = \frac { c _ { \theta } c _ { \mathcal { O } } ^ { 2 } \Gamma _ { \text{b} } } { \sigma ^ { 2 } A _ { 0 } } \left ( \frac { v _ { \text{Am} } } { v _ { \text{Ab} } } \right ) ^ { ( 1 - c _ { 1 } ) / 2 } ( \sigma w ^ { \sigma } \ln w - w ^ { \sigma } + 1 ), \quad \ \ ( 4. 2 8 ) \\ \nu. \quad 1 \pm c _ { \sim }$$

$$w = \frac { v _ { \mathrm A } } { v _ { \text{Am} } } \quad \sigma = \frac { 1 + c _ { \mathrm O } } { 2 },$$

g HFb is the value of g HF rms , at r = r b , and c θ = τ (ph) v /τ θ .

4.5. Turbulent-heating rate

The turbulent-heating rate in our model is

$$Q = \frac { \supset } { 2 } \left [ \gamma _ { \text{HF} } ^ { + } ( z _ { \text{HF}, \text{rms} } ^ { + } ) ^ { 2 } + \gamma _ { \text{LF} } ^ { + } ( z _ { \text{LF}, \text{rms} } ^ { + } ) ^ { 2 } + \gamma ^ { - } ( z _ { \text{rms} } ^ { - } ) ^ { 2 } \right ],$$

where

$$\gamma ^ { - } = \gamma ^ { - } _ { L F } + \gamma ^ { - } _ { H F }$$

is the cascade rate of the outer-scale z -fluctuations, and γ -LF ( γ -HF ) is the contribution to γ -from interactions between z -fluctuations and LF (HF) z + fluctuations. To allow for either weakly turbulent ( χ -LF &lt; 1 ) or strongly turbulent ( χ -LF /greaterorequalslant 1 ) shearing of z -fluctuations by z + LF fluctuations, we set

$$\gamma _ { L F } ^ { - } = \frac { c _ { \text{diss} } z _ { L F, \text{rms} } ^ { + } A } { L _ { \perp } } \times \begin{cases} \chi _ { L F } ^ { - } & \text{ if } \chi _ { L F } ^ { - } \leq 1 \\ 1 & \text{ if } \chi _ { L F } ^ { - } > 1 \end{cases}.$$

In analogy to (4.9), we define the critical-balance parameter for the shearing of z -fluctuations by z + HF fluctuations to be

$$\chi _ { \text{HF} } ^ { - } = \frac { z _ { \text{HF}, \text{rms} } ^ { + } L _ { r, \text{HF} } ^ { + } } { L _ { \perp } v _ { \text{A} } },$$

where we have omitted the factor of A , because we take z -to be aligned with z + LF but not with z + HF . We then set

$$\gamma _ { \text{HF} } ^ { - } = \frac { c _ { \text{diss} } z _ { \text{HF}, \text{rms} } ^ { + } } { L _ { \perp } } \times \begin{cases} \ \chi _ { \text{HF} } ^ { - } & \text{ if } \chi _ { \text{HF} } ^ { - } \leq 1 \\ 1 & \text{ if } \chi _ { \text{HF} } ^ { - } > 1 \end{cases}.$$

## 4.6. Comparison with simulation results

To compare our model with one of our numerical simulations, we treat L ⊥ ( r b ) = L box ( r b ) / 3 , L + r, LF ( r b ) , L + r, HF ( r b ) , and z + rms ( r 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

$$\frac { g _ { \text{HFb} } ^ { 2 } } { g _ { \text{LFb} } ^ { 2 } } = 1$$

in all our model solutions. We take L + r, HF ( r b ) in our simulations to be the radial separation ∆r at which C r , ∆r ( b ) = 1 / 2 , where

$$C ( r, \Delta r ) & = \frac { \langle g ( \tilde { x }, \tilde { y }, r, t ) \cdot g ( \tilde { x }, \tilde { y }, r + \Delta r, t ) \rangle } { \langle | g ( \tilde { x }, \tilde { y }, r, t ) | ^ { 2 } \rangle } \quad & ( 4. 3 6 ) \\ \text{ial autocorrelation function of the $a$ functions and $\tilde{$r$ and $\tilde{$r$ are}$}$$

We take the free parameters c diss , c -O , τ θ , A 0 , and L ∇ 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.

is the radial autocorrelation function of the g fluctuations, and ˜ x and ˜ y are defined following (3.28). On the other hand, because figure 4 shows that the LF fluctuations are energetically dominant at r = r max , we define L + r, LF ( r max ) to be the value of ∆r at which C r ( max , -∆r ) = 1 2 / . Applying (4.21), we then set L + r, LF ( r b ) = L + r, LF ( r max )( U b + v Ab ) / U [ ( r max ) + v A ( r max )] = 0 886 . L + LF ( r max ) . The values of L box ( r b ) , L + r, LF ( r b ) , L + r, HF ( r b ) , and z + rms ( r b ) in our three simulations are listed in table 3.

Figures 3 through 6 show the radial profiles of z + rms , z -rms , z + HF rms , , z + LF rms , , sin θ , 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 + HF rms , /z + LF rms , and

Table 3. Boundary Conditions in Our Analytic Model for Matching Runs 1 through 3

| Quantity                                                        | Run 1                                                             | Run 2                                                            | Run 3                                                           |
|-----------------------------------------------------------------|-------------------------------------------------------------------|------------------------------------------------------------------|-----------------------------------------------------------------|
| z + rms ( r b ) L + r, HF ( r b ) L + r, LF ( r b ) L ⊥ ( r b ) | 61 km/s 0 . 015 R /circledot 0 . 071 R /circledot 1 . 4 × 10 3 km | 55 km/s 0 . 015 R /circledot 0 . 27 R /circledot 1 . 4 × 10 3 km | 41 km/s 0 . 10 R /circledot 0 . 35 R /circledot 5 . 3 × 10 3 km |

z + rms is the r.m.s. amplitude of the outward-propagating Elsasser variable, L + r, LF is the radial correlation length of the low-frequency outward-propagating HeinemannOlbert variable g LF , L + r, HF is the radial correlation length of the high-frequency outward-propagating Heinemann-Olbert variable g HF , L ⊥ is the perpendicular outer scale (see (3.23)), and r b is the radius of the coronal base defined in (3.10).

sin θ decrease with increasing r , particularly in Run 2. The radial decrease in z + HF /z + LF is consistent with our expectation that high-frequency z + fluctuations cascade and dissipate more rapidly than low-frequency z + fluctuations, because high-frequency z + fluctuations are not aligned with z -. In our model, the radial decrease in sin θ is related both to the comparatively rapid cascade of the unaligned high-frequency z + fluctuations and the fact that the low-frequency z + fluctuations become increasingly aligned with z -as they interact nonlinearly with z -. We note that the decrease in sin θ coincides with an increase in z -rms for the reasons described following (4.14). The model reproduces the z ± rms 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 &gt; 3 R /circledot . The most notable failing of the model is that z -rms = Q = 0 at r = r m , because our estimate of z -rms is proportional to the local value of d v A / d , r which vanishes at r = r m . A more realistic model would account for the fact that z -fluctuations propagate a finite distance before cascading and dissipating, which would smooth out the profiles of z -rms and Q in the vicinity of r = r 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 ± rms 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. Previous studies of the Elsasser power spectra in MHD turbulence

In this section, we review previous studies of the Elsasser power spectra in homogeneous RMHD turbulence and reflection-driven MHD turbulence. The reader already familiar with this literature may wish to skip directly to Section 6. We follow the convention of describing the turbulent z ± fluctuations as collections of AW packets, using λ to denote the length scale of a wave packet measured perpendicular to the magnetic field, l ± λ to denote the correlation length measured along the magnetic field of z ± wave packets with

Table 4. Best-Fit Free Parameters in our Analytic Model

| Parameter                | Value                                         |
|--------------------------|-----------------------------------------------|
| c diss c - O τ θ A 0 L ∇ | 0 . 2 1 . 8 3 . 2 min 0.6 0 . 15 R /circledot |

The quantity c diss is a coefficient appearing in the cascade/damping rates γ + HF and γ + LF ((4.5) and (4.7)), c -O is a coefficient in our estimate of z -rms (see (4.14) through (4.16)), τ θ and A 0 are constants appearing in our estimate of the alignment angle ((4.20)), and L ∇ is a length scale that affects the degree to which z -fluctuations produced by non-WKB reflection at r &gt; r m /similarequal 1 7 . R /circledot cancel out the z -fluctuations produced by non-WKB reflection at r &lt; r m (see (4.14) through (4.17)).

perpendicular scale λ , and δz ± λ to denote the amplitude of wave packets at scale λ - i.e., the r.m.s. increment in z ± across a distance λ perpendicular to the magnetic field.

## 5.1. Balanced, homogeneous RMHD turbulence

In 'balanced turbulence,' the statistical properties of z + and z -fluctuations are identical. In particular,

$$\delta z _ { \lambda } ^ { + } = \delta z _ { \lambda } ^ { - } \quad l _ { \lambda } ^ { + } = l _ { \lambda } ^ { - },$$

and the cross helicity (the difference between the energies per unit mass of z + and z -fluctuations) is zero. In homogeneous RMHD turbulence, the strongest nonlinear interactions are local in scale, meaning that δz ± λ fluctuations are cascaded primarily by z ∓ fluctuations at perpendicular scales comparable to λ . To understand how a δz ± λ 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. 2007). This slice 'collides' with a series of counterpropagating δz ∓ λ wave packets. Each collision has a duration of

$$t _ { \lambda } ^ { \pm } & \sim \frac { l _ { \lambda } ^ { \mp } } { v _ { \mathrm A } }.$$

The instantaneous rate at which δz ∓ λ wave packets shear δz ± λ wave packets is ∼ δz ∓ λ /λ . During a single collision, the aforementioned 'slice' of the δz ± λ fluctuation undergoes a fractional distortion of order (Goldreich &amp; Sridhar 1995; Goldreich &amp; Sridhar 1997; Lithwick et al. 2007)

$$\chi _ { \lambda } ^ { \pm } = \frac { \delta z _ { \lambda } ^ { \mp } l _ { \lambda } ^ { \mp } } { \lambda v _ { A } }.$$

## 5.1.1. Weak balanced turbulence

If χ ± λ /lessmuch 1 , a δz ± λ wave packet undergoes only a small fractional change during each collision, and the turbulence is weak. Ng &amp; Bhattacharjee (1996, 1997) and

Goldreich &amp; Sridhar (1997) 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 δz ± λ wave packet is ∼ N 1 / 2 χ ± λ . After N ∼ ( χ ± -λ ) 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

$$\tau _ { \lambda } ^ { \pm } & \sim \left ( \chi _ { \lambda } ^ { \pm } \right ) ^ { - 2 } t _ { \lambda } ^ { \pm } \sim \frac { \lambda ^ { 2 } v _ { A } } { ( \delta z _ { \lambda } ^ { \mp } ) ^ { 2 } l _ { \lambda } ^ { \mp } }. \\ \ e \, \delta z _ { \lambda } ^ { + } \text{ nor } \delta z _ { \lambda } ^ { - } \text{ wave packet is altered significantly during any single}$$

Because neither the δz + λ nor δz -λ wave packet is altered significantly during any single collision, the leading and trailing edges of a δz ± λ 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 λ (Shebalin et al. 1983); i.e.,

$$l _ { \lambda } ^ { \pm } & \in \lambda ^ { 0 }. \\ &. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad.$$

In the inertial range, the z ± energy-cascade rate (per unit mass), /epsilon1 ± , is independent of scale:

$$\epsilon ^ { \pm } \sim \frac { ( \delta z _ { \lambda } ^ { \pm } ) ^ { 2 } } { \tau _ { \lambda } ^ { \pm } } \in \lambda ^ { 0 }. \quad \quad \quad ( 5. 6 )$$

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

$$\begin{matrix} & \delta z _ { \lambda } ^ { \pm } \in \lambda ^ { 1 / 2 }. & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & \\ & \cdot & \cdot & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & && & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & &  & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & &
 & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & \ & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & - & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & &. & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & &   & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & . & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & &
 \end \ & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & 
 & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & &.. & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & &. \ & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & &.
 & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & &
 { & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & } & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & \\ & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & = & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & &$$

The scaling of the one-dimensional power spectrum of the z ± fluctuations, denoted E ± ( k ⊥ ) , follows from the relation

$$k _ { \perp } E ^ { \pm } ( k _ { \perp } ) \sim \left ( \delta z _ { \lambda } ^ { \pm } \right ) ^ { 2 } \Big | _ { \lambda = k _ { \perp } ^ { - 1 } }, \\ \text{ponent of the wave vector perpendicular to the background magnetic} \quad, \quad, \quad.$$

∣ where k ⊥ is the component of the wave vector perpendicular to the background magnetic field. Equations (5.7) and (5.8) imply that

$$\begin{array} {$$

The scaling in (5.9) has been found in direct numerical simulations (Perez &amp; Boldyrev 2008) as well as in exact solutions to the weak-turbulence wave kinetic equations for incompressible MHD turbulence (Galtier et al. 2000). It is worth noting, however, that in weak-turbulence theory all AWs are cascaded by k ‖ = 0 modes, where k ‖ is the wave vector component along 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 &amp; Perez 2009; Schekochihin et al. 2012; Meyrand et al. 2015), as discussed further in Section 5.2.1.

## 5.1.2. Strong balanced turbulence

If χ ± λ /greaterorsimilar 1 , then each slice of a δz ± λ wave packet is strongly distorted during a single collision, the turbulence is strong, and z ± energy at scale λ cascades to smaller scales on the time scale

$$\tau _ { \lambda } ^ { \pm } & \sim \frac { \lambda }$$

leading to a scale-independent energy-cascade rate

$$\epsilon ^ { \pm } \sim \frac { ( \delta z _ { \$$

Equations (5.1) and (5.11) imply that

$$\delta z _ { \lambda } ^ { \pm } \in \lambda ^ { 1 /$$

which implies via (5.8) that

$$E ^ { \pm } ( k _ { \perp } ) \in k _ { \per$$

Goldreich &amp; Sridhar (1995) 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.,

$$\chi _ { \lambda } ^ { \pm } \sim 1. \quad \quad \quad$$

Numerical simulations confirm that this 'critical-balance' conjecture describes strong RMHD turbulence not only on average (Cho &amp; Vishniac 2000), but structure by structure (Mallet et al. 2015). Together, (5.12) and (5.14) imply that

$$l _ { \lambda } ^ { \pm } \in \lambda ^ { 2 / 3 }.$$

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 -3 / 2 ⊥ , because of scale-dependent dynamic alignment (Boldyrev 2005, 2006; Mason et al. 2008; Perez et al. 2012) and/or intermittency (Maron &amp; Goldreich 2001; Chandran et al. 2015 b ; Mallet &amp; Schekochihin 2017). On the other hand, Beresnyak (2012, 2014) argued for a scaling closer to k -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 Mallet et al. (2017 b a , ) and Loureiro &amp; Boldyrev (2017 b a , ), who investigated the disruption of sheet-like structures in RMHD turbulence by the tearing instability and magnetic reconnection (see also Pucci &amp; Velli 2014; Pucci et al. 2018; Vech et al. 2018).

## 5.2. Imbalanced RMHD turbulence in homogeneous plasmas

In 'imbalanced turbulence,' one of the Elsasser variables, say z + , has a substantially higher r.m.s. amplitude than the other:

$$z _ { r m s } ^ { + } > z _ { r m s } ^ { - }.$$

Equation (5.16) includes the highly imbalanced case, in which z + rms /greatermuch z -rms , as well as moderately imbalanced turbulence, in which, e.g., z + rms /similarequal 2 z -rms .

## 5.2.1. Weak imbalanced turbulence

When (5.16) is satisfied and

$$\chi ^ { + } _ { \lambda } \ll 1 \quad \chi ^ { - } _ {$$

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

$$\alpha ^ { + } + \alpha ^ { - } = 4,$$

where

$$E ^ { \pm } ( k _ { \perp } ) \in k _ { \perp } ^ { - \alpha ^ { \pm } }, & & ( 5. 1 9 ) \\ \intertext { co version of } \, ( 3 \, 3 \cdots \, I \text{$$

the homogeneous-turbulence version of (3.35). Lithwick &amp; Goldreich (2003) argued that in weak incompressible MHD turbulence, the spectra are 'pinned' at the dissipation

wavenumber k ⊥ d , with E + ( k ⊥ d ) = E -( k ⊥ d ) , and that the more energetic Elsasser variable has the steeper inertial-range power spectrum. Boldyrev &amp; Perez (2009) espoused a different picture, in which a 'condensate' of magnetic fluctuations at k ‖ = 0 dominates the energy cascade, leading to a state in which α + = α -= 2 . Schekochihin et al. (2012) developed a theory accounting for both weakly turbulent AWs with nonzero k ‖ and 2D modes with k ‖ = 0 , and found that α + = α -= 2 for the weakly turbulent modes and α + = α -= 1 for the 2D modes in the imbalanced case.

## 5.2.2. Strong imbalanced turbulence

When (5.16) is satisfied and χ + λ or χ -λ is /greaterorsimilar 1 , the turbulence is considered strong. A number of authors have developed models of strong imbalanced MHD turbulence (e.g., Beresnyak &amp; Lazarian 2008; Chandran 2008a; Beresnyak &amp; Lazarian 2009; Perez &amp; Boldyrev 2009, 2010; Podesta &amp; Bhattacharjee 2010). Here we focus on the study by Lithwick et al. (2007) (hereafter LGS), who explored an assumption about the forcing of outer-scale 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

$$\chi _ { \lambda } ^ { - } \gtrsim 1.$$

Equation (5.20) implies that δz -λ fluctuations are sheared on a time scale λ/δz + λ that is comparable to or less than the time l + λ /v A for a slice of a δz -λ wave packet to pass through a counter-propagating δz + λ wave packet. The cascade time scale for δz -λ wave packets is therefore

$$\tau _ { \lambda } ^ { - } \sim \frac { \lambda } { \delta z _ { \lambda } ^ { + } }. \quad \ \ \ \ \ ( 5. 2 1 ) \\. \quad \ \. \quad. \quad. \quad. \quad. \$$

LGS argued that, since a δz -λ wave packet cascades after it has propagated along the background magnetic field for a distance ∼ v A τ -λ , the parallel correlation length of the δz -λ wave packet is

$$l _ { \lambda } ^ { - } \sim v _ { A } \tau _ { \lambda } ^ { - } \sim \frac { v _ { \Lambda } \lambda } { \delta z _ { \lambda } ^ { + } }. \quad \quad \quad$$

LGS further argued that, since δz + λ wave packets separated by a distance l -λ along the magnetic field are sheared by uncorrelated δz -λ fluctuations,

$$l _ { \lambda } ^ { + } \sim l _ { \lambda } ^ { - }. & & ( 5. 2 3 ) \\. &. &. &. &.$$

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

$$\chi _ { \lambda } ^ { - } \sim 1 \quad \chi _ { \lambda } ^ { + } \sim \frac { \delta z _ { \lambda } ^ { - } } { \delta z _ { \lambda } ^ { + } } < 1.$$

The apparent implication of the second half of (5.24), particularly when δz -λ /δz + λ /lessmuch 1 , is that δz + λ wave packets cascade in a weakly turbulent manner, through multiple, uncorrelated collisions with δz -λ wave packets, each of which leads to a small fractional change in the δz + λ wave packet of order χ + λ (see Section 5.1.1). LGS argued, however, that each δz + λ wave packet is in fact sheared coherently throughout its lifetime, even when χ + λ ∼ δz -λ /δz + λ /lessmuch 1 . To establish this conclusion, LGS considered the ' z + frame,' which moves with z + fluctuations at speed v A along B 0 relative to the background plasma. They then proposed a thought experiment in which the amplitude of z -is infinitesimal, so that z -has negligible effect upon z + . The z + vector field is then time-independent in the z + frame. If the z + fluctuations are initialized with a power-law spectrum spanning

the entire inertial range, and if z -fluctuations are continuously injected at the outer scale with an arbitrarily long coherence time T in the z + frame, then the z -fluctuations will cascade to small scales and set up not just a statistical steady state, but an actual steady state in the z + frame in which the z -vector field is independent of time. This latter conclusion follows because z -is nonlinearly distorted by z + , which is constant in time in the z + frame. The δz -λ fluctuations encountered by a δz + λ wave packet are therefore coherent for an arbitrarily long time, and in particular for a time much longer than the crossing time

$$t _ { \text{cross}, \lambda } ^ { + } \sim \frac { l _ { \lambda } ^ { - } } { v _ { A } }$$

required for a slice of the δz + λ wave packet to propagate through a δz -λ wave packet.

Building upon this thought experiment, LGS proceeded to consider the more realistic case in which δz -λ is finite, but still small compared to δz + λ at all λ . They made a key assumption, which we call the 'coherence assumption,' that the coherence time T (at a fixed position in the z + frame) of the forcing of outer-scale z -fluctuations is at least as long as the 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 δz -λ fluctuations encountered by a δz + λ wave packet is the variation of the z + vector field, not the crossing of counter-propagating wave packets, and the δz + λ wave packet is sheared coherently throughout its lifetime. The z + cascade time scale at scale λ then becomes

$$\tau _ { \lambda } ^ { + } & \sim \frac { \lambda } { \delta z _ { \lambda } ^ { - } },$$

and

$$\epsilon ^ { + } \sim \frac { ( \delta z _ { \lambda } ^ { + } ) ^ { 2 } } { \tau _ { \lambda } ^ { + } } \sim \frac { ( \delta z _ { \lambda } ^ { + } ) ^ { 2 } \delta z _ { \lambda } ^ { - } } { \lambda }.$$

Because of (5.24), τ -λ ∼ λ/δz + λ , and

$$\epsilon ^ { - } \sim \frac { ( \delta z _ { \lambda } ^ { - } ) ^ { 2 } } { \tau _ { \lambda } ^ { - } } \sim \frac { ( \delta z _ { \lambda } ^ { - } ) ^ { 2 } \delta z _ { \lambda } ^ { + } } { \lambda }.$$

Setting /epsilon1 ± ∝ λ 0 , LGS combined Equations (5.27) and (5.28) to obtain

$$\delta z _ { \lambda } ^ { + } \in \delta z _ { \lambda } ^ { - } \in \lambda ^ { 1 / 3 },$$

which, via (5.8), implies that

$$E ^ { \pm } ( k _ { \perp } ) \in k _ { \perp } ^ { - 5 / 3 }.$$

## 5.3. Anomalous coherence in reflection-driven MHD turbulence

Velli et al. (1989) (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 ± ( k ) . They divided the z ± fluctuations into 'primary' and 'secondary' components, where the primary components of z ± had the usual phase velocities of U ± v A . The secondary components of z ± were driven modes produced by the reflection of z ∓ fluctuations and as a consequence had phase velocities of U ∓ v A . VGM considered the super-Alfvénic solar wind at r &gt; r A and took z -to be dominated by secondary fluctuations. VGM estimated the r.m.s. amplitude of the secondary component of z -at

scale k -1 , which we denote z -k, s , to be

$$z _ { k, s } ^ { - } & \sim \frac { z _ { k, p } ^ { + } } { k v _ { A } \tau _ { r } },$$

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

$$\epsilon ^ { + } \sim k z _ { k, s } ^ { - } ( z _ { k, p } ^ { + } ) ^ { 2 } \sim \frac { ( z _ { k, p } ^ { + } ) ^ { 3 } } { v _ { \mathrm A } \tau _ { r } } \quad \cdots \quad \cdots \quad \cdots \cdots \cdots \cdots \cdots \cdots$$

and took /epsilon1 + to be independent of k , obtaining z + k, p ∝ k 0 . Equations (5.8) and (5.31) then yield

$$E ^ { + } ( k ) \in k ^ { - 1 } \quad E ^ { - } ( k ) \in k ^ { - 3 }. \quad \text{(5.33)} \\ \text{macro } & \text{the VCM modal with the } \text{$1\, C\,$\,$modal$ decoded in $1000 $5$} 9$$

It is useful to compare the VGM model with the LGS model discussed in Section 5.2.2. In both models, the z -fluctuations are anomalously coherent in the reference frame of the z + fluctuations. In the LGS model, this coherence is introduced via the 'coherence assumption' discussed in Section 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 z -fluctuations that are produced as secondary z -fluctuations cascade to small scales. In the LGS model, these tertiary z -fluctuations are anomalously coherent in the z + reference frame and drive the Elsasser spectra towards a k -5 / 3 ⊥ scaling rather than a k -1 scaling.

## 5.4. Inverse cascade in reflection-driven MHD turbulence

van Ballegooijen &amp; Asgari-Targhi (2017) 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 /epsilon1 ± ( k ⊥ , r, t ) at which nonlinear interactions transfer z ± energy from perpendicular wavenumbers less than k ⊥ to perpendicular wavenumbers greater than k ⊥ (their Equation (17) divided by ρ , with R → L box / 2 , f ± ,k → φ ± ,k , and a → k ⊥ L box / 2 ),

$$\epsilon ^ { \pm } ( k _ { \perp }, r, t ) = \frac { 1 } { [ L _ { \text{box} } ( r ) ] ^ { 2 } } \sum _ { k _ { \per$$

where φ ± k is the Fourier transform (in x and y ) of the Elsasser stream function φ ± (defined such that z ± = ∇ φ ± × B 0 /B 0 ), and M lji is a dimensionless mode-coupling coefficient that depends upon k ⊥ l , k ⊥ j , and k ⊥ i , but not upon the mode amplitudes. They found that /epsilon1 + ( k ⊥ , r, t ) became negative across a broad range of k ⊥ within a modest range of radii just larger than the radius (or radii) at which d v A / d r changes signs - e.g., just beyond the Alfvén-speed maximum at r = r 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 m surface at r = r 1 and one just outside the r = r m surface at r = r 2 , such that | d v A / d r | was the same at the two radii. They noted that, because z + rms /greatermuch z -rms , z -fluctuations

cascade much more rapidly than z + fluctuations. There thus exists a range of values of r 2 -r 1 for which the time ∆t 12 required for z + fluctuations to propagate from r 1 to r 2 is small compared to the outer-scale z + cascade time scale but large compared to the outer-scale z -cascade time scale. For values of r 2 -r 1 in this range,

$$& \phi _ { +, k } ( r _ { 2 }, t ) \simeq \phi _ { +, k } ( r _ { 1 }, t - \Delta t _ { 1 2 } ) & & ( 5. 3 5 )$$

$$\phi _ { -, k } ( r _ { 2 }, t ) \simeq - \phi _ { -, k } ( r _ { 1 }, t - \Delta t _ { 1 2 } ), \quad \quad \quad$$

where ∆t 12 is the 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 z + fluctuations during their transit from r 1 to r 2 . Equation (5.36) follows from the change in sign of d v A / d r at r = r m . † Because φ -,k changes sign and φ + ,k remains almost unchanged, /epsilon1 + 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 /epsilon1 + to be positive at r 1 . At r 2 , /epsilon1 + thus becomes negative, indicating an inverse cascade.

van Ballegooijen &amp; Asgari-Targhi (2017) found that as the z + fluctuations propagate farther beyond r = r m , they gradually adjust to the new value of 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 ⊥ (at which the nonlinear time is short) and later at small k ⊥ . In one of their simulations, there is an inverse cascade of z + energy throughout the region between the Alfvén-speed maximum at 1 4 . R /circledot and an outer radius of r = 2 5 . R /circledot . In a second simulation, there is an inverse cascade between the Alfvén-speed maximum at 1 6 . R /circledot and an outer radius of 4 R /circledot . Since the outer-scale z + cascade time is comparable to the time required for v A to change by a factor of 2 in the z + reference frame (Dmitruk et al. 2002; Chandran &amp; Hollweg 2009), the above results indicate that the inverse cascade persists (in the z + reference frame) for a time comparable to the outer-scale z + energy-cascade time scale.

Although /epsilon1 + became negative between r /similarequal r m and r /similarequal 2 r m in the numerical simulations of van Ballegooijen &amp; Asgari-Targhi (2017), the energy-dissipation rate (computed from the dissipation terms added to the governing equations) decreased by only a factor of /similarequal 2 within the inverse-cascade region. The reason for this is that the direct-cascade region at r &lt; r m had already 'done the work' of transporting z + energy to large k ⊥ , and the inverse cascade between r /similarequal r m and r /similarequal 2 r m was unable to completely evacuate the highk ⊥ part of the spectrum.

## 6. The Elsasser power spectra in our numerical simulations

In our numerical simulations, α + and α -approach /similarequal 3 / 2 as r increases to values /greaterorsimilar r A , as illustrated in figure 7. These spectral indices are broadly consistent with the LGS model of strong imbalanced turbulence. As discussed in Section 5.2.2, the central assumption of the LGS model is the 'coherence assumption' - that outer-scale z -fluctuations are injected in a manner that remains coherent over the lifetime of the outer-scale z + fluctuations when viewed in the ' z + reference frame,' which moves along B 0 at the same velocity ( U + v A ) as the z + fluctuations. It is difficult, at least for us, to justify this assumption with any generality for homogeneous RMHD turbulence. However, the coherence assumption is often satisfied in reflection-driven MHD turbulence,

- † Because the z -fluctuations cascade very rapidly near r = r m , (2.20) can be solved approximately by balancing the last term on the left-hand side and the first term on the right-hand side. In this approximation, changing the sign of d v A / d r changes the sign of f .

because the outer-scale z -fluctuations are produced by the reflection of outer-scale z + fluctuations, and by definition these z + fluctuations remain coherent in the z + reference frame throughout their lifetimes. A second requirement of the LGS model is that χ -λ /greaterorsimilar 1 . This requirement is marginally satisfied at r /greaterorsimilar r A in all three simulations, as we will document in greater detail in a separate publication. The LGS model thus provides a credible explanation for the Elsasser power spectra at r /greaterorsimilar r A in Runs 1 through 3. The discrepancy between the predicted α ± = 5 3 / scaling and the measured α ± /similarequal 3 / 2 scaling may result from some combination of intermittency and scale-dependent dynamic alignment, as in homogeneous RMHD turbulence (see Section 5.1.2).

As discussed in Section 3.10 (see also Velli 1993; Réville et al. 2018), the steep Alfvénspeed gradient in the upper chromosphere acts as a high-pass filter. Highk ⊥ z + fluctuations, which have large nonlinear frequencies and hence large linear frequencies (Goldreich &amp; Sridhar 1995), can propagate through this region with minimal reflection. In contrast, lowk ⊥ z + fluctuations undergo strong non-WKB reflection as they propagate from the lower chromosphere to the transition region. This selective transmission accounts for the very small value of α + just above the transition region in Runs 1 and 2. The z + spectrum in Run 3 does not flatten in the same way, presumably because the nonlinear time scale is larger than in Runs 1 and 2 because of the larger value of L ⊥ , causing all the z + fluctuations in Run 3 to undergo significant reflection in the upper chromosphere.

As discussed in Section 5.4, van Ballegooijen &amp; Asgari-Targhi (2017) showed that the z + fluctuations undergo a transient inverse cascade at r m /lessorsimilar r /lessorsimilar 2 r m , where r m is the location of the Alfvén-speed maximum ( 1 71 . R /circledot in our simulations). This inverse cascade results from the change in sign of d v A / d r at r = r m , which reverses the direction of the fast-cascading z -fluctuations, which in turn reverses the sign of /epsilon1 + in (5.34). The tendency for z -fluctuations to reverse direction at r = r m destroys the anomalous coherence of the z -fluctuations in the z + reference frame near r = r m , making the LGS model inapplicable. We do not have a detailed theory for how the spectra should scale between r = r m and r = 2 r m in the presence of this inverse cascade, but the simulation results indicates that the z + spectrum flattens significantly in this region relative to the LGS prediction.

## 7. Other parameter regimes and lack of universality

One of the principal sources of uncertainty in modeling MHD turbulence in the solarwind acceleration region concerns the dominant length scales and time scales of the AWs launched by the Sun. For the correlation lengths and correlation times that we have considered in this work, the two-component analytic model developed in Section 4 reasonably approximates our simulation results, and the Elsasser power spectra in our simulations evolve, at least approximately, towards the scalings of the LGS model at r /greaterorsimilar r A . However, we have also carried out another simulation with higher-frequency photospheric forcing and the same perpendicular correlation length as in Run 3. This additional simulation is not well described by either our two-component model or the LGS model. For example, at very large r , the z + power spectrum evolves towards a k -1 ⊥ scaling, albeit at radii for which δB /similarequal B 0 . We will describe this simulation in more detail in a future publication, but we mention it now to caution the reader that the picture we have developed in this paper does not apply universally for all combinations of correlation times and correlation lengths at the photosphere.

## 8. Phase Mixing

By focusing on transverse, non-compressive fluctuations and neglecting density fluctuations, we neglect 'phase mixing' (Heyvaerts &amp; Priest 1983), by which we mean the process in which an initially planar AW phase front becomes corrugated as it propagates through a medium in which v A (or U ) varies across the magnetic field. This corrugation corresponds to a transfer of fluctuation energy to larger k ⊥ . Phase mixing could provide the additional heating that seems to be needed (see figure 6) to power the fast solar wind at r /lessorsimilar 1 3 . R /circledot over and above the heating provided by reflection-driven MHD turbulence. Observations of comet Lovejoy show that the density varies by a factor of ∼ 6 over distances of a few thousand km perpendicular to B 0 at r = 1 3 . R /circledot in both closedfield regions and open-field regions (Raymond et al. 2014). On the other hand, Helios radio occultation data indicate that the fractional density variations are /similarequal 0 1 . -0 2 . at 5 R /circledot &lt; r &lt; 20 R /circledot (Hollweg et al. 2010). We conjecture that the transition from large δn/n 0 at r /similarequal 1 3 . R /circledot to small δn/n 0 at r /greaterorsimilar 5 R /circledot results from mixing of density fluctuations by the non-compressive component of the turbulence, which acts to reduce δn/n 0 as plasma flows away from the Sun. The limited radial extent of the largeδn/n 0 region suggests that most of the phase mixing occurs close to the Sun. Moreover, since phase mixing is more effective for AWs with larger parallel wavenumbers and frequencies, phase mixing at r /lessorsimilar 5 R /circledot may act as a low-pass filter, by preferentially removing high-frequency AW fluctuation energy. Future investigations of reflection-driven MHD turbulence that account for phase mixing will be important for developing a more complete understanding of solar-wind turbulence and its role in the origin of the solar wind.

## 9. Conclusion

We have carried out three direct numerical simulations of reflection-driven MHD turbulence within a narrow magnetic flux tube that extends from the photosphere, through the chromosphere, through a coronal hole, and out to a maximum heliocentric distance of 21 R /circledot . Our simulations assume fixed, observationally motivated profiles for ρ , U , and B 0 and solve only for the non-compressive, transverse components of the fluctuating magnetic field and velocity. In each simulation, the turbulence is driven by an imposed, randomly evolving, photospheric velocity field that has a single characteristic time scale and length scale. Because outward-propagating AWs undergo strong reflection at the transition region, there is an approximately equal mix of z + and z -fluctuations in the chromosphere, and vigorous turbulence develops within the chromosphere (van Ballegooijen et al. 2011). As a result, the waves that escape into the corona have a broad spectrum of wavenumbers and frequencies. In the corona and solar wind, outwardpropagating z + fluctuations undergo partial non-WKB reflection, thereby generating inward-propagating z -fluctuations, but z + rms /greatermuch z -rms .

In order to explain the radial profiles of z ± rms and the turbulent-heating rate in our simulations, we have developed an analytic model of reflection-driven MHD turbulence that relies on the following conjectures: (i) the Sun launches two populations of z + fluctuations into the corona, a short-radial-correlation-length (HF) population and a longradial-correlation-length (LF) population; (ii) non-WKB reflection of LF z + fluctuations is the dominant source of z -fluctuations; (iii) LF z + fluctuations become aligned with z -at r &gt; r m , where r m is defined in (3.12), causing LF z + fluctuations to cascade and dissipate more slowly than HF z + fluctuations; (iv) the change in sign of d v A / d r at r = r m leads to a reduction in z -rms at r &lt; r m ; and (v) z -fluctuations are anomalously coherent in a reference frame that moves outward with the z + fluctuations, because the

z -fluctuations are produced by the outward-propagating z + fluctuations via non-WKB reflection (Velli et al. 1989; Verdini et al. 2009).

To compare our analytic model and numerical results, we determine the inner boundary conditions in our model by setting the quantities listed in the left column of table 3 equal to their measured or inferred values at the coronal base in our simulations. We then vary the five free parameters in our model (see table 4) to maximize the agreement between the model and simulations, using a single set of free-parameter values to match all three simulations. The resulting best-fit profiles of z ± rms and Q in our model agree reasonably well with our numerical results. The turbulent-heating rate in our simulations is also comparable to the turbulent-heating rate in the solar-wind model of Chandran et al. (2011) at r /greaterorsimilar 1 3 . R /circledot , which agreed with a number of observational constraints. This suggests that MHD turbulence can account for much of the heating that occurs in the fast solar wind.

The inertial-range Elsasser power spectra in our simulations vary with radius. In the lower chromosphere, the spectral indices α + and α -(defined in (3.35)) are /similarequal 3 / 2 , consistent with theories of balanced RMHD turbulence (Section 5.1). In Runs 1 and 2, α + drops with increasing r in the upper chromosphere, reaching values less than 1 just above the transition region. We attribute this spectral flattening to the steep Alfvén-speed gradient in the upper chromosphere, which acts as a high-pass filter (Velli 1993; Réville et al. 2018), as discussed in Section 3.10. Much farther from the Sun, at r /greaterorsimilar 10 R /circledot , α + and α -are reasonably close to 3/2 in all three runs, in approximate agreement with the LGS model of strong imbalanced turbulence, which is reviewed in Section 5.2.2. However, at smaller radii, between r /similarequal r m = 1 7 . R /circledot and r /similarequal 2 r m , α + hovers near unity in Runs 1 and 2. We attribute this latter behavior to a disruption of the anomalous coherence of inertial-range z -fluctuations in the z + reference frame. This disruption is caused by the sign change in d v A / d r at r = r m , which, as shown by van Ballegooijen &amp; Asgari-Targhi (2017), leads to an inverse cascade of z + energy in this same region (Section 5.4).

As mentioned in Section 7, we have carried out additional, as-yet-unpublished, numerical simulations similar to the ones we report here, but with different photospheric boundary conditions. For some values of the correlation length and correlation time of the photospheric velocity field, the fluctuations at r /greaterorsimilar 10 R /circledot conform to neither the analytic model of Section 4 nor the LGS model described in Section 5.2.2. Determining how the properties of non-compressive turbulence at r /greaterorsimilar 10 R /circledot depend upon the photospheric boundary conditions remains an open problem. Further work is also needed to determine how compressive and non-compressive fluctuations interact and evolve as they propagate away from the Sun and also to investigate the role of non-transverse (e.g., spherically polarized) fluctuations (see, e.g., Vasquez &amp; Hollweg 1996; Horbury et al. 2018; Squire et al. 2019).

Observations have led to a detailed picture of solar-wind turbulence at r /similarequal 1 au (e.g., Belcher &amp; Davis 1971; Matthaeus &amp; Goldstein 1982; Bruno &amp; Carbone 2005; Podesta et al. 2007; Horbury et al. 2008; Chen et al. 2012; Wicks et al. 2013a). With the recent launch of NASA's Parker Solar Probe (Fox et al. 2016), it will soon become possible to measure velocity and density fluctuations (Kasper et al. 2016) as well as electric-field and magnetic-field fluctuations (Bale et al. 2016) at heliocentric distances as small as 9 8 . R /circledot . Such measurements will provide critical tests for numerical and theoretical models such as the ones we have presented here.

We thank M. Asgari-Targhi, A. Schekochihin, and A. van Ballegooijen for helpful discussions. We also thank the three reviewers for their comments and suggestions,

which helped improve the manuscript. This work was supported in part by NASA grants NNX11AJ37G, NNX15AI80, NNX16AG81G, NNX16AH92G, NNX17AI18G, and 80NSSC19K0829, NASA grant NNN06AA01C to the Parker Solar Probe FIELDS Experiment, and NSF grant PHY-1500041. High-performance-computing resources were provided by the Argonne Leadership Computing Facility (ALCF) at Argonne National Laboratory, which is supported by the Office of Science of the U.S. Department of Energy under contract DE-AC02-06CH11357. The ALCF resources were granted under INCITE projects from 2012 to 2014. High-performance computing resources were also provided by the Texas Advanced Computing Center under the NSF-XSEDE Project TG-ATM100031.

## REFERENCES

Bale, S. D., Goetz, K., Harvey, P. R., Turin, P., Bonnell, J. W., Dudok de Wit, T., Ergun, R. E., MacDowall, R. J., Pulupa, M., Andre, M., Bolton, M., Bougeret, J.-L., Bowen, T. A., Burgess, D., Cattell, C. A., Chandran, B. D. G., Chaston, C. C., Chen, C. H. K., Choi, M. K., Connerney, J. E., Cranmer, S., Diaz-Aguado, M., Donakowski, W., Drake, J. F., Farrell, W. M., Fergeau, P., Fermin, J., Fischer, J., Fox, N., Glaser, D., Goldstein, M., Gordon, D., Hanson, E., Harris, S. E., Hayes, L. M., Hinze, J. J., Hollweg, J. V., Horbury, T. S., Howard, R. A., Hoxie, V., Jannet, G., Karlsson, M., Kasper, J. C., Kellogg, P. J., Kien, M., Klimchuk, J. A., Krasnoselskikh, V. V., Krucker, S., Lynch, J. J., Maksimovic, M., Malaspina, D. M., Marker, S., Martin, P., Martinez-Oliveros, J., McCauley, J., McComas, D. J., McDonald, T., Meyer-Vernet, N., Moncuquet, M., Monson, S. J., Mozer, F. S., Murphy, S. D., Odom, J., Oliverson, R., Olson, J., Parker, E. N., Pankow, D., Phan, T., Quataert, E., Quinn, T., Ruplin, S. W., Salem, C., Seitz, D., Sheppard, D. A., Siy, A., Stevens, K., Summers, D., Szabo, A., Timofeeva, M., Vaivads, A., Velli, M., Yehle, A., Werthimer, D. &amp; Wygant, J. R. 2016 The FIELDS Instrument Suite for Solar Probe Plus. Measuring the Coronal Plasma and Magnetic Field, Plasma Waves and Turbulence, and Radio Signatures of Solar Transients. Sp. Sci. Rev. 204 , 49-82.

Banerjee, D., Teriaca, L., Doyle, J. G. &amp; Wilhelm, K. 1998 Broadening of SI VIII lines observed in the solar polar coronal holes. Astron. Astrophys. 339 , 208-214.

Barnes, A. 1966 Collisionless damping of hydromagnetic waves. Physics of Fluids 9 , 1483-1495.

Belcher, J. W. &amp; Davis, Jr., L. 1971 Large-amplitude Alfvén waves in the interplanetary medium, 2. J. Geophys. Res. 76 , 3534-3563.

Beresnyak, A. 2012 Basic properties of magnetohydrodynamic turbulence in the inertial range. Mon. Not. R. Astron. Soc. 422 , 3495-3502, arXiv: 1111.5329.

Beresnyak, A. 2014 Spectra of Strong Magnetohydrodynamic Turbulence from High-resolution Simulations. Astrophys. J. Lett. 784 , L20, arXiv: 1401.4177.

Beresnyak, A. &amp; Lazarian, A. 2008 Strong Imbalanced Turbulence. Astrophys. J. 682 , 1070-1075, arXiv: 0709.0554.

Beresnyak, A. &amp; Lazarian, A. 2009 Structure of Stationary Strong Imbalanced Turbulence. Astrophys. J. 702 , 460-471, arXiv: 0904.2574.

Boldyrev, S. 2005 On the Spectrum of Magnetohydrodynamic Turbulence. Astrophys. J. Lett. 626 , L37-L40, arXiv: arXiv:astro-ph/0503053.

Boldyrev, S. 2006 Spectrum of Magnetohydrodynamic Turbulence. Physical Review Letters 96 (11), 115002-+, arXiv: arXiv:astro-ph/0511290.

Boldyrev, S. &amp; Perez, J. C. 2009 Spectrum of Weak Magnetohydrodynamic Turbulence. Physical Review Letters 103 (22), 225001, arXiv: 0907.4475.

Boldyrev, S., Perez, J. C., Borovsky, J. E. &amp; Podesta, J. J. 2011 Spectral Scaling Laws in Magnetohydrodynamic Turbulence Simulations and in the Solar Wind. Astrophys. J. Lett. 741 , L19, arXiv: 1106.0700.

Bruno, R. &amp; Carbone, V. 2005 The Solar Wind as a Turbulence Laboratory. Living Reviews in Solar Physics 2 , 4-+.

Canuto, C., Hussaini, M., Quarteroni, A. &amp; Zang, T. 1988 Spectral Methods in Fluid Dynamics . Springer-Verlag.

Chandran, B. D. G. 2005 Weak Compressible Magnetohydrodynamic Turbulence in the Solar Corona. Phys. Rev. Lett. 95 (26), 265004-+, arXiv: arXiv:astro-ph/0511586.

Chandran, B. D. G. 2008 Weakly Turbulent Magnetohydrodynamic Waves in Compressible Lowβ Plasmas. Phys. Rev. Lett. 101 (23), 235004-+, arXiv: 0810.5360.

Chandran, B. D. G. 2008a Strong Anisotropic MHD Turbulence with Cross Helicity. Astrophys. J. 685 , 646-658, arXiv: 0801.4903.

Chandran, B. D. G. 2018 Parametric instability, inverse cascade and the range of solar-wind turbulence. Journal of Plasma Physics 84 , 905840106.

Chandran, B. D. G., Dennis, T. J., Quataert, E. &amp; Bale, S. D. 2011 Incorporating Kinetic Physics into a Two-fluid Solar-wind Model with Temperature Anisotropy and Low-frequency Alfvén-wave Turbulence. Astrophys. J. 743 , 197, arXiv: 1110.3029.

Chandran, B. D. G. &amp; Hollweg, J. V. 2009 Alfvén Wave Reflection and Turbulent Heating in the Solar Wind from 1 Solar Radius to 1 AU: An Analytical Treatment. Astrophys. J. 707 , 1659-1667, arXiv: 0911.1068.

Chandran, B. D. G., Perez, J. C., Verscharen, D., Klein, K. G. &amp; Mallet, A. 2015 a On the Conservation of Cross Helicity and Wave Action in Solar-wind Models with Non-WKB Alfvén Wave Reflection. Astrophys. J. 811 , 50, arXiv: 1509.01135.

Chandran, B. D. G., Schekochihin, A. A. &amp; Mallet, A. 2015 b Intermittency and Alignment in Strong RMHD Turbulence. Astrophys. J. 807 , 39, arXiv: 1403.6354.

Chen, C. H. K., Mallet, A., Schekochihin, A. A., Horbury, T. S., Wicks, R. T. &amp; Bale, S. D. 2012 Three-dimensional Structure of Solar Wind Turbulence. Astrophys. J. 758 , 120, arXiv: 1109.2558.

Cho, J. &amp; Lazarian, A. 2003 Compressible magnetohydrodynamic turbulence: Mode coupling, scaling relations, anisotropy, viscosity-damped regime and astrophysical implications. Mon. Not. R. Astron. Soc. 345 , 325-339.

Cho, J. &amp; Vishniac, E. T. 2000 The Anisotropy of Magnetohydrodynamic Alfvénic Turbulence. Astrophys. J. 539 , 273-282, arXiv: arXiv:astro-ph/0003403.

Cohen, R. H. &amp; Dewar, R. L. 1974 On the backscatter instability of solar wind Alfven waves. J. Geophys. Res. 79 , 4174-4178.

Cranmer, S. R. &amp; van Ballegooijen, A. A. 2005 On the generation, propagation, and reflection of Alfvén waves from the solar photosphere to the distant heliosphere. Astrophys. J. Suppl. 156 , 265-293.

Cranmer, S. R., van Ballegooijen, A. A. &amp; Edgar, R. J. 2007 Self-consistent Coronal Heating and Solar Wind Acceleration from Anisotropic Magnetohydrodynamic Turbulence. Astrophys. J. Suppl. 171 , 520-551, arXiv: arXiv:astro-ph/0703333.

De Pontieu, B., McIntosh, S. W., Carlsson, M., Hansteen, V. H., Tarbell, T. D., Schrijver, C. J., Title, A. M., Shine, R. A., Tsuneta, S., Katsukawa, Y., Ichimoto, K., Suematsu, Y., Shimizu, T. &amp; Nagata, S. 2007 Chromospheric Alfvénic Waves Strong Enough to Power the Solar Wind. Science 318 , 1574-7.

Dmitruk, P. &amp; Matthaeus, W. H. 2003 Low-Frequency Waves and Turbulence in an Open Magnetic Region: Timescales and Heating Efficiency. Astrophys. J. 597 , 1097-1105.

Dmitruk, P., Matthaeus, W. H., Milano, L. J., Oughton, S., Zank, G. P. &amp; Mullan, D. J. 2002 Coronal heating distribution due to low-frequency, wave-driven turbulence. Astrophys. J. 575 , 571-577.

Esser, R., Fineschi, S., Dobrzycka, D., Habbal, S. R., Edgar, R. J., Raymond, J. C., Kohl, J. L. &amp; Guhathakurta, M. 1999 Plasma properties in coronal holes derived from measurements of minor ion spectral lines and polarized white light intensity. Astrophys. J. Lett. 510 , L63-L67.

Feldman, W. C., Habbal, S. R., Hoogeveen, G. &amp; Wang, Y. 1997 Experimental constraints on pulsed and steady state models of the solar wind near the Sun. J. Geophys. Res. 102 , 26905-26918.

Fox, N. J., Velli, M. C., Bale, S. D., Decker, R., Driesman, A., Howard, R. A., Kasper, J. C., Kinnison, J., Kusterer, M., Lario, D., Lockwood, M. K., McComas, D. J., Raouafi, N. E. &amp; Szabo, A. 2016 The Solar Probe Plus Mission: Humanity's First Visit to Our Star. Sp. Sci. Rev. 204 , 7-48.

Galeev, A. A. &amp; Oraevskii, V. N. 1963 The Stability of Alfvén Waves. Soviet Physics Doklady 7 , 988.

Galtier, S., Nazarenko, S. V., Newell, A. C. &amp; Pouquet, A. 2000 A weak turbulence theory for incompressible magnetohydrodynamics. Journal of Plasma Physics 63 , 447488.

Goldreich, P. &amp; Sridhar, S. 1995 Toward a theory of interstellar turbulence. 2: Strong alfvenic turbulence. Astrophys. J. 438 , 763-775.

Goldreich, P. &amp; Sridhar, S. 1997 Magnetohydrodynamic turbulence revisited. Astrophys. J. 485 , 680-688.

Heinemann, M. &amp; Olbert, S. 1980 Non-WKB Alfven waves in the solar wind. J. Geophys. Res. 85 , 1311-1327.

Heyvaerts, J. &amp; Priest, E. R. 1983 Coronal heating by phase-mixed shear Alfven waves. Astron. Astrophys. 117 , 220-234.

Hollweg, J. V., Cranmer, S. R. &amp; Chandran, B. D. G. 2010 Coronal Faraday Rotation Fluctuations and a Wave/Turbulence-driven Model of the Solar Wind. Astrophys. J. 722 , 1495-1503.

Hollweg, J. V. &amp; Isenberg, P. A. 2002 Generation of the fast solar wind: A review with emphasis on the resonant cyclotron interaction. Journal of Geophysical Research (Space Physics) 107 , 1147-+.

Hollweg, J. V. &amp; Isenberg, P. A. 2007 Reflection of Alfvén waves in the corona and solar wind: An impulse function approach. Journal of Geophysical Research (Space Physics) 112 , 8102-+.

Horbury, T. S., Forman, M. &amp; Oughton, S. 2008 Anisotropic Scaling of Magnetohydrodynamic Turbulence. Physical Review Letters 101 (17), 175005-+, arXiv: 0807.3713.

Horbury, T. S., Matteini, L. &amp; Stansby, D. 2018 Short, large-amplitude speed enhancements in the near-Sun fast solar wind. Mon. Not. R. Astron. Soc. 478 , 19801986.

Kasper, J. C., Abiad, R., Austin, G., Balat-Pichelin, M., Bale, S. D., Belcher, J. W., Berg, P., Bergner, H., Berthomier, M., Bookbinder, J., Brodu, E., Caldwell, D., Case, A. W., Chandran, B. D. G., Cheimets, P., Cirtain, J. W., Cranmer, S. R., Curtis, D. W., Daigneau, P., Dalton, G., Dasgupta, B., DeTomaso, D., Diaz-Aguado, M., Djordjevic, B., Donaskowski, B., Effinger, M., Florinski, V., Fox, N., Freeman, M., Gallagher, D., Gary, S. P., Gauron, T., Gates, R., Goldstein, M., Golub, L., Gordon, D. A., Gurnee, R., Guth, G., Halekas, J., Hatch, K., Heerikuisen, J., Ho, G., Hu, Q., Johnson, G., Jordan, S. P., Korreck, K. E., Larson, D., Lazarus, A. J., Li, G., Livi, R., Ludlam, M., Maksimovic, M., McFadden, J. P., Marchant, W., Maruca, B. A., McComas, D. J., Messina, L., Mercer, T., Park, S., Peddie, A. M., Pogorelov, N., Reinhart, M. J., Richardson, J. D., Robinson, M., Rosen, I., Skoug, R. M., Slagle, A., Steinberg, J. T., Stevens, M. L., Szabo, A., Taylor, E. R., Tiu, C., Turin, P., Velli, M., Webb, G., Whittlesey, P., Wright, K., Wu, S. T. &amp; Zank, G. 2016 Solar Wind Electrons Alphas and Protons (SWEAP) Investigation: Design of the Solar Wind and Coronal Plasma Instrument Suite for Solar Probe Plus. Sp. Sci. Rev. 204 , 131-186.

Lithwick, Y. &amp; Goldreich, P. 2003 Imbalanced weak magnetohydrodynamic turbulence. Astrophys. J. 582 , 1220-1240, arXiv: astro-ph/0208046.

Lithwick, Y., Goldreich, P. &amp; Sridhar, S. 2007 Imbalanced strong MHD turbulence. Astrophys. J. 655 , 269-274, arXiv: astro-ph/0607243.

Loureiro, N. F. &amp; Boldyrev, S. 2017 a Collisionless Reconnection in Magnetohydrodynamic and Kinetic Turbulence. Astrophys. J. 850 , 182.

Loureiro, N. F. &amp; Boldyrev, S. 2017 b Role of Magnetic Reconnection in Magnetohydrodynamic Turbulence. Physical Review Letters 118 (24), 245101.

Luo, Q. &amp; Melrose, D. 2006 Anisotropic weak turbulence of Alfvén waves in collisionless astrophysical plasmas. Mon. Not. R. Astron. Soc. 368 , 1151-1158, arXiv: arXiv:astroph/0602295.

Mallet, A. &amp; Schekochihin, A. A. 2017 A statistical model of three-dimensional anisotropy

and intermittency in strong Alfvénic turbulence. Mon. Not. R. Astron. Soc. 466 , 39183927, arXiv: 1606.00466.

Mallet, A., Schekochihin, A. A. &amp; Chandran, B. D. G. 2015 Refined critical balance in strong Alfvénic turbulence. Mon. Not. R. Astron. Soc. 449 , L77-L81, arXiv: 1406.5658.

Mallet, A., Schekochihin, A. A. &amp; Chandran, B. D. G. 2017 a Disruption of Alfvénic turbulence by magnetic reconnection in a collisionless plasma. Journal of Plasma Physics 83 , 905830609.

Mallet, A., Schekochihin, A. A. &amp; Chandran, B. D. G. 2017 b Disruption of sheet-like structures in Alfvénic turbulence by magnetic reconnection. Mon. Not. R. Astron. Soc. 468 , 4862-4871, arXiv: 1612.07604.

Maron, J. &amp; Goldreich, P. 2001 Simulations of Incompressible Magnetohydrodynamic Turbulence. Astrophys. J. 554 , 1175-1196, arXiv: arXiv:astro-ph/0012491.

Mason, J., Cattaneo, F. &amp; Boldyrev, S. 2008 Numerical measurements of the spectrum in magnetohydrodynamic turbulence. Phys. Rev. E. 77 (3), 036403-+, arXiv: 0706.2003.

Matthaeus, W. H. &amp; Goldstein, M. L.

1982 Measurement of the rugged invariants of magnetohydrodynamic turbulence in the solar wind.

J. Geophys. Res.

87

,

6011-6028.

Meyrand, R., Kanekar, A., Dorland, W. &amp; Schekochihin, A. A. 2019 Fluidization of collisionless plasma turbulence. Proceedings of the National Academy of Science 116 , 1185-1194, arXiv: 1808.04284.

Meyrand, R., Kiyani, K. H. &amp; Galtier, S. 2015 Weak magnetohydrodynamic turbulence and intermittency. Journal of Fluid Mechanics 770 , R1.

Müller, W. &amp; Grappin, R. 2005 Spectral Energy Dynamics in Magnetohydrodynamic Turbulence. Physical Review Letters 95 (11), 114502-+, arXiv: arXiv:physics/0509019.

Ng, C. S. &amp; Bhattacharjee, A. 1996 Interaction of Shear-Alfven Wave Packets: Implication for Weak Magnetohydrodynamic Turbulence in Astrophysical Plasmas. Astrophys. J. 465 , 845-+.

Ng, C. S. &amp; Bhattacharjee, A. 1997 Scaling of anisotropic spectra due to the weak interaction of shear-Alfvén wave packets. Physics of Plasmas 4 , 605-610.

Perez, J. C. &amp; Boldyrev, S. 2008 On Weak and Strong Magnetohydrodynamic Turbulence. Astrophys. J. Lett. 672 , L61-L64, arXiv: arXiv:0712.2086.

Perez, J. C. &amp; Boldyrev, S. 2009 Role of Cross-Helicity in Magnetohydrodynamic Turbulence. Physical Review Letters 102 (2), 025003-+, arXiv: 0807.2635.

Perez, J. C. &amp; Boldyrev, S. 2010 Strong magnetohydrodynamic turbulence with cross helicity. Physics of Plasmas 17 (5), 055903, arXiv: 1004.3798.

Perez, J. C. &amp; Chandran, B. D. G. 2013 Direct Numerical Simulations of Reflection-Driven, Reduced MHD Turbulence from the Sun to the Alfven Critical Point. Astrophys. J. 776 , 124, arXiv: 1308.4046.

Perez, J. C., Mason, J., Boldyrev, S. &amp; Cattaneo, F. 2012 On the Energy Spectrum of Strong Magnetohydrodynamic Turbulence. Physical Review X 2 (4), 041005, arXiv: 1209.2011.

Podesta, J. J. &amp; Bhattacharjee, A. 2010 Theory of Incompressible Magnetohydrodynamic Turbulence with Scale-dependent Alignment and Cross-helicity. Astrophys. J. 718 , 11511157, arXiv: 0903.5041.

Podesta, J. J., Roberts, D. A. &amp; Goldstein, M. L. 2007 Spectral Exponents of Kinetic and Magnetic Energy Spectra in Solar Wind Turbulence. Astrophys. J. 664 , 543-548.

Pucci, F. &amp; Velli, M. 2014 Reconnection of Quasi-singular Current Sheets: The 'Ideal' Tearing Mode. Astrophys. J. Lett. 780 , L19.

Pucci, F., Velli, M., Tenerani, A. &amp; Del Sarto, D. 2018 Onset of fast 'ideal' tearing in thin current sheets: Dependence on the equilibrium current profile. Physics of Plasmas 25 (3), 032113, arXiv: 1801.08412.

Raymond, J. C., McCauley, P. I., Cranmer, S. R. &amp; Downs, C. 2014 The Solar Corona as Probed by Comet Lovejoy (C/2011 W3). Astrophys. J. 788 , 152, arXiv: 1405.1639.

Réville, V., Tenerani, A. &amp; Velli, M. 2018 Parametric Decay and the Origin of the Lowfrequency Alfvénic Spectrum of the Solar Wind. Astrophys. J. 866 , 38.

Richardson, R. S. &amp; Schwarzschild, M. 1950 On the Turbulent Velocities of Solar Granules. Astrophys. J. 111 , 351.

Sagdeev, R. Z. &amp; Galeev, A. A. 1969 Nonlinear Plasma Theory .

Schekochihin, A. A., Cowley, S. C., Dorland, W., Hammett, G. W., Howes, G. G., Quataert, E. &amp; Tatsuno, T. 2009 Astrophysical Gyrokinetics: Kinetic and Fluid Turbulent Cascades in Magnetized Weakly Collisional Plasmas. Astrophys. J. Suppl. 182 , 310-377, arXiv: 0704.0044.

Schekochihin, A. A., Nazarenko, S. V. &amp; Yousef, T. A. 2012 Weak Alfvén-wave turbulence revisited. Phys. Rev. E. 85 (3), 036406, arXiv: 1110.6682.

Schekochihin, A. A., Parker, J. T., Highcock, E. G., Dellar, P. J., Dorland, W. &amp; Hammett, G. W. 2016 Phase mixing versus nonlinear advection in drift-kinetic plasma turbulence. Journal of Plasma Physics 82 (2), 905820212, arXiv: 1508.05988.

Shebalin, J. V., Matthaeus, W. &amp; Montgomery, D. 1983 Anisotropy in MHD turbulence due to a mean magnetic field. Journal of Plasma Physics 29 , 525.

Shoda, M., Suzuki, T. K., Asgari-Targhi, M. &amp; Yokoyama, T. 2019 Three-dimensional Simulation of the Fast Solar Wind Driven by Compressible Magnetohydrodynamic Turbulence. The Astrophysical Journal 880 (1), L2, arXiv: 1905.11685.

Squire, J., Schekochihin, A. A., Quataert, E. &amp; Kunz, M. W. 2019 Magneto-immutable turbulence in weakly collisional plasmas. Journal of Plasma Physics 85 (1), 905850114, arXiv: 1811.12421.

Tenerani, A., Velli, M. &amp; Hellinger, P. 2017 The Parametric Instability of Alfvén Waves: Effects of Temperature Anisotropy. Astrophys. J. 851 , 99.

Tu, C. &amp; Marsch, E. 1995 MHD structures, waves and turbulence in the solar wind: Observations and theories. Space Science Reviews 73 , 1-210.

Usmanov, A. V., Goldstein, M. L. &amp; Matthaeus, W. H. 2014 Three-fluid, Threedimensional Magnetohydrodynamic Solar Wind Model with Eddy Viscosity and Turbulent Resistivity. Astrophys. J. 788 , 43.

van Ballegooijen, A. A. &amp; Asgari-Targhi, M. 2016 Heating and Acceleration of the Fast Solar Wind by Alfvén Wave Turbulence. Astrophys. J. 821 , 106, arXiv: 1602.06883.

van Ballegooijen, A. A. &amp; Asgari-Targhi, M. 2017 Direct and Inverse Cascades in the Acceleration Region of the Fast Solar Wind. Astrophys. J. 835 , 10, arXiv: 1612.02501.

van Ballegooijen, A. A., Asgari-Targhi, M., Cranmer, S. R. &amp; DeLuca, E. E. 2011 Heating of the Solar Chromosphere and Corona by Alfvén Wave Turbulence. Astrophys. J. 736 , 3-+, arXiv: 1105.0402.

van der Holst, B., Sokolov, I. V., Meng, X., Jin, M., Manchester, IV, W. B., Tóth, G. &amp; Gombosi, T. I. 2014 Alfvén Wave Solar Model (AWSoM): Coronal Heating. Astrophys. J. 782 , 81, arXiv: 1311.4093.

Vasquez, B. J. &amp; Hollweg, J. V. 1996 Formation of arc-shaped Alfvén waves and rotational discontinuities from oblique linearly polarized wave trains. J. Geophys. Res. 101 , 1352713540.

Vech, D., Mallet, A., Klein, K. G. &amp; Kasper, J. C. 2018 Magnetic Reconnection May Control the Ion-scale Spectral Break of Solar Wind Turbulence. Astrophys. J. Lett. 855 , L27, arXiv: 1803.00065.

Velli, M. 1993 On the propagation of ideal, linear Alfven waves in radially stratified stellar atmospheres and winds. Astron. Astrophys. 270 , 304-314.

Velli, M., Grappin, R. &amp; Mangeney, A. 1989 Turbulent cascade of incompressible unidirectional Alfven waves in the interplanetary medium. Physical Review Letters 63 , 1807-1810.

Verdini, A., Grappin, R., Pinto, R. &amp; Velli, M. 2012 On the Origin of the 1/f Spectrum in the Solar Wind Magnetic Field. Astrophys. J. Lett. 750 , L33, arXiv: 1203.6219.

Verdini, A. &amp; Velli, M. 2007 Alfvén waves and turbulence in the solar atmosphere and solar wind. Astrophys. J. 662 , 669-676, arXiv: arXiv:astro-ph/0702205.

Verdini, A., Velli, M. &amp; Buchlin, E. 2009 Turbulence in the Sub-Alfvénic Solar Wind Driven by Reflection of Low-Frequency Alfvén Waves. Astrophys. J. Lett. 700 , L39-L42, arXiv: 0905.2618.

Verdini, A., Velli, M., Matthaeus, W. H., Oughton, S. &amp; Dmitruk, P. 2010 A Turbulence-Driven Model for Heating and Acceleration of the Fast Wind in Coronal Holes. Astrophys. J. Lett. 708 , L116-L120, arXiv: 0911.5221.

Wicks, R. T., Mallet, A., Horbury, T. S., Chen, C. H. K., Schekochihin, A. A. &amp;

Mitchell, J. J. 2013a Alignment and Scaling of Large-Scale Fluctuations in the Solar Wind. Physical Review Letters 110 (2), 025003, arXiv: 1209.5362.

Yoon, P. H. &amp; Fang, T.-M. 2009 Proton heating by parallel Alfvén wave cascade. Physics of Plasmas 16 (6), 062314.

Zank, G. P., Adhikari, L., Hunana, P., Tiwari, S. K., Moore, R., Shiota, D., Bruno, R. &amp; Telloni, D. 2018 Theory and Transport of Nearly Incompressible Magnetohydrodynamic Turbulence. IV. Solar Coronal Turbulence. Astrophys. J. 854 , 32.

Zhou, Y. &amp; Matthaeus, W. H. 1989 Non-WKB evolution of solar wind fluctuations: turbulence modeling approach. Geophys. Res. Lett. 16 , 755-758.

A

Zhou, Y. &amp; Matthaeus, W. H. 1990 Transport and turbulence modeling of solar wind fluctuations. J. Geophys. Res. 95 , 10291-10311.
Benjamin D. G. Chandran Jean C. Perez
Reflection-driven MHD turbulence in the solar atmosphere and solar wind
We present 3D numerical simulations and an analytic model of
reflection-driven MHD turbulence in the solar wind. Our simulations describe
transverse, non-compressive MHD fluctuations within a narrow magnetic flux tube
that extends from the photosphere out to a heliocentric distance $r$ of 21
solar radii $(R_s)$. We launch outward-propagating "$z^+$ fluctuations" into
the simulation domain by imposing a randomly evolving photospheric velocity
field. As these fluctuations propagate away from the Sun, they undergo partial
reflection, producing inward-propagating "$z^-$ fluctuations."
Counter-propagating fluctuations subsequently interact, causing fluctuation
energy to cascade to small scales and dissipate. Our analytic model
incorporates alignment, allows for strongly or weakly turbulent nonlinear
interactions, and divides the $z^+$ fluctuations into two populations with
different characteristic radial correlation lengths. The inertial-range power
spectra in our simulations evolve toward a $k_\perp^{-3/2}$ scaling at $r>10
R_s$, where $k_\perp$ is the wave-vector component perpendicular to the
background magnetic field. In two of our simulations, the $z^+$ power spectra
are much flatter between the coronal base and $r \simeq 4 R_s$. We argue that
these spectral scalings are caused by: (1) high-pass filtering in the upper
chromosphere; (2) the anomalous coherence of inertial-range $z^-$ fluctuations
in a reference frame propagating outwards with the $z^+$ fluctuations; and (3)
the change in the sign of the radial derivative of the Alfv\'en speed at $r=r_m
\simeq 1.7 R_s$, which disrupts this anomalous coherence between $r=r_m$ and
$r\simeq 2r_m$. At $r>1.3 R_s$, the turbulent heating rate in our simulations
is comparable to the heating rate in a previously developed solar-wind model
that agreed with a number of observational constraints.