2.1. Magnetic geometry and numerical domain

We limit our model to heliocentric distances less thanFootnote ³

(2.1) \begin{equation} r_{\textrm {max}} = 72 R_{\odot }. \end{equation}

This enables us to neglect the effects of solar rotation on the direction of the background magnetic field $\boldsymbol{B}$ and to treat $\boldsymbol{B}$ as approximately radial. We then analyse the solar-wind outflow within a narrow magnetic flux tube centred on a radial magnetic-field line. We take the solar-wind outflow velocity $\boldsymbol{U}$ to be everywhere parallel to the magnetic field $\boldsymbol{B}$ , setting

(2.2) \begin{equation} \boldsymbol{U} = U \boldsymbol{\hat {b}}, \end{equation}

where $\boldsymbol{\hat {b}} = \boldsymbol{B}/B$ . The magnetic-field-strength profile $B(r)$ is arbitrary, but fixed in time. Although we allow the magnetic field to undergo super-radial expansion (Kopp & Holzer Reference Kopp and Holzer1976), the radially oriented magnetic flux tube that we consider is sufficiently narrow that we may approximate

(2.3) \begin{equation} \boldsymbol{B} \boldsymbol{\cdot} \boldsymbol{\nabla } \rightarrow B \frac {\partial }{\partial r}. \end{equation}

The cross-sectional area of the flow $a(r)$ is related to $B(r)$ by magnetic flux conservation,

(2.4) \begin{equation} \frac {a(r)}{a(r_{\textrm {b}})}\, = \, \frac {B(r_{\textrm {b}})}{B(r)}, \end{equation}

where $r_{\textrm {b}}$ is the minimum heliocentric distance in the model, which we take to lie somewhere in the upper transition region, just beneath the corona.

2.2. Moments of the Vlasov equation

Following a number of authors (e.g. Kulsrud Reference Kulsrud1983; Snyder, Hammett & Dorland Reference Snyder, Hammett and Dorland1997; Sharma et al. Reference Sharma, Hammett, Quataert and Stone2006), we take moments of the Vlasov equation for a proton--electron plasma in the limit in which the plasma frequency and cyclotron frequency are much larger than any other frequency of interest. In this limit, the plasma is quasi-neutral, and the proton and electron distribution functions are gyrotropic. We view all moments of these distribution functions (e.g. the proton density $n$ , which equals the electron density) as the sum of a background value and a fluctuation, where the background value of a quantity is its average over the cross-section of the flux tube at fixed $r$ . The background values are thus functions of $r$ and $t$ alone.

The resulting equations allow for proton temperature anisotropy, but we assume an isotropic electron temperature $T_{\textrm {e}}$ and determine the radial component of the electron heat flux $q_{\textrm {e}}$ using a simple closure (§2.4). In contrast, $q_{\perp \textrm {p}}$ and $q_{\parallel \textrm {p}}$ , the radial components of the perpendicular and parallel proton heat fluxes, respectively, evolve dynamically in a way that depends upon three different types of fourth velocity moments of the proton distribution function. To close the proton fluid equations, we take these fourth velocity moments to have the same values that they would have in a bi-Maxwellian plasma with the same $T_{\perp \textrm {p}}$ and $T_{\parallel \textrm {p}}$ as the actual plasma. We incorporate reflection-driven AW turbulence by treating the dominant, outward-propagating AWs as a separate fluid with energy density $\mathcal {E}_{\textrm {w}}$ that interacts with the plasma via a wave-pressure force and turbulent heating.

The eight dependent variables in our model ( $n$ , $U$ , $T_{\textrm {e}}$ , $T_{\perp \textrm {p}}$ , $T_{\parallel \textrm {p}}$ , $q_{\perp \mathrm{p}}$ , $q_{\parallel \textrm {p}}$ and $\mathcal {E}_{\textrm {w}}$ ) satisfy the following eight equations (Chandran et al. Reference Chandran, Dennis, Quataert and Bale2011):

(2.5) \begin{equation} \frac {\mathrm{d} n}{\mathrm{d} t}\, =\, -\, \frac {n}{a} \frac {\partial }{\partial r} \left ( a U\right )\!, \end{equation}

(2.6) \begin{equation} \frac {\mathrm{d} U}{\mathrm{d}t} \, = \, -\, \frac {k_{\mathrm{B}}}{\rho } \frac {\partial }{\partial r}\left [ n(T_{\mathrm{e}} + T_{\parallel \textrm {p}})\right ] + \frac {k_{\textrm {B}}(T_{\perp \textrm {p}} - T_{\parallel \textrm {p}})}{m_{\mathrm{p}} a} \frac {\partial a}{\partial r} - \frac {G M_{\odot }}{r^2} - \frac {1}{2\rho } \frac {\partial {\mathcal E}_{\mathrm{w}}}{\partial r}, \end{equation}

(2.7) \begin{equation} \frac {3}{2}n^{5/3} k_{\mathrm{B}} \frac {\mathrm{d}}{\mathrm{d}t}\left (\frac {T_{\mathrm{e}}}{n^{2/3}}\right ) \, =\, Q_{\mathrm{e}} - \frac {1}{a} \frac {\partial }{\partial r} (a q_{\mathrm{e}}) + 3\nu _{\mathrm{pe}} n k_{\mathrm{B}} (T_{\mathrm{p}} - T_{\mathrm{e}}) - n^2 \varLambda _{\mathrm{rad}}(T_{\mathrm{e}}), \end{equation}

(2.8) \begin{align} B n k_{\mathrm{B}} \frac {\mathrm{d}}{\mathrm{d}t}\left ( \frac {T_{\perp \textrm {p}}}{B}\right ) &= Q_{\perp \textrm {p}} - \frac {1}{a^2}\frac {\partial }{\partial r}\left (a^2 q_{\perp \textrm {p}}\right )\nonumber\\&\quad + \frac {1}{3} \nu _{\mathrm{p}} n k_{\mathrm{B}}(T_{\parallel \textrm {p}} - T_{\perp \textrm {p}}) + 2\nu _{\mathrm{pe}} n k_{\mathrm{B}} (T_{\mathrm{e}} - T_{\perp \textrm {p}}), \end{align}

(2.9) \begin{align} \frac {n^3 k_{\mathrm{B}}}{2B^2} \frac {\mathrm{d}}{\mathrm{d}t}\left (\frac {B^2 T_{\parallel \textrm {p}}}{n^2}\right ) &= Q_{\parallel \textrm {p}} - \frac {1}{a} \frac {\partial }{\partial r} ( a q_{\parallel \textrm {p}})\nonumber\\&\quad + \frac {q_{\perp \textrm {p}}}{a}\frac {\partial a}{\partial r} + \frac {1}{3} \nu _{\mathrm{p}} n k_{\mathrm{B}}(T_{\perp \textrm {p}} - T_{\parallel \textrm {p}}) + \nu _{\mathrm{pe}} n k_{\mathrm{B}}(T_{\mathrm{e}} - T_{\parallel \textrm {p}}),\end{align}

(2.10) \begin{equation} n^2 \frac {\mathrm{d}}{\mathrm{d}t}\left (\frac {q_{\perp \textrm {p}}}{n^2}\right ) \,=\, - \, \frac {n k_{\mathrm{B}}^2 T_{\parallel \textrm {p}}}{m_{\mathrm{p}}} \frac {\partial T_{\perp \textrm {p}}}{\partial r} + \frac {n k_{\mathrm{B}}^2 T_{\perp \textrm {p}}(T_{\perp \textrm {p}} - T_{\parallel \textrm p})}{m_{\mathrm{p}} a} \frac {\partial a}{\partial r} - \, \nu _{\mathrm{p}} q_{\perp \textrm {p}}, \end{equation}

(2.11) \begin{equation} \frac {n^4}{B^3} \frac {\mathrm{d}}{\mathrm{d}t}\left (\frac {B^3 q_{\parallel \textrm {p}}}{n^4}\right )\, =\, - \frac {3 n k_{\mathrm{B}}^2 T_{\parallel \textrm {p}}}{2m_{\mathrm{p}}} \frac {\partial T_{\parallel \textrm {p}}}{\partial r} - \nu _{\mathrm{p}} q_{\parallel \textrm {p}} \end{equation}

and (Dewar Reference Dewar1970)

(2.12) \begin{equation} \frac {\partial {\mathcal E}_{\mathrm{w}} }{\partial t} + \frac {1}{a} \frac {\partial }{\partial r}\left [a (U + v_{\mathrm{A}}){\mathcal E}_{\mathrm{w}}\right ] + \frac {{\mathcal E}_{\mathrm{w}}}{2a}\frac {\partial }{\partial r}(a U) = - Q . \end{equation}

Here $v_{\mathrm{A}} = B/\sqrt {4 \mathrm{\pi } \rho }$ , $\rho = m_{\mathrm{p}} n$ is the mass density, $M_{\odot }$ is the mass of the Sun and

(2.13) \begin{equation} \frac {\mathrm{d}}{\mathrm{d}t} = \frac {\partial }{\partial t} + U \frac {\partial }{\partial r} \end{equation}

is the time derivative in the plasma frame. The quantities $Q_{\mathrm{e}}$ , $Q_{\perp \textrm {p}}$ and $Q_{\parallel \textrm {p}}$ are, respectively, the electron heating rate, the perpendicular proton heating rate and the parallel proton heating rate resulting from the dissipation of AW turbulence (see §3), and

(2.14) \begin{equation} Q = Q_{\mathrm{e}} + Q_{\perp \textrm {p}} + Q_{\parallel \textrm {p}} \end{equation}

is the total turbulent heating rate. The quantity

(2.15) \begin{equation} \nu _{\mathrm{pe}} = \frac { 4 \sqrt {2 \mathrm{\pi } m_{\mathrm{e}}} \,e^4 n \ln \varLambda }{3 m_{\mathrm{p}} (k_{\mathrm{B}}T_{\mathrm{e}})^{3/2}} \end{equation}

is the Coulomb collision frequency for the exchange of energy between protons and electrons (Schunk Reference Schunk1975), where $m_{\mathrm{p}}$ and $m_{\mathrm{e}}$ are the proton and electron masses, respectively, and $\ln \varLambda$ is the Coulomb logarithm, which we set equal to 23. The proton--proton collision frequency in our model is given by

(2.16) \begin{equation} \nu _{\mathrm{p}} = \nu _{\mathrm{pp, C}} + \nu _{\mathrm{inst}}, \end{equation}

where

(2.17) \begin{equation} \nu _{\mathrm{pp, C}} = \frac { 4 \sqrt {\mathrm{\pi }} \, e^4 n \ln \varLambda }{3 \sqrt {m_{\mathrm{p}}} (k_{\mathrm{B}}T_{\mathrm{p}})^{3/2}} \end{equation}

is the proton--proton Coulomb collision frequency (Schunk Reference Schunk1975),

(2.18) \begin{equation} T_{\mathrm{p}} = \frac {2T_{\perp \textrm {p}} + T_{\parallel \textrm {p}}}{3}, \end{equation}

and $\nu _{\mathrm{inst}}$ is a model temperature-isotropisation rate resulting from small-scale plasma fluctuations that grow when the proton temperature anisotropy exceeds the threshold of the oblique firehose or mirror instability (see §2.5).

The quantity $\varLambda _{\mathrm{rad}}(T_{\mathrm{e}})$ on the right-hand side of (2.7) is the optically thin radiative loss function, which we take to have the same value as the quantity $\varLambda (T)$ in the single-fluid model of Cranmer et al. (Reference Cranmer, van Ballegooijen and Edgar2007). To determine this latter function, we read off the value of $\varLambda (T)$ from figure 1 of Cranmer et al. (Reference Cranmer, van Ballegooijen and Edgar2007) at 35 evenly spaced values of $\log _{10}[T/( \mathrm{1\;K})]$ ranging from $3.6$ to $7.0$ . We then approximate $\log _{10}\varLambda (T)$ as a continuous and piecewise linear function of $\log _{10}T$ that passes through these 35 ordered pairs of $(T, \varLambda (T))$ values. We plot the resulting value of $\varLambda _{\mathrm{rad}}(T_{\mathrm{e}})$ in figure 1.

Figure 1. The optically thin radiative loss function appearing in (2.7) and (2.19).

2.3. Total energy equation

Multiplying (2.6) by $\rho U$ and adding the resulting equation to the sum of (2.7)–(2.9) and (2.12) yields

(2.19) \begin{equation} \frac {\partial {\mathcal E}_{\mathrm{tot}} }{\partial t} + \frac {1}{a} \frac {\partial }{\partial r}(a F_{\mathrm{tot}}) = - n^2 \varLambda _{\mathrm{rad}}(T_{\mathrm{e}}), \end{equation}

where

(2.20) \begin{equation} {\mathcal E}_{\mathrm{tot}} = \frac {\rho U^2}{2} - \frac {G M_{\odot }\rho }{r} + n k_{\mathrm{B}} \left ( \frac {3T_{\mathrm{e}}}{2} + T_{\perp \textrm {p}} + \frac {T_{\parallel \textrm {p}}}{2} \right ) + {\mathcal E}_{\mathrm{w}} \end{equation}

is the total energy density and

(2.21) \begin{align} F_{\mathrm{tot}} & = \frac { \rho U^3}{2} - \frac {U G M_{\odot }\rho }{r} + U n k_{\mathrm{B}} \left ( \frac {5 T_{\mathrm{e}}}{2} + T_{\perp \textrm {p}} + \frac { 3T_{\parallel \textrm {p}}}{2}\right ),\nonumber\\ &\quad + \, q_{\mathrm{e}} + q_{\perp \textrm {p}} + q_{\parallel \textrm {p}} + \left ( \frac {3 U }{2} + v_{\mathrm{A}} \right ){\mathcal E}_{\mathrm{w}} \end{align}

is the total energy flux. In steady state, and at heights $\gtrsim 0.1 R_{\odot }$ above the photosphere where radiative cooling is negligible, $aF_{\mathrm{tot}}$ becomes independent of $r$ .

2.4. Electron heat flux

In the low corona, the Coulomb mean free path is much smaller than the temperature scale height, and the radial component of the electron heat flux is well approximated by the Spitzer formula (Spitzer & Härm Reference Spitzer and Härm1953)

(2.22) \begin{equation} q_{\mathrm{e, S}} = - \alpha _{\mathrm{S}} T_{\mathrm{e}}^{5/2} \frac {\partial T_{\mathrm{e}}}{\partial r}, \end{equation}

where

(2.23) \begin{equation} \alpha _{\mathrm{S}} = \frac {1.84 \times 10^{-5} }{\ln \varLambda } \mbox{ erg} \mbox{ s}^{-1} \mbox{ K}^{-7/2} \mbox{ cm}^{-1}. \end{equation}

We approximate the radial component of the electron heat flux far outside the corona using Hollweg’s collisionless heat-flux formula (Hollweg Reference Hollweg1974, Reference Hollweg1976)

(2.24) \begin{equation} q_{\mathrm{e, H}} = \frac {3}{2} \alpha _{\mathrm{H}} U n k_{\mathrm{B}} T_{\mathrm{e}}, \end{equation}

where $\alpha _{\mathrm{H}}$ is a dimensionless constant that we treat as a free parameter. To obtain a continuous transition between these two heat-flux regimes, we set

(2.25) \begin{equation} q_{\mathrm{e}} = \frac {\left (r/r_{\mathrm{H}}\right )^4 q_{\mathrm{e, H}}}{1 + \left (r/r_{\mathrm{H}}\right )^4} \,+\,\frac {q_{\mathrm{e, S}}}{1 + \left (r/r_{\mathrm{H}}\right )^4}, \end{equation}

where $r_{\mathrm{H}}$ is a free parameter.

2.5. Temperature isotropisation from firehose and mirror instabilities

The proton temperature-anisotropy ratio $R = T_{\perp \textrm {p}}/T_{\parallel \textrm {p}}$ in the solar wind is observed to be limited from above by the mirror instability threshold

(2.26) \begin{equation} R_{\mathrm{m}} = 1 + 0.77(\beta _{\parallel \textrm {p}} + 0.016)^{-0.76} \end{equation}

and from below by the oblique firehose instability threshold

(2.27) \begin{equation} R_{\mathrm{f}} = 1 - 1.4(\beta _{\parallel \textrm {p}} + 0.11)^{-1} \end{equation}

(Kasper, Lazarus & Gary Reference Kasper, Lazarus and Gary2002; Hellinger et al. Reference Hellinger, Trávníček, Kasper and Lazarus2006; Bale et al. Reference Bale, Kasper, Howes, Quataert, Salem and Sundkvist2009),Footnote ⁴ where

(2.28) \begin{equation} \beta _{\parallel \textrm {p}} = \frac {8\mathrm{\pi } n k_{\mathrm{B}} T_{\parallel \textrm p}}{B^2}. \end{equation}

This suggests that when the plasma crosses one of these thresholds, short-wavelength fluctuations grow and cause pitch-angle scattering that drives the plasma back into the stable region of parameter space. We incorporate such instability-induced pitch-angle scattering into our model through the term $\nu _{\mathrm{inst}}$ in (2.16), where

(2.29) \begin{equation} \nu _{\mathrm{inst}} = \nu _0 \exp \left [ \frac {12(R - R_{\mathrm{m}})}{R_{\mathrm{m}}}\right ] + \nu _0 \exp \left [ \frac {12(\overline { R}_{\mathrm{f}} - R)}{ \overline {R}_{\mathrm{f}}}\right ], \end{equation}

$\nu _0 = 0.02\sqrt {G M_{\odot }/R_{\odot }^3}$ and $\overline { R}_{\mathrm{f}} = \max ( R_{\mathrm{f}}, 10^{-6})$ (cf. Sharma et al. Reference Sharma, Hammett, Quataert and Stone2006; Chandran et al. Reference Chandran, Dennis, Quataert and Bale2011).

3.1. The root-mean-square amplitudes of the turbulent fluctuations

We describe the turbulent fluctuations using the Elsasser variables (Elsasser Reference Elsasser1950)

(3.2) \begin{equation} \delta \boldsymbol{z}^\pm = \delta \boldsymbol{v} \pm \frac {\delta \boldsymbol{B}}{\sqrt {4 \mathrm{\pi } m_{\mathrm{p}} n}}. \end{equation}

We take the background magnetic field $\boldsymbol{B}$ to point toward the Sun, and therefore, $\delta \boldsymbol{z}^+$ ( $\delta \boldsymbol{z}^-$ ) represents Alfvénic fluctuations propagating away from (toward) the Sun in the plasma frame. We define the perpendicular and parallel outer scales $L_\perp (r)$ and $L_\parallel ^\pm (r)$ to be the length scales that make the dominant contribution to the fluctuation energy.Footnote ⁵ Although $L_\perp$ can be different (and evolve with $r$ differently) for $\delta \boldsymbol{z}^+$ and $\delta \boldsymbol{z}^-$ (Meyrand et al. Reference Meyrand, Squire, Mallet and Chandran2025), for simplicity, we take it to be the same.

The energy density of outward-propagating AWs is ${\mathcal E}_{\mathrm{w}} = m_{\mathrm{p}}n \left (\delta z^+ _{\mathrm{rms}}\right )^2/4$ , where $\delta z^\pm _{\mathrm{rms}}$ is the root-mean-square (r.m.s.) amplitude of $\delta \boldsymbol{z}^\pm$ . Equivalently,

(3.3) \begin{equation} \delta z^+_{\mathrm{rms}} = 2 \left (\frac {{\mathcal E}_{\mathrm{w}}}{m_{\mathrm{p}} n}\right )^{1/2}. \end{equation}

At $r\lesssim 20 R_{\odot }$ , where most of the heating and acceleration of the solar wind take place, the turbulence is highly imbalanced (e.g. Verdini & Velli Reference Verdini and Velli2007; Perez & Chandran Reference Perez and Chandran2013; Chen et al. Reference Chen2020; McIntyre et al. Reference McIntyre, Chen, Squire, Meyrand and Simon2024), meaning that

(3.4) \begin{equation} \delta z^-_{\mathrm{rms}} \ll \delta z^+_{\mathrm{rms}}. \end{equation}

We derive our turbulent heating model in this highly imbalanced limit. Equations (3.2) and (3.4) imply that the r.m.s. fluctuating velocity is approximately

(3.5) \begin{equation} \delta v_{\mathrm{rms}} = \frac {1}{2} \delta z^+_{\mathrm{rms}}. \end{equation}

We estimate $\delta z^-_{\mathrm{rms}}$ by balancing the rate at which $\delta \boldsymbol{z}^-$ is produced by reflection against the rate at which $\delta \boldsymbol{z}^-$ cascades and dissipates, which yields

(3.6) \begin{equation} \delta z^-_{\mathrm{rms}} = \frac {L_{\perp } (U + v_{\mathrm{A}})}{v_{\mathrm{A}}} \,\left | \frac {\partial v_{\mathrm{A}}}{\partial r}\right |. \end{equation}

Equation (3.6) was originally derived in the strong-turbulence regime by Dmitruk et al. (Reference Dmitruk, Matthaeus, Milano, Oughton, Zank and Mullan2002) for the case without background flow and by Chandran & Hollweg (Reference Chandran and Hollweg2009) for the case with background flow. Subsequently, Chandran & Perez (Reference Chandran and Perez2019) showed that (3.6) is also approximately valid when $\delta \boldsymbol{z}^-$ fluctuations are in the weak-turbulence regime. We discuss the weak- and strong-turbulence regimes further in §§3.3 and 3.4

In the numerical examples that we present in §4, $v_{\mathrm{A}}$ varies so rapidly near $r=r_{\mathrm{b}}$ (the minimum heliocentric distance in the model, approximately 2000 km above the photosphere) that the value of $\delta z^-_{\mathrm{rms}}$ in (3.6) can exceed $\delta z^+_{\mathrm{rms}}$ , violating the assumptions under which (3.6) was derived. To avoid this, we replace $|\partial v_{\mathrm{A}}/\partial r|$ in (3.6) with

(3.7) \begin{equation} \left | \frac {\partial v_{\mathrm{A}}}{\partial r} \right |_{\mathrm{eff}}= \mathrm{min}\left ( \left | \frac {\partial v_{\mathrm{A}}}{\partial r} \right |_{\mathrm{eff}}, \frac { c_{\mathrm{refl}}\delta z^+_{\mathrm{rms}}}{L_\perp }\right ), \end{equation}

where the constant $c_{\mathrm{refl}}$ is an adjustable free parameter. Equation (3.7), in combination with (3.6), prevents $\delta z^-_{\mathrm{rms}}/\delta z^+_{\mathrm{rms}}$ from exceeding $c_{\mathrm{refl}} (U+v_{\mathrm{A}})/v_{\mathrm{A}}$ .

For simplicity, we make the approximation that $L_\perp ^{-1}$ and $\left (L_\parallel ^+\right )^{-1}$ scale with radius in the same way as the perpendicular and parallel wavenumbers of outward-propagating AWs in the WKB (Wentzel–Kramers–Brillouin) approximation in a steady-state solar wind with a radial magnetic field. This yields

(3.8) \begin{equation} L_\perp (r) = L_{\perp \textrm kG} \left [\frac {1 \mbox{ kG}}{B(r)}\right ]^{1/2} \qquad L_\parallel ^+(r,t) = \left [U(r,t)+v_{\mathrm{A}}(r,t)\right ] \tau _{\mathrm{b}}, \end{equation}

where $L_{\perp \textrm kG}$ is a free parameter and

(3.9) \begin{equation} \tau _{\mathrm{b}} = \frac {L_\perp (r_{\mathrm{b}})}{\delta v_{\mathrm{rms}}(r_{\mathrm{b}})} \end{equation}

is the approximate correlation time of the fluctuations at $r=r_{\mathrm{b}}$ .

3.2. The energy cascade in magnetohydrodynamics

We view the fluctuations as a collection of $\delta \boldsymbol{z}^\pm$ wave packets with length scales $\lambda$ and $l^\pm _\lambda$ measured perpendicular and parallel to the magnetic field, respectively, and amplitudes $\delta z^\pm _\lambda$ . We take outer-scale wave packets to satisfy the equations

(3.10) \begin{equation} \delta z^\pm _{L_\perp } = \delta z^\pm _{\mathrm{rms}}, \qquad l^\pm _{L_\perp } = L^\pm _\parallel . \end{equation}

Each $\delta \boldsymbol{z}^\pm$ wave packet propagates at velocity $\mp v_{\mathrm{A}} \boldsymbol{\hat {b}}$ while being distorted by nonlinear interactions with counter-propagating wave packets. As pointed out by Lithwick, Goldreich & Sridhar (Reference Lithwick, Goldreich and Sridhar2007), it is useful to consider the evolution of a ‘slice’ of a $\delta \boldsymbol{z}^\pm$ wave packet at perpendicular scale $\lambda$ , that is, a cross-section of the $\delta \boldsymbol{z}^\pm$ wave packet in the plane perpendicular to $\boldsymbol{B}$ . The time required for this slice to propagate through a counter-propagating $\delta \boldsymbol{z}^\mp$ wave packet of perpendicular scale $\lambda$ is approximately

(3.11) \begin{equation} \tau _{\mathrm{lin}, \lambda }^\pm = \frac {l_\lambda ^\mp }{v_{\textrm {A}}}. \end{equation}

The time scale associated with the shearing of the $\delta \boldsymbol{z}^\pm$ wave packet by a $\delta \boldsymbol{z}^\mp$ wave packet at scale $\lambda$ is

(3.12) \begin{equation} \tau _{\mathrm{nl}, \lambda }^\pm = \frac {\lambda }{ \delta z^\mp _\lambda }. \end{equation}

The critical-balance parameter is (Goldreich & Sridhar Reference Goldreich and Sridhar1995)

(3.13) \begin{equation} \chi ^\pm _{\lambda } = \frac {\tau _{\mathrm{lin,} \lambda }^\mp }{\tau _{\textrm {nl,} \lambda }^\mp } = \frac {\delta z^\pm _\lambda l_\lambda ^\pm }{\lambda v_{\mathrm{A}}}, \end{equation}

where we have adopted the $\pm$ labelling convention of Lithwick et al. (Reference Lithwick, Goldreich and Sridhar2007) (hereafter LGS07) for $\chi ^\pm _\lambda$ . When $\chi ^+_\lambda \gtrsim 1$ , a $\delta \boldsymbol{z}^-$ wave packet at scale $\lambda$ is strongly sheared and distorted before it can propagate through a $\delta \boldsymbol{z}^+$ wave packet at scale $\lambda$ . In this strong-turbulence regime, the time required for a $\delta \boldsymbol{z}^-$ wave packet at perpendicular scale $\lambda$ to pass its energy on to smaller scales is approximately

(3.14) \begin{equation} \tau _{\mathrm{casc,}\lambda }^- = \frac {\lambda }{\delta z^+_\lambda } \qquad (\textrm {if} \,\chi ^+_\lambda \geqslant 1). \end{equation}

On the other hand, when $\chi ^+_\lambda \ll 1$ , each slice of a $\delta z^-$ wave packet at scale $\lambda$ passes through a counter-propagating $\delta z^+$ wave packet at scale $\lambda$ before the slice is strongly distorted. In this weak-turbulence regime, the effects on the $\delta z^-$ wave packet of consecutive ‘collisions’ with $\delta z^+$ wave packets at scale $\lambda$ add incoherently, like the steps in a random walk. A $\delta \boldsymbol{z}^-$ wave packet undergoes a fractional distortion ${\simeq} \chi ^+_\lambda$ during a single collision and is thus strongly distorted after $\sim \left (\chi ^+_\lambda \right )^{-2}$ collisions. As the duration of a single collision is $\sim l_\lambda ^+/v_{\mathrm{A}}$ , the energy of a $\delta \boldsymbol{z}^-$ wave packet at scale $\lambda$ cascades to smaller scales on the time scale (Kraichnan Reference Kraichnan1965; Ng & Bhattacharjee Reference Ng and Bhattacharjee1996, Reference Ng and Bhattacharjee1997; Goldreich & Sridhar Reference Goldreich and Sridhar1997)

(3.15) \begin{equation} \tau _{\textrm {casc,} \lambda }^- = \left (\chi ^+_\lambda \right )^{-2} \frac {l_\lambda ^+}{v_{\mathrm{A}}} = \left (\chi ^+_\lambda \right )^{-1} \frac {\lambda }{\delta z^+_\lambda } \qquad (\textrm {if}\, \chi ^+_\lambda \lt 1). \end{equation}

The rate at which $\delta \boldsymbol{z}^\pm$ energy cascades to smaller scales at scale $\lambda$ is

(3.16) \begin{equation} \epsilon ^\pm _\lambda = \frac {m_{\mathrm{p}} n \left (\delta z^\pm _\lambda \right )^2}{4 \tau ^\pm _{{\mathrm{casc,}} \lambda }}. \end{equation}

We take $\epsilon ^\pm _\lambda$ to be independent of $\lambda$ for $\rho _{\mathrm{p}} \lesssim \lambda \leqslant L_\perp$ , and we use the abbreviated notation $\epsilon ^\pm$ to denote the value of $\epsilon ^\pm _\lambda$ anywhere in this scale range. Applying (3.10) and (3.14) through (3.16) to the outer-scale $\delta \boldsymbol{z}^-$ fluctuations, we find that

(3.17) \begin{equation} \epsilon ^- = \frac {m_{\mathrm{p}}n (\delta z^-_{\mathrm{rms}})^2 \delta z^+_{\mathrm{rms}}}{4L_\perp } \times \mathrm{min}\left (\chi ^+_{L_\perp }, 1 \right ). \end{equation}

Equations (3.6), (3.13) and (3.17) imply that $\epsilon ^- \propto L_\perp$ when $\chi ^+_{L_\perp } \geqslant 1$ and $\epsilon ^- \propto L_\perp ^0$ when $\chi ^+_{L_\perp } \lt 1$ . In both cases, doubling $L_\perp$ doubles $\delta z^-_{\mathrm{rms}}$ in (3.6), and in the strong-turbulence limit this doubles $\epsilon ^-$ in (3.17). However, in the weak-turbulence limit, doubling $L_\perp$ without changing $\delta z^+_{\mathrm{rms}}$ , $v_{\mathrm{A}}$ or $L^+_\parallel$ cuts $\chi ^+_{L_\perp }$ in half and leaves $\epsilon ^-$ in (3.17) unchanged.

Because $\delta \boldsymbol{z}^-$ fluctuations are generated by the reflection of $\delta \boldsymbol{z}^+$ fluctuations and are also sheared by $\delta \boldsymbol{z}^+$ fluctuations, $\delta \boldsymbol{z}^-$ fluctuations are highly coherent in the ` $\delta \boldsymbol{z}^+$ reference frame,’ which propagates away from the Sun at the group velocity of the $\delta \boldsymbol{z}^+$ fluctuations (Velli, Grappin & Mangeney Reference Velli, Grappin and Mangeney1989). In particular, in the $\delta \boldsymbol{z}^+$ reference frame, the $\delta \boldsymbol{z}^-$ fluctuations at scale $\lambda$ remain coherent until the $\delta \boldsymbol{z}^+$ fluctuations at scale $\lambda$ evolve appreciably (LGS07). Lithwick et al. (Reference Lithwick, Goldreich and Sridhar2007) made this last argument for the strong-turbulence regime in which $\chi _\lambda ^+\sim 1$ , but the essence of their argument applies equally well to the weak-turbulence regime of reflection-driven turbulence, in which $\chi ^+_\lambda \lt 1$ . In both regimes, if $\delta \boldsymbol{z}^-$ is infinitesimal so that $\delta \boldsymbol{z}^+$ does not evolve via nonlinear interactions, then (neglecting the radial inhomogeneity of the background) $\delta \boldsymbol{z}^+$ is constant in the $\delta \boldsymbol{z}^+$ frame. As a consequence, $\delta \boldsymbol{z}^-$ is constant in the $\delta \boldsymbol{z}^+$ frame, because $\delta \boldsymbol{z}^-$ is both generated by $\delta \boldsymbol{z}^+$ (via reflection) and cascaded by $\delta \boldsymbol{z}^+$ . If, in a thought experiment, one now increases the amplitude of $\delta \boldsymbol{z}^-$ so that $\delta \boldsymbol{z}^+$ evolves via nonlinear interactions, the question arises: How long does one have to wait before the changes in $\delta \boldsymbol{z}^+$ give rise to changes in $\delta \boldsymbol{z}^-$ in the $\delta \boldsymbol{z}^+$ frame? Noting that the nonlinear evolution time of $\delta \boldsymbol{z}^+$ fluctuations is an increasing function of $\lambda$ , LGS07 argued that the $\delta \boldsymbol{z}^-$ fluctuations at scale $\lambda$ do not evolve in the $\delta \boldsymbol{z}^+$ frame until the $\delta \boldsymbol{z}^+$ fluctuations at comparable scales have evolved. Prior to such a time, $\delta \boldsymbol{z}^+$ fluctuations at scales ${\ll} \lambda$ will have changed appreciably, but such changes will have little effect upon the $\delta \boldsymbol{z}^-$ fluctuations at scale $\lambda$ , which are sensitive primarily to the reflection of $\delta \boldsymbol{z}^+$ fluctuations at scales $\gtrsim \lambda$ and the way that such $\delta \boldsymbol{z}^+$ fluctuations cause $\delta \boldsymbol{z}^-$ fluctuations at scales $\gtrsim \lambda$ to cascade. Therefore, a $\delta \boldsymbol{z}^+$ fluctuation at scale $\lambda$ is sheared coherently by $\delta \boldsymbol{z}^-$ fluctuations at scale $\lambda$ throughout the lifetime of the $\delta \boldsymbol{z}^+$ fluctuation, and the energy-cascade time scale of $\delta z_\lambda ^+$ fluctuations is approximately

(3.18) \begin{equation} \tau _{\mathrm{casc,} \lambda }^+ = \frac {\lambda }{\delta z^-_{\lambda }}, \end{equation}

regardless of the value of $\chi ^-_\lambda$ . We refer the reader to LGS07 for a more detailed discussion of this turbulence phenomenology as it plays out in the strong-turbulence regime.

Applying (3.16) and (3.18) to fluctuations at the outer scale, and making use of (3.6) and (3.7), we find the magnetohydrodynamic (MHD) prediction for the rate at which $\delta \boldsymbol{z}^+$ energy cascades to smaller scales (cf. Dmitruk et al. Reference Dmitruk, Matthaeus, Milano, Oughton, Zank and Mullan2002; Chandran & Hollweg Reference Chandran and Hollweg2009):

(3.19) \begin{equation} \epsilon ^+ = \frac {m_{\mathrm{p}} n}{4}\left (\frac {U + v_{\mathrm{A}}}{v_{\mathrm{A}}}\right ) \left | \frac {\partial v_{\mathrm{A}}}{\partial r}\right |_{\mathrm{eff}} \left (\delta z^+_{\mathrm{rms}}\right )^2. \end{equation}

Just as $\epsilon ^-$ is independent of $L_\perp$ when $\chi ^+_{L_\perp } \lt 1$ in the sense described following (3.17), $\epsilon ^+$ is independent of $L_\perp$ for fixed $U$ , $v_{\mathrm{A}}$ , $\partial v_{\mathrm{A}}/\partial r$ and $\delta z^+_{\mathrm{rms}}$ .

3.3. The possibility of weak turbulence at the large-scale end of the inertial range

When $\chi _{L_\perp }^+ \lt 1$ , $\delta \boldsymbol{z}^-$ fluctuations with $\lambda \sim L_\perp$ are weakly turbulent. In this subsection we describe how we model weak turbulence at the large-scale end of the inertial range when it arises. In the numerical examples in §4, weak turbulence is limited to regions outside the low corona and to a comparatively narrow range of scales near $L_\perp$ for reasons that are explained towards the end of this subsection. As a consequence, the weak-turbulence scalings that we derive here do not have a large impact on our numerical solutions in §4. Nevertheless, we include this discussion for completeness and because, to our knowledge, there are no existing models for the scale dependence of $\delta z^\pm _\lambda$ and $l^\pm _\lambda$ in weak anisotropic reflection-driven turbulence.

To obtain weak-turbulence scalings, we consider the asymptotic weak-turbulence limit, in which

(3.20) \begin{equation} \chi ^+_\lambda \ll 1, \end{equation}

and we take $l^-_\lambda$ to be the distance a $\delta \boldsymbol{z}^-$ wave packet at perpendicular scale $\lambda$ can propagate along the magnetic field during its cascade time scale. It then follows from (3.15) that

(3.21) \begin{equation} l_\lambda ^- = l_\lambda ^+ \left ( \chi ^+_\lambda \right )^{-2}. \end{equation}

(Setting $\lambda = L_\perp$ in (3.21) implies, via (3.10), that $L_\parallel ^- = L_\parallel ^+ \left (\chi ^+_{L_\perp }\right )^{-2}$ when $\chi ^+_{L_\perp } \ll 1$ .) Because $l_\lambda ^- \gg l_\lambda ^+$ , the shearing of $\delta \boldsymbol{z}^+$ wave packets at perpendicular scale $\lambda$ by $\delta \boldsymbol{z}^-_\lambda$ wave packets at perpendicular scale $\lambda$ does not reduce the parallel correlation length of the $\delta \boldsymbol{z}^+$ wave packets at scale $\lambda$ , and thus, we set

(3.22) \begin{equation} l^+_\lambda \propto \lambda ^0. \end{equation}

As mentioned following (3.16), in the inertial range, the energy-cascade rate is independent of $\lambda$ . Upon substituting (3.15) into (3.16) and taking $\epsilon ^-_\lambda$ to be independent of $\lambda$ , we find that

(3.23) \begin{equation} \epsilon ^-_\lambda = \left (\frac {m_{\mathrm{p}} n}{4}\right ) \frac { \left (\delta z^-_\lambda \right )^2 \left (\delta z^+_\lambda \right )^2 l_\lambda ^+}{\lambda ^2 v_{\mathrm{A}} } \propto \lambda ^0. \end{equation}

Likewise, after substituting (3.18) into (3.16)Footnote ⁶ and taking $\epsilon ^+$ to be independent of $\lambda$ , we obtain

(3.24) \begin{equation} \epsilon ^+_\lambda = \left (\frac {m_{\mathrm{p}} n}{4}\right ) \frac { \left (\delta z^+_\lambda \right )^2 \delta z^-_\lambda }{\lambda } \propto \lambda ^0. \end{equation}

Dividing (3.23) by (3.24) and making use of (3.22) yields

(3.25) \begin{equation} \delta z^-_\lambda \propto \lambda , \end{equation}

which corresponds to a perpendicular $\boldsymbol{z}^-$ power spectrum $E^-(k_\perp ) \propto k_\perp ^{-3}$ . Upon substituting (3.25) into either (3.23) or (3.24), we find that

(3.26) \begin{equation} \delta z^+_\lambda \propto \lambda ^0, \end{equation}

which corresponds to a perpendicular $\boldsymbol{z}^+$ power spectrum $E^+(k_\perp ) \propto k_\perp ^{-1}$ . It then follows from (3.13), (3.22) and (3.26) that

(3.27) \begin{equation} \chi ^+_\lambda \propto \lambda ^{-1}. \end{equation}

Equation (3.27) implies that if the turbulence starts out in the weak-turbulence regime at the outer scale with $\chi ^+_{L_\perp } \lt 1$ , then the turbulence transitions to the strong-turbulence regime at $\lambda = \chi ^+_{L_\perp } L_\perp$ . For arbitrary values of $\chi ^+_{L_\perp }$ , the scale

(3.28) \begin{equation} \lambda _{\mathrm{str}} = L_\perp \times \mathrm{min}\left (\chi ^+_{L_\perp }, 1\right ) \end{equation}

then marks the large-scale end of the strongly turbulent part of the inertial range.

3.4. Intermittent, strong MHD turbulence at $\rho _{\mathrm{p}} \lt \lambda \lt \lambda _{\mathrm{str}}$

Parker Solar Probe observations reveal that strong turbulence in the near-Sun solar wind is intermittent (Sioulas et al. Reference Sioulas2024). In intermittent turbulence, the majority of the fluctuation energy at scale $\lambda$ is concentrated into a fraction of the volume that decreases as $\lambda$ decreases. Equivalently, as $\lambda$ decreases, the probability distribution function (PDF) of the fluctuation amplitudes at scale $\lambda$ broadens, and the tail of this distribution accounts for an increasing fraction of the fluctuation energy at scale $\lambda$ . Accounting for intermittency is important in solar-wind models because it affects how the turbulent heating power is apportioned between protons and electrons, and between parallel and perpendicular proton heating (Mallet et al. Reference Mallet, Klein, Chand, Hoppock, Bowen, Salem and Bale2019), as we will describe in detail in §3.5.

To model intermittent turbulence, we need a mathematical model for the PDF and how it varies with $\lambda$ . For this, we adopt the model of Chandran et al. (Reference Chandran, Sioulas, Bale, Bowen, David, Meyrand and Yerger2025), which agrees with a number of PSP observations. In the remainder of this subsection, we summarise the key ideas of this model and the equations from the model that we use in §3.5 to evaluate the turbulent heating rate.

Chandran et al. (Reference Chandran, Sioulas, Bale, Bowen, David, Meyrand and Yerger2025) considered the Elsasser increments of scale $\lambda$ at position $\boldsymbol{x}$ and time $t$ ,

(3.29) \begin{equation} \Delta \boldsymbol{z}^\pm _\lambda \left (\boldsymbol{x}, \boldsymbol{\hat {s}}, t\right )= \boldsymbol{z}^\pm \left (\boldsymbol{x} + 0.5 \lambda \boldsymbol{\hat {s}}, t \right ) - \boldsymbol{z}^\pm \left (\boldsymbol{x} - 0.5 \lambda \boldsymbol{\hat {s}}, t \right ), \end{equation}

where $\boldsymbol{\hat {s}}$ is a unit vector perpendicular to $\boldsymbol{B}$ . They then defined the characteristic amplitude of the $\delta \boldsymbol{z}^\pm$ structure of perpendicular scale $\lambda$ at position $\boldsymbol{x}$ and time $t$ to be

(3.30) \begin{equation} \delta z^\pm _\lambda (\boldsymbol{x},t) = \frac {1}{2\mathrm{\pi }} \int _0^{2\mathrm{\pi }} \mathrm{d}\theta \, \left |\Delta \boldsymbol{z}^\pm _\lambda \left (\boldsymbol{x}, \boldsymbol{\hat {s}}, t\right )\right |, \end{equation}

where the angle $\theta$ specifies the direction of $\boldsymbol{\hat {s}}$ within the plane perpendicular to $\boldsymbol{B}$ . Rather than considering the detailed $\boldsymbol{x}$ and $t$ dependence of $\delta z^\pm _\lambda (\boldsymbol{x},t)$ , Chandran et al. (Reference Chandran, Sioulas, Bale, Bowen, David, Meyrand and Yerger2025) viewed the fluctuations at each $\lambda$ as a statistical ensemble and treated $\delta z^\pm _\lambda (\boldsymbol{x},t)$ (henceforth simply $\delta z^\pm _\lambda$ ) as a random variable given by

(3.31) \begin{equation} \delta z_\lambda ^\pm = \overline { z}^\pm \beta ^q, \end{equation}

where $\beta$ is a constant $\in (0,1)$ , $q$ is a random integer with a Poisson distribution

(3.32) \begin{equation} P(q) = \frac {{\rm e}^{-\mu } \mu ^q}{q!} \end{equation}

and $\mu$ is the scale-dependent mean of $q$ . Chandran et al. (Reference Chandran, Sioulas, Bale, Bowen, David, Meyrand and Yerger2025) treated $\overline { z}^\pm$ in (3.31) as a scale-independent random number with an arbitrary distribution, and we simply set

(3.33) \begin{equation} \overline { z}^\pm = \delta z^\pm _{\mathrm{rms}}. \end{equation}

As in previous studies (Grauer, Krug & Marliani Reference Grauer, Krug and Marliani1994; Politano & Pouquet Reference Politano and Pouquet1995; Mallet & Schekochihin Reference Mallet and Schekochihin2017), Chandran et al. (Reference Chandran, Sioulas, Bale, Bowen, David, Meyrand and Yerger2025) determined $\mu$ by assuming that the most intense fluctuations (those with $q=0$ ) are sheet-like with a volume-filling factor $P(0) \propto \lambda$ , obtaining

(3.34) \begin{equation} \mu = A + \ln \left (\frac {\lambda _{\mathrm{str}}}{\lambda }\right ), \end{equation}

where $A$ is a constant that determines the breadth of the amplitude distribution at the largest strongly turbulent scale. We set $A=0$ , which, in conjunction with (3.33), amounts to taking the turbulence to be characterised by a single amplitude (rather than a broad distribution) at scale $\lambda _{\mathrm{str}}$ . Qualitatively, as $\lambda$ decreases, $\mu$ increases, causing the distribution of $q$ to broaden, which in turn broadens the PDF of the fluctuation amplitudes.

Chandran et al. (Reference Chandran, Sioulas, Bale, Bowen, David, Meyrand and Yerger2025) solved for $\beta$ by taking the average $\delta \boldsymbol{z}^+$ cascade rate to be independent of $\lambda$ , obtaining

(3.35) \begin{equation} \beta = 2 W_0(1/2) = 0.7035, \end{equation}

where $W_0$ is the Lambert $W$ function. To obtain (3.35), Chandran et al. (Reference Chandran, Sioulas, Bale, Bowen, David, Meyrand and Yerger2025) assumed that $\epsilon ^+$ is dominated by the $\delta z^+_\lambda$ fluctuations in the tail of the distribution (small $q$ values), for which $\delta z^+_\lambda \gt w^+_\lambda$ , where

(3.36) \begin{equation} w^\pm _{\lambda } = \overline { z}^\pm \beta ^\mu = \overline { z}^\pm \left (\frac {\lambda }{\lambda _{\mathrm{str}}}\right )^{-\ln \beta } \end{equation}

is the approximate median value of $\delta z^\pm _\lambda$ . They further argued that, where $\delta z^+_\lambda \gt w^+_\lambda$ , $\delta z^-_\lambda$ can be reasonably approximated as

(3.37) \begin{equation} \delta z^-_\lambda = \frac {w^-_\lambda w^+_\lambda }{\delta z^+_\lambda }. \end{equation}

In (3.37), $\delta z^-_\lambda$ is inversely proportional to $\delta z^+_\lambda$ in the tail of the $\delta z^+_\lambda$ distribution because of the strong shearing (and rapid cascading) experienced by $\delta z^-_\lambda$ where $\delta z^+_\lambda$ is large. The average $\delta \boldsymbol{z}^+$ cascade power at scale $\lambda$ is then approximately

(3.38) \begin{equation} \left \langle \epsilon ^+_\lambda \right \rangle = \sum _{q=0}^\infty P(q) \frac {m_{\mathrm{p}} n \left (\delta z^+_\lambda \right )^2 \delta z^-_\lambda }{4 \lambda } = \sum _{q=0}^\infty P(q) \epsilon ^+_{\lambda ,q}, \end{equation}

where

(3.39) \begin{equation} \epsilon ^+_{\lambda , q} = \frac { m_{\mathrm{p}} n \left (\overline { z}^+\right )^2 \overline { z}^-}{4 \lambda } \left ( \frac {\lambda }{\lambda _{\mathrm{str}}}\right )^{-2\ln \beta } \beta ^q . \end{equation}

Equations (3.38) and (3.39) quantify how much each part of the PDF of $\delta z^+_\lambda$ contributes to $\left \langle \epsilon ^+_\lambda \right \rangle$ . To determine the parallel length scales of the AW packets, Chandran et al. (Reference Chandran, Sioulas, Bale, Bowen, David, Meyrand and Yerger2025) applied the arguments of LGS07 to intermittent turbulence, setting

(3.40) \begin{equation} l^+_\lambda = l^-_\lambda = \frac {v_{\mathrm{A}} \lambda }{\delta z^+_\lambda } \end{equation}

throughout the $\delta z^+_\lambda$ distribution.

3.5. The turbulent heating rate

We make the simplifying assumption that all dissipation takes place either at $\lambda \sim \rho _{\mathrm{p}}$ (proton-gyroscale dissipation) or at $\lambda \ll \rho _{\mathrm{p}}$ (sub-proton-scale dissipation). We note that $\tau _{\mathrm{casc,}\lambda }^- \ll \tau _{\mathrm{casc,} \lambda }^+$ and assume that $\tau _{\mathrm{casc,}\lambda }^-$ is much smaller than the dissipation time scale of $\delta \boldsymbol{z}^-$ fluctuations at $\lambda \sim \rho _{\mathrm{p}}$ . Given this assumption, $\epsilon ^-$ cascades past $\lambda \sim \rho _{\mathrm{p}}$ and dissipates almost entirely at $\lambda \ll \rho _{\mathrm{p}}$ , heating only the electrons. We take the perpendicular proton heating rate $Q_{\perp \textrm {p}}$ to be the rate at which stochastic proton heating (e.g. McChesney, Stern & Bellan Reference McChesney, Stern and Bellan1987) removes energy from $\delta \boldsymbol{z}^+$ fluctuations at $\lambda \sim \rho _{\mathrm{p}}$ , and we equate $Q_{\parallel \textrm {p}}$ with the rate at which proton Landau damping (LD) and transit-time damping (TTD) remove energy from $\delta \boldsymbol{z}^+$ fluctuations at $\lambda \sim \rho _{\mathrm{p}}$ . We define $Q_{\mathrm{e1}}$ to be the rate at which electron LD and TTD remove energy from $\delta \boldsymbol{z}^+$ fluctuations at $\lambda \sim \rho _{\mathrm{p}}$ and take the total electron heating rate to be

(3.41) \begin{equation} Q_{\mathrm{e}} = Q_{\mathrm{e1}} + Q_{\mathrm{e, sub-proton}}, \end{equation}

where

(3.42) \begin{equation} Q_{\mathrm{e, sub-proton}} = \epsilon ^- + \epsilon ^+ - Q_{\perp \textrm {p}} - Q_{\parallel \textrm {p}} - Q_{\mathrm{e1}} \end{equation}

is the rate at which energy cascades past the proton gyroradius scale to $\lambda \ll \rho _{\mathrm{p}}$ .

As a preliminary step towards evaluating $Q_{\perp \textrm {p}}$ , $Q_{\parallel \textrm {p}}$ and $Q_{\mathrm{e1}}$ , we relate $\lambda$ and $l^+_\lambda$ to characteristic perpendicular and parallel wavenumbers, $k_\perp$ and $k_\parallel (k_\perp )$ , using the relations suggested by figure 1 of Huang et al. (Reference Huang2023):

(3.43) \begin{equation} k_\perp =\frac {\mathrm{\pi }}{\lambda }, \qquad k_\parallel (k_\perp ) = \frac {\mathrm{\pi }}{l_\lambda ^\pm }. \end{equation}

As described further in the following, we evaluate the proton-gyroscale heating rates using the properties of the turbulence at

(3.44) \begin{equation} k_\perp \rho _{\mathrm{p}} = 1 \;\; \;\;\longleftrightarrow \;\;\;\; \lambda = \mathrm{\pi } \rho _{\mathrm{p}}. \end{equation}

Equation (3.43) and the critical-balance relation (3.40) imply that

(3.45) \begin{equation} k_\parallel (\rho _{\mathrm{p}}^{-1}) v_{\mathrm{A}} = \frac {\mathrm{\pi } v_{\mathrm{A}}}{l_{\mathrm{\pi } \rho _{\mathrm{p}}}} = \frac {\mathrm{\pi } \delta z^+_{\mathrm{p}}}{\mathrm{\pi } \rho _{\mathrm{p}}} = \frac {\delta z^+_{\mathrm{p}}}{\rho _{\mathrm{p}}}, \end{equation}

where $\delta z^+_{\mathrm{p}}$ is the value of $\delta z^+_\lambda$ at $\lambda = \mathrm{\pi } \rho _{\mathrm{p}}$ . Although fluctuations at $k_\perp \rho _{\mathrm{p}} = 1$ are at the transition between MHD scales and kinetic scales, we approximate these fluctuations using results from MHD. In particular, we assume that the PDF of $\delta z^+_{\mathrm{p}}$ is given by (3.31), (3.32) and (3.34) with $\lambda = \mathrm{\pi } \rho _{\mathrm{p}}$ , that the energy density of these fluctuations is $m_{\mathrm{p}} n \left (\delta z^+_{\mathrm{p}}\right )^2/4$ , and that the amplitude of the fluctuating electron fluid velocity (or $\boldsymbol{E}\boldsymbol{\times } \boldsymbol{B}$ velocity) associated with the $\delta z^+_{\mathrm{p}}$ fluctuations at $k_\perp \rho _{\mathrm{p}} =1$ is given by

(3.46) \begin{equation} \delta v_{\mathrm{p}} = \frac {1}{2} \delta z^+_{\mathrm{p}}. \end{equation}

We define $\gamma _{\mathrm{p}}(k_\perp )$ and $\gamma _{\mathrm{e}}(k_\perp )$ to be the rates at which (kinetic) AWs at perpendicular wavenumber $k_\perp$ and parallel wavenumber $k_\parallel (k_\perp )$ lose energy to linear LD/TTD via interactions with protons and electrons, respectively. We adopt the analytic expressions for $\gamma _{\mathrm{p}}(k_\perp )$ and $\gamma _{\mathrm{e}}(k_\perp )$ derived by Howes et al. (Reference Howes, Cowley, Dorland, Hammett, Quataert and Schekochihin2006) within the gyrokinetics approximation for the case in which $m_{\mathrm{e}} T_{\mathrm{p}}/(m_{\mathrm{p}} T_{\mathrm{e}}) \ll \beta _{\mathrm{p}} \ll 1$ , simplified to the limit in which $k_\perp \rho _{\mathrm{e}} \ll 1$ , where $\rho _{\mathrm{e}}$ is the electron gyroradius:

(3.47) \begin{equation} \gamma _{\mathrm{e}}(k_\perp ) = \frac {\mathrm{\pi }^{1/2}}{2} \left (\frac { \alpha _{\mathrm{p}}^2 m_{\mathrm{e}} T_{\mathrm{e}}}{m_{\mathrm{p}} T_{\mathrm{p}} \beta _{\mathrm{p}}}\right )^{1/2} \left | k_\parallel (k_\perp ) \right | v_{\mathrm{A}} \end{equation}

and

(3.48) \begin{equation} \gamma _{\mathrm{p}}(k_\perp ) = \frac {\mathrm{\pi }^{1/2}}{2} \frac {\alpha _{\mathrm{p}} T_{\mathrm{e}}^2}{ \beta _{\mathrm{p}}^{1/2} T_{\mathrm{p}}^2} e^{- \alpha _{\mathrm{p}} - \overline {\omega }^2 \beta _{\mathrm{p}}^{-1}} I_0(\alpha _{\mathrm{p}}) \left |k_\parallel (k_\perp )\right | v_{\mathrm{A}}, \end{equation}

where

(3.49) \begin{equation} \alpha _{\mathrm{p}} = \frac {1}{2} \left (k_\perp \rho _{\mathrm{p}}\right )^2, \qquad \beta _{\mathrm{p}} = \frac {8 \mathrm{\pi } n k_{\mathrm{B}}T_{\mathrm{p}}}{B^2}, \end{equation}

$T_{\mathrm{p}}$ is defined in (2.18), $I_0(x)$ is the modified Bessel function of the first kind of order 0 and $\overline {\omega } = \omega /(k_\parallel v_{\mathrm{A}})$ . Like $\delta z^+_{\mathrm{p}}$ , $\gamma _{\mathrm{e}}(\rho_{\rm p}^{-1})$ and $\gamma _{\parallel \textrm {p}}(\rho_{\rm p}^{-1})$ are functions of $q$ via (3.45).

Using these linear damping rates, we define hypothetical rates $\Gamma _{\mathrm{e},q}$ and $\Gamma _{\parallel \mathrm{p},q}$ at which electrons and protons, respectively, would drain energy from $\delta z^+_{ \lambda }$ fluctuations at $\lambda = \mathrm{\pi } \rho _{\mathrm{p}}$ (i.e. $\delta z^+_{\mathrm{p}}$ fluctuations) with a given value of $q$ in (3.31) via LD and TTD if the fluctuation amplitudes were unaffected by this energy transfer:

(3.50) \begin{equation} \Gamma _{\mathrm{e}, q} = \frac {c_{\mathrm{LD}}}{2}\gamma _{\mathrm{e}}\left (\rho _{\mathrm{p}}^{-1}\right ) m_{\mathrm{p}} n \left (\delta z^+_{\mathrm{p}}\right )^2 \end{equation}

and

(3.51) \begin{equation} \Gamma _{\parallel \mathrm{p}, q} = \frac {1}{2} \gamma _{\parallel \textrm p}\left (\rho _{\mathrm{p}}^{-1}\right ) m_{\mathrm{p}} n \left (\delta z^+_{\mathrm{p}}\right )^2,\end{equation}

where the constant $c_{\mathrm{LD}}$ is a free parameter that enables us to modify the effective electron LD rate.Footnote ⁷ In writing (3.50) and (3.51), we have accounted for the fact that the energy damping rate is twice the amplitude damping rate. For perpendicular proton heating, we define $\Gamma _{\perp \textrm {p},q}$ to be the hypothetical rate at which $\delta z^+_{\mathrm{p}}$ fluctuations would cause stochastic perpendicular proton heating if there were no reduction of $\delta z^+_{\mathrm{p}}$ by dissipation at $k_\perp \rho _{\mathrm{p}} \sim 1$ . Following Chandran et al. (Reference Chandran, Li, Rogers, Quataert and Germaschewski2010), we set

(3.52) \begin{equation} \Gamma _{\perp \mathrm{p},q} = \frac {m_{\mathrm{p}} n \,\delta v_{\mathrm{p}}^3}{\rho _{\mathrm{p}}} \exp \left ( - \frac {c_2 v_{\perp \textrm th, p}}{\delta v_{\mathrm{p}}}\right ), \end{equation}

where $c_2$ is an adjustable parameter, and we have set $c_1=1$ in equation (25) of Chandran et al. (Reference Chandran, Li, Rogers, Quataert and Germaschewski2010) to reduce the number of free parameters in our model.Footnote ⁸

To go beyond the hypothetical heating rates in (3.50) through (3.52), we need to account for the feedback of the damping/dissipation process on the proton-gyroscale fluctuations, which prevents the proton-gyroscale turbulent heating rate from exceeding the rate at which energy cascades through the inertial range. We do so by taking the contribution of the $\delta z^+_\lambda$ fluctuations at $\lambda = \pi \rho _{\mathrm{p}}$ with a given value of $q$ in (3.31) to $Q_{\perp \textrm {p}}$ , $Q_{\parallel \textrm p}$ and $Q_{\mathrm{e1}}$ to be, respectively,

(3.53) \begin{equation} Q_{\perp \mathrm{p}, q} = \frac {\Gamma _{\perp \mathrm{p}, q} \epsilon ^+_{\mathrm{p},q}}{\epsilon ^+_{\mathrm{p},q} + \Gamma _{\perp \mathrm{p}, q} + \Gamma _{\parallel \mathrm{p},q} + \Gamma _{\mathrm{e}, q}}, \end{equation}

(3.54) \begin{equation} Q_{\parallel \mathrm{p}, q} = \frac {\Gamma _{\parallel \mathrm{p}, q} \epsilon ^+_{\mathrm{p},q}}{\epsilon ^+_{\mathrm{p},q} + \Gamma _{\perp \mathrm{p}, q} + \Gamma _{\parallel \mathrm{p},q} + \Gamma _{\mathrm{e},q}} \end{equation}

and

(3.55) \begin{equation} Q_{\mathrm{e}1, q} = \frac {\Gamma _{\mathrm{e},q} \epsilon ^+_{\mathrm{p},q}}{\epsilon ^+_{\mathrm{p},q} + \Gamma _{\perp \mathrm{p}, q} + \Gamma _{\parallel \mathrm{p},q} + \Gamma _{\mathrm{e},q}}, \end{equation}

where $\epsilon ^+_{\mathrm{p},q}$ is the value of $\epsilon ^+_{\lambda ,q}$ in (3.39) with $\lambda = \mathrm{\pi } \rho _{\mathrm{p}}$ , which is the rate at which the energy cascade replenishes the energy of the fluctuations with a given value of $q$ at $\lambda = \mathrm{\pi } \rho _{\mathrm{p}}$ . If $\Gamma _{\perp \mathrm{p},q} + \Gamma _{\parallel \mathrm{p}, q} + \Gamma _{\mathrm{e},q}$ is $\ll \epsilon ^+_{\mathrm{p},q}$ , then the energy drained from fluctuations at $k_\perp \rho _{\mathrm{p}} \sim 1$ has a negligible effect on their amplitudes, and (3.53) through (3.55) imply that $Q_{\perp \mathrm{p},q} \simeq \Gamma _{\perp \mathrm{p},q}$ , $Q_{\parallel \mathrm{p},q} \simeq \Gamma _{\parallel \mathrm{p},q}$ and $Q_{\mathrm{e}1,q} \simeq \Gamma _{\mathrm{e},q}$ . Conversely, if $\Gamma _{\perp \mathrm{p},q} + \Gamma _{\parallel \mathrm{p}, q} + \Gamma _{\mathrm{e},q} \gg \epsilon ^+_{\mathrm{p}, q}$ then the energy drained from the proton-gyroscale fluctuations strongly modifies their amplitudes and $Q_{\perp \mathrm{p},q} \ll \Gamma _{\perp \mathrm{p},q}$ , $Q_{\parallel \mathrm{p},q} \ll \Gamma _{\parallel \mathrm{p},q}$ and $Q_{\mathrm{e}1,q} \ll \Gamma _{\mathrm{e},q}$ . In this regime, (3.53) through (3.55) imply that $Q_{\perp \mathrm{p},q} + Q_{\parallel \mathrm{p},q} + Q_{\mathrm{e1}} \simeq \epsilon ^+_q$ , $Q_{\perp \mathrm{p},q}/Q_{\parallel \mathrm{p},q} = \Gamma _{\perp \mathrm{p},q} / \Gamma _{\parallel \mathrm{p},q}$ and $Q_{\perp \mathrm{p},q}/Q_{\mathrm{e},q} =\Gamma _{\perp \mathrm{p},q} / \Gamma _{\mathrm{e},q}$ . The total heating rates from the full distribution of proton-gyroscale fluctuations are

(3.56) \begin{equation} Q_{\perp \textrm {p}} = \sum _{q=0}^\infty P(q) Q_{\perp \mathrm{p},q}, \qquad Q_{\parallel \textrm {p}} = \sum _{q=0}^\infty P(q) Q_{\parallel \mathrm{p},q}, \qquad Q_{\mathrm{e1 }} = \sum _{q=0}^\infty P(q) Q_{\mathrm{e}1,q}. \end{equation}

We emphasise that the turbulent heating described by (3.50) through (3.56) results entirely from fluctuations at $\lambda = \pi \rho _{\mathrm{p}}$ . We neglect turbulent heating from larger-scale fluctuations and model the turbulent heating from smaller-scale fluctuations via (3.42).

4.1. Magnetic-field model and radial grid

All the examples presented in this section are based on the same axisymmetric magnetic-field model, which is designed to emulate the solar and heliospheric magnetic field near the minimum of the solar cycle. Following Cranmer et al. (Reference Cranmer, van Ballegooijen and Edgar2007), we start with the dipole-plus-quadrupole-plus-current-sheet model of Banaszkiewicz et al. (Reference Banaszkiewicz, Axford and McKenzie1998), adopting the latter authors’ favoured parameters of $K=1$ , $M=1.789$ , $a_1 = 1.538$ and $Q=1.5$ . We then add to this model the low-solar-atmosphere magnetic-field model of Hackenberg et al. (Reference Hackenberg, Marsch and Mann2000) using the parameters from their figure 1: $L=3\times 10^9 \mbox{ cm}$ , $d= 3.4\times 10^7 \mbox{ cm}$ and $B_{\mathrm{max}} = 1.5 \mbox{ kG}$ . We illustrate the magnetic-field lines in this Banaszkiewicz--Hackenberg (B-H) model in figure 2(a).

Our solar-wind model applies to narrow magnetic flux tubes directed radially outwards from the Sun and does not take into account the curvature of magnetic-field lines. In order to approximate the B–H magnetic-field structure with our radial magnetic flux tubes, we compute $B$ as a function of arc length $s$ in the B-H model along a number of magnetic-field lines that intersect the Sun at spherical polar angles $\theta _{\odot }$ and that reach heliolatitudes $\Theta _{60 R_{\odot }}$ at a heliocentric distance $r=60 R_{\odot }$ , where the heliolatitude is $90^\circ$ minus the spherical polar angle. For example, $\theta _{\odot }=0$ and $\Theta _{60 R_{\odot }} = 90^\circ$ both correspond to the radial magnetic-field line connected to the Sun’s north pole. We also compute the function $r(s)$ , the heliocentric distance as a function of arc length, along these field lines as well as the inverse of this function, $s(r)$ . We then compose the function $B(r) = B(s(r))$ for each of these field lines and solve (2.5) through (2.12) within radial magnetic flux tubes that have the same $B(r)$ profiles. Figure 2(b) shows the $B(r)$ profiles for three different values of $\theta _{\odot }$ , which correspond to the three values of $\Theta _{60 R_{\odot }}$ listed in the figure. In the figures to follow, we label our model solutions by the $\Theta _{60 R_{\odot }}$ value of the B–H magnetic-field line that our model attempts to emulate, but we omit the $60 R_{\odot }$ subscript to simplify the figure legends.

To construct our numerical grid of $r$ values, we first take $N=1002$ points that are evenly spaced in the variable $\ln ((s-R_{\odot })/R_{\odot })$ between a minimum value of $\ln (2000 \mbox{ km}/ R_{\odot })$ and a maximum value of $\ln (71)$ . These $N$ points correspond to $N$ arclengths $s_j$ , with $j$ ranging from 1 to $N$ , which we convert into $N$ heliocentric distances $r_j$ by setting $r_j = r(s_j)$ . The resulting $r_j$ values range from approximately $R_{\odot } + 2000 \mbox{ km}$ to $72 R_{\odot }$ . The spacing between neighbouring grid points increases from ${\simeq} 20 \mbox{ km}$ near the inner edge of the grid to ${\simeq} 5\times 10^5 \mbox{ km}$ near the outer edge of the grid.

Figure 2. (a) Magnetic-field lines in the axisymmetric, solar-minimum, B-H magnetic-field model (Banaszkiewicz et al. Reference Banaszkiewicz, Axford and McKenzie1998; Hackenberg et al. Reference Hackenberg, Marsch and Mann2000). (b) The magnetic-field strength $B$ as a function of heliocentric distance $r$ along three different B–H magnetic-field lines that intersect the Sun at spherical polar angle $\theta _{\odot }$ and that connect to heliolatitudes $\Theta _{60 R_{\odot }}$ at $r=60 R_{\odot }$ . We use these $B(r)$ profiles for the radial magnetic flux tubes in our 1-D solar-wind model to emulate magnetic-field lines in the 2-D B–H model. The vertical line shows the minimum radius $r_{\mathrm{b}}$ in our solar-wind model.

4.2. Parameter values

We focus on the set of parameters listed in table 1. We follow van Ballegooijen & Asgari-Targhi (Reference van Ballegooijen and Asgari-Targhi2016) in viewing the open magnetic flux in the heliosphere as originating primarily in magnetic flux tubes that, at the photosphere, have diameters ${\simeq} 100 \mbox{ km}$ and magnetic-field strengths ${\simeq} 1 \mbox{ kG}$ . These flux tubes cover approximately $ 1\,\%$ of the photosphere but expand and merge in the low corona, where $B\simeq 10 \mbox{ G}$ and the flux tubes fill most of the volume. By flux conservation, the diameters of the flux tubes are ${\simeq} 10^3 \mbox{ km}$ in the low corona. Photospheric motions at the base of these flux tubes launch waves that propagate outwards through the chromosphere and into the corona. Because the motions of different flux tubes are uncorrelated at the photosphere, the motions of different flux tubes are largely uncorrelated in the corona. We therefore set $L_{\perp \textrm kG} = 100 \mbox{ km}$ . It follows from (3.8) that this choice of $L_{\perp \textrm kG}$ corresponds to a perpendicular correlation length $L_\perp (r)$ that is $10^3 \mbox{ km}$ in the low corona where $B\simeq 10 \mbox{ G}$ . We note that a very similar value of $L_{\perp \textrm kG}$ ( $ 90.9 \mbox{ km}$ ) in the AW-driven solar-wind model of Cranmer et al. (Reference Cranmer, van Ballegooijen and Edgar2007) led to good agreement with measurements of the Faraday rotation of linearly polarised radio transmissions from the Helios spacecraft (Hollweg, Cranmer & Chandran Reference Hollweg, Cranmer and Chandran2010).

Table 1. Parameter Values.

We determine the remaining parameters approximately as follows. We choose $c_{\mathrm{refl}}$ to match observations of $T_{\mathrm{e}}$ in the low corona, and we then adjust $c_{\mathrm{LD}}$ and $c_2$ together to match observational constraints on $U(r)$ at large $r$ and large $\Theta$ from the Ulysses spacecraft (see §4.3), as well as observational constraints on $T_{\perp \textrm {p}}$ at large $r$ from PSP’s tenth perihelion encounter (E10) from 18 November 2021 until 21 November 2021 (see §4.4). We then adjust $\alpha _{\mathrm{H}}$ and $r_{\mathrm{H}}$ to improve slightly the agreement between the model and observational constraints on $T_{\mathrm{e}}$ from PSP during E10 (see §4.4).

Although it is our variation of the parameters $c_{\mathrm{refl}}$ , $c_{\mathrm{LD}}$ , $c_2$ , $\alpha _{\mathrm{H}}$ and $r_{\mathrm{H}}$ that enables our model to match the observations described in the previous paragraph, our model agrees with other observations without the need for further parameter adjustment. These observations include the latitudinal dependence of $U$ at large $r$ and the radial profiles of $n$ , $\delta v_{\mathrm{rms}}$ and $T_{\parallel \textrm {p}}$ seen by PSP during E10, as we discuss in §§4.3 and 4.4.

4.3. Comparison with Ulysses measurements of $U_\infty (\Theta )$ near solar minimum

Because the radiative cooling rate is $\propto n^2$ , radiative cooling is only important in our model in the low corona.Footnote ⁹ In steady state outside the low corona, both $a(r) F_{\mathrm{tot}}(r)$ and $a(r) \rho (r) U(r)$ are independent of $r$ , and thus, $F_{\mathrm{tot}}(r)/ [\rho (r) U(r)]$ is independent of $r$ , where $F_{\mathrm{tot}}$ is the total energy flux defined in (2.21). Asymptotically far from the Sun, almost all of the solar-wind energy flux is in the form of bulk-flow kinetic energy, and the wind speed is approximately

(4.3) \begin{equation} U_{\infty } = \left [\frac {2 F_{\mathrm{tot}}(r)}{\rho (r) U(r)}\right ]^{1/2}, \end{equation}

where the right-hand side of (4.3) can be evaluated at any $r$ far outside the low corona. In figure 3(a) we plot $U_\infty$ as a function of heliolatitude $\Theta$ in our model as well as measurements of $U$ from the first polar orbit of the Ulysses spacecraft in 1994–1995, near the minimum of the solar cycle. Notably, our model reproduces the observed decrease in the outflow velocity as $|\Theta |$ drops to ${\lesssim}20^\circ$ . The reason for this reduction is that, at $|\Theta | \lesssim 20^\circ$ , there is a large drop in $B$ as $r$ increases from ${\simeq}1.2 R_{\odot }$ to ${\simeq}1.6 R_{\odot }$ , as illustrated in figure 2. This drop in $B$ leads to a large Alfvén-speed gradient, which enhances the turbulent heating rata via (3.19). This enhanced heating rate increases the temperature, which increases the density scale height, causing more mass to be loaded into the corona out to the sonic critical point (at a heliocentric distance of a few $R_{\odot }$ ). This, in turn, increases the mass loss rate of the solar wind (Hansteen & Velli Reference Hansteen and Velli2012), causing the total energy budget to be shared by more particles, which reduces $U_\infty$ (cf. Leer & Holzer Reference Leer and Holzer1980; Chandran Reference Chandran2021).

Figure 3. (a) The asymptotic solar-wind speed $U_\infty$ defined in (4.3) as a function of heliolatitude $\Theta$ and the solar-wind speed measured by the Ulysses spacecraft during its first polar orbit in 1994 and 1995 (Goldstein et al. Reference Goldstein, Neugebauer, Phillips, Bame, Gosling, McComas, Wang, Sheeley and Suess1996). (b) The transverse pressure at $r= 68 R_{\odot }$ in our model as a function of $\Theta$ .

Cranmer et al. (Reference Cranmer, van Ballegooijen and Edgar2007) and Lionello et al. (Reference Lionello, Velli, Downs, Linker and Mikić2014a ) developed 1-D AW-driven solar-wind models and used the same Banaszkiewicz et al. (Reference Banaszkiewicz, Axford and McKenzie1998) magnetic-field model that we have used to model solar-wind streams at different heliolatitudes. In the study by Lionello et al. (Reference Lionello, Velli, Downs, Linker and Mikić2014a ), the drop in $U_\infty (\Theta )$ and increase in thermal pressure near the ecliptic were more sharply peaked at small $\Theta$ than in figure 3. In the study by Cranmer et al. (Reference Cranmer, van Ballegooijen and Edgar2007), the drop in $U_\infty (\Theta )$ and change in thermal pressure were likewise more sharply peaked at small $\Theta$ than in figure 3, but the thermal pressure decreased at small $\Theta$ . Lionello et al. (Reference Lionello, Velli, Downs, Linker and Mikić2014a ) pointed out that pressure differentials across the magnetic field lead to a loss of transverse force balance, and that, in a more realistic two-dimensional (2-D) or 3-D model, the high-pressure region would expand into the low-pressure region, thereby modifying the magnetic-field profiles. To obtain an order-of-magnitude estimate of the importance of this effect, we note that the increment in $\Theta$ through which a high-pressure region can expand during the time $\sim r/U$ it takes the wind to reach heliocentric distance $r$ is approximately $\Delta \Theta = (r/U) \times (\Delta c_{\mathrm{s}} / r) = \Delta c_{\mathrm{s}} / U$ radians, where $\Delta c_{\mathrm{s}}$ is the change in sound speed over the heliolatitude interval $\Delta \Theta$ . An upper limit on $\Delta \Theta$ is obtained by setting $\Delta c_{\mathrm{s}} = c_{\mathrm{s}}$ , in which case $\Delta \Theta = c_{\mathrm{s}}/U = 1/M$ , where $M$ is the sonic Mach number. As $M\gtrsim 10$ at large $r$ , this upper limit on $\Delta \Theta$ is $\lesssim 5.7^\circ$ . However, the pressure differential over $5.7^\circ$ in figure 3 is only a fraction of the full pressure, and so the actual value of $\Delta \Theta$ is significantly less than $5.7^\circ$ . We thus conjecture that the meridional expansion of the slow-wind region that would arise in a 2-D or 3-D version of our model would be comparatively modest.

4.4. Comparison with PSP E10 observations

During the approach phase of PSP’s tenth perihelion encounter with the Sun (E10), the spacecraft travelled approximately radially in the Sun’s rotating frame for ${\simeq}20 R_{\odot }$ , remaining within a single solar-wind stream that was magnetically connected to a small low-latitude coronal hole (Davis et al. Reference Davis, Chandran, Bowen, Badman, de Wit, Chen, Bale, Huang, Sioulas and Velli2023). This ‘fast radial scan’ provided a good opportunity to measure the radial variation of the solar-wind properties within the same magnetic flux tube, and we therefore compare our model results with PSP E10 data. A problem with doing so, however, is that the low-latitude coronal hole from which this wind stream emerged was a non-axisymmetric structure that is not captured by the axisymmetric B–H magnetic-field model. To compare our model with PSP E10 observations, we assume that the $B(r)$ profile along the magnetic flux tube encountered by PSP is similar to the $B(r)$ profile along the magnetic-field lines in our model along which $U_\infty \simeq 550 \,\mbox{ km} \mbox{ s}^{-1}$ (a value slightly above the proton outflow velocity measured by PSP during E10). We then compare PSP E10 observations with our model solution along one of these field lines: the one that reaches heliolatitude $\Theta = 20.9^\circ$ at large $r$ , for which $U_\infty = 541 \mbox{ km} \mbox{ s}^{-1}$ .

Figure 4. Radial profiles of several quantities in our model solution that reaches heliolatitude $\Theta = 20.9^\circ$ at large $r$ : $n$ is the proton number density, $U$ is the solar-wind outflow velocity, $v_{\mathrm{A}}$ is the Alfvén speed, $\delta v_{\mathrm{rms}}$ is the r.m.s. fluctuating velocity, $T_{\perp \textrm {p}}$ and $T_{\parallel \textrm {p}}$ are the perpendicular and parallel proton temperatures, $T_{\mathrm{e}}$ is the electron temperature, $\chi ^+_{L_\perp }$ is the outer-scale critical-balance parameter defined in (3.13), $\beta _{\parallel \textrm {p}}$ is the parallel proton beta and $\lambda _{\mathrm{mfp}}$ is the electron Coulomb mean free path. The vertical dotted line in panel (b) is described in the text.

In figure 4 we plot our numerical solutions for $n$ , $U$ , $T_{\perp \textrm {p}}$ , $T_{\parallel \textrm {p}}$ , $T_{\mathrm{e}}$ , $\beta _{\parallel \textrm {p}}$ , $\chi ^+_{L_\perp }$ and the ratio of the electron Coulomb mean free path

(4.4) \begin{equation} \lambda _{\mathrm{mfp}} \equiv \frac {\sqrt {k_{\mathrm{B}} T_{\mathrm{e}}/m_{\mathrm{e}}}}{\nu _{\mathrm{e}}} \end{equation}

to $r$ , where

(4.5) \begin{equation} \nu _{\mathrm{e}} = 2.9 \times 10^{-6} \left (\frac {n}{1 \mbox{ cm}^{-3}}\right ) \left (\frac {k_{\mathrm{B}}T_{\mathrm{e}}}{1 \mbox{ eV}}\right )^{-3/2} \ln \varLambda \mbox{ s}^{-1} \end{equation}

is the electron collision frequency and $\ln \varLambda$ is the Coulomb logarithm defined following (2.15) (Book Reference Book1983). The vertical dotted line in figure 4(b) indicates the radius inside of which the Alfvén-speed gradient is so large that $\delta z^-_{\mathrm{rms}}$ is limited to the value $c_{\mathrm{refl}} (U+v_{\mathrm{A}})/v_{\mathrm{A}}$ via (3.6) and (3.7). The PSP $n$ and $T_{\mathrm{e}}$ data in figure 4 are taken from the SPAN-E instrument (Whittlesey et al. Reference Whittlesey2020) of the Solar Wind Electrons, Alphas and Protons (SWEAP) instrument suite (Kasper et al. Reference Kasper2016), and the velocity and proton-temperature data are taken from SWEAP’s SPAN-I instrument (Livi et al. Reference Livi2022). The overall calibration of the SPAN-E density measurements is obtained by comparing with quasithermal-noise measurements of the electron density (Moncuquet et al. Reference Moncuquet2020) from the PSP FIELDS instrument suite (Bale et al. Reference Bale2016). We also plot in figure 4 remote observations of $T_{\perp \textrm {p}}$ from UVCS (Esser et al. Reference Esser, Fineschi, Dobrzycka, Habbal, Edgar, Raymond, Kohl and Guhathakurta1999), $T_{\mathrm{e}}$ from SUMER (Landi Reference Landi2008), $\delta v_{\mathrm{rms}}$ from Hinode (De Pontieu et al. Reference De Pontieu2007) and $n$ from white-light scattering in coronal holes (Allen Reference Allen1973). We note, however, that the UVCS, SUMER and white-light-scattering observations are for large polar coronal holes (the source regions of fast solar wind with $U_\infty \simeq 750 \mbox{ km} \mbox{ s}^{-1}$ ), which may differ from the small low-latitude coronal hole that produced the wind stream observed by PSP during E10.

Our model solution is in broad, but not complete, agreement with these observations. The most notable differences are that (1) $T_{\perp \textrm {p}}$ in the model is smaller than the (polar-coronal-hole) value observed by UVCS at $1.3 R_{\odot } \lt r \lt 1.8 R_{\odot }$ ; and (2) $T_{\mathrm{e}}$ in the model is higher than in the PSP observations, particularly near $r=15 R_{\odot }$ .

Beyond the comparison with observations, figure 4(f) shows that (1) $\delta \boldsymbol{z}^-$ fluctuations transition out of the strong-turbulence regime as $r$ increases past ${\simeq} 2 R_{\odot }$ ; (2) the low corona is highly collisional, with $\lambda _{\mathrm{mfp}} \ll r$ , whereas the outer corona and solar wind are only weakly collisional, with $\lambda _{\mathrm{mfp}} \gtrsim r$ ; and (3) $\beta _{\parallel \textrm {p}}$ is $\lesssim 0.1$ at $r\lt 20 R_{\odot }$ but increases with increasing $r$ , reaching a value ${\simeq} 1$ by $r=50 R_{\odot }$ .

Figure 5. (a) Heating and cascade rates divided by the mass density $\rho$ : $Q$ is the total heating rate; $Q_{\mathrm{e1}}$ ( $Q_{\parallel \textrm {p}}$ ) is the rate at which electron (proton) LD/TTD drain energy from fluctuations at $k_\perp \rho _{\mathrm{p}} \sim 1$ ; $\epsilon _{\textrm {sub-proton}}$ is the energy-cascade rate at $k_\perp \rho _{\mathrm{p}} \gt 1$ ; and $Q_{\perp \textrm {p}}$ is the proton stochastic heating rate. (b) The total energy flux $F_{\mathrm{tot}}(r)$ defined in (2.21) multiplied by the cross-sectional area of the flow $a(r)$ normalised to the value of this product at the inner radius of the model $r_{\mathrm{b}}$ .

Figure 5 shows several diagnostics of the heating and energy flow in this solution. Figure 5(a) shows that electron LD at $\lambda \sim \rho _{\mathrm{p}}$ is the dominant heating mechanism at $r\lesssim 2 R_{\odot }$ , that stochastic proton heating is the dominant heating mechanism at $r \gtrsim 2 R_{\odot }$ , and that dissipation at sub-proton scales and parallel proton heating are subdominant at all $r$ . The reason that parallel proton heating is negligible at $r\lesssim 15 R_{\odot }$ is that it requires a Landau resonance between the waves and the ions, in which the parallel ion velocity (along $\boldsymbol{B}$ ) matches the parallel wave phase velocity. However, the parallel thermal speed is $\beta _{\parallel \textrm {p}} ^{1/2} v_{\mathrm{A}}$ , and $\beta _{\parallel \textrm {p}} \lt 0.1$ at $r\lesssim 15 R_{\odot }$ , as shown in figure 4(f), implying that the large majority of the protons are much slower than the AWs and kinetic Alfvén waves (KAWs), whose parallel phase velocities are $\geqslant v_{\mathrm{A}}$ . The reason that the stochastic proton heating rate $Q_{\perp \textrm {p}}$ becomes relatively more important in comparison to the parallel electron heating rate as $r$ increases above ${\simeq} 2 R_{\odot }$ is that $Q_{\perp \textrm {p}}$ is a rapidly increasing function of $\delta v_{\mathrm{rms}}/v_{\mathrm{A}}$ because of the exponential in (3.52). As illustrated in figure 4(b), $\delta v_{\mathrm{rms}}/v_{\mathrm{A}}$ increases monotonically from $\lt 10^{-2}$ in the low corona to ${\simeq} 1/4$ at $r=10 R_{\odot }$ .

A comparison between figure 5(a) and figure 4(b) shows that the heating rate vanishes at the local extrema of $v_{\mathrm{A}}(r)$ . This feature of our model results from the fact that we have estimated $\delta z^-_{\mathrm{rms}}$ in (3.6) using only the local value of $|\partial v_{\mathrm{A}}/\partial r|$ , which determines the local rate of reflection. In coronal holes and the solar wind, after $\delta \boldsymbol{z}^-$ fluctuations are produced by reflection, they propagate some distance before cascading and dissipating, and $\delta z^-_{\mathrm{rms}}$ should not vanish anywhere. It would be possible to smooth out the $\delta z^-_{\mathrm{rms}}$ and $Q$ profiles, but we have opted not to do so, as the vanishing of $Q$ at these locations does not lead to sharp features in the temperature profiles and because smoothing would introduce an additional free parameter (the smoothing length).

In figure 5(b) we plot the product of the total energy flux $F_{\mathrm{tot}}(r)$ and the area of the flow $a(r)$ . This product should be constant in steady state in the absence of radiative cooling, as per (2.19). The fact that $a(r) F_{\mathrm{tot}}(r)$ is approximately constant at $r\gtrsim 1.2 R_{\odot }$ in figure 5 shows that radiative cooling plays a role only in the low solar corona. The sharp drop in $a(r) F_{\mathrm{tot}}(r)$ at the left of the plot represents a decrease in the total energy flux of ${\simeq} 3 \%$ over a distance of ${\simeq} 130 \mbox{ km}$ . Although this abrupt drop has the appearance of a possible artefact associated with the boundary condition and the neighbouring grid cell, the drop is in fact resolved across ${\simeq} 7$ grid cells and coincides with the large drop in the density (by a factor of ${\simeq} 2.5$ ) over this same radial interval. We interpret this abrupt drop as an approximation of the physical cooling that occurs near the top of the transition region.

Figure 6. (a) The quantity $P(q)$ is the PDF of the integer $q$ in the equation for the proton-gyroscale fluctuation amplitude $\delta z^+_{\mathrm{p}} = (0.7035)^q \delta z^+_{\mathrm{rms}}$ (see §3.4), where $\delta z^+_{\mathrm{p}}$ is the value of $\delta z^+_\lambda$ at $\lambda = \mathrm{\pi }\rho _{\mathrm{p}}$ . The top and bottom axes show $q$ and the corresponding value of $\delta z^+_{\mathrm{p}}/ \delta z^+_{\mathrm{rms}}$ . (b) The quantities  $P(q) Q_{\perp \mathrm{p},q}$ and $P(q) Q_{\mathrm{e1},q}$ show, respectively, the contribution to the perpendicular proton heating rate and electron heating rate from each part of the distribution of fluctuation amplitudes at $\lambda = \mathrm{\pi } \rho _{\mathrm{p}}$ . The quantity $P(q) (\epsilon ^+_q - Q_{\perp \mathrm{p},q} - Q_{\mathrm{e1},q})$ is the amount of the $\boldsymbol{z}^+$ cascade power at $\lambda = \mathrm{\pi } \rho _{\mathrm{p}}$ that escapes dissipation at $\lambda = \mathrm{\pi } \rho _{\mathrm{p}}$ and cascades to $\lambda \ll \rho _{\mathrm{p}}$ as a function of $q$ . Both panels correspond to $r=2R_{\odot }$ along the model magnetic-field line that reaches heliolatitude $\Theta = 20.9^\circ$ at large $r$ .

Figure 6(a) shows the PDF of the integer $q$ in (3.31), evaulated for $\delta z^+_\lambda$ at $\lambda = \pi \rho_{\rm p}$ (i.e., the PDF of $\log_\beta(\delta z^+_{\rm p}/\overline{z}^+)$ ). Figure 6(b) shows the contributions of the different parts of the $\delta z^+_{\mathrm{p}}$ distribution to $Q_{\perp \textrm {p}}$ and $Q_{\mathrm{e1}}$ (which are, respectively, the stochastic proton heating rate and the rate at which electrons are heated by LD/TTD of Alfvénic fluctuations at the proton gyroscale). As these panels show, stochastic proton heating results mainly from the large-amplitude end of the distribution, whereas $Q_{\mathrm{e1}}$ is spread out more evenly across the $\delta z^+_{\mathrm{p}}$ distribution. Figure 6(b) also shows that $Q_{\perp \mathrm{p},q} + Q_{\mathrm{e}1, q} \simeq \epsilon ^+_q$ at $q\leqslant 7$ , which indicates that almost all of the cascade power going into the strongest fluctuations at $\lambda \simeq \rho _{\mathrm{p}}$ is dissipated by stochastic heating and electron LD. Damping at $\lambda \sim \rho _{\mathrm{p}}$ thus truncates the tail of the $\delta z^+_{\mathrm{p}}$ distribution, reducing the degree of intermittency at $\lambda \lesssim \rho _{\mathrm{p}}$ .

4.5. Energy flux ratios

In figure 7 we plot the individual components of the total energy flux $F_{\mathrm{tot}}$ divided by $F_{\mathrm{tot}}$ , where $F_{\mathrm{tot}}$ is defined in (2.21). To label these components, we have rewritten (2.21) as

(4.6) \begin{equation} F_{\mathrm{tot}} = F_U + F_{\mathrm{g}} + F_{\mathrm{enth,e}} + F_{\mathrm{enth,p}} + q_{\mathrm{e}} + q_{\mathrm{p}} + F_{\mathrm{w}}, \end{equation}

where

(4.7) \begin{equation} F_U = \frac { \rho U^3}{2} \end{equation}

is the bulk-flow kinetic energy flux,

(4.8) \begin{equation} F_{\mathrm{g}} = - \frac { G M_{\odot }\rho U}{r} \end{equation}

is the gravitational potential-energy flux,

(4.9) \begin{equation} F_{\mathrm{enth,e}} = \frac {5}{2}n k_{\mathrm{B}} T_{\mathrm{e}}U \end{equation}

is the electron enthalpy flux,

(4.10) \begin{equation} F_{\mathrm{enth,p}} = n k_{\mathrm{B}} \left ( T_{\perp \textrm {p}} + \frac { 3}{2}T_{\parallel \textrm {p}}\right )U \end{equation}

is the proton enthalpy flux, $q_{\mathrm{e}}$ is the electron heat flux,

(4.11) \begin{equation} q_{\mathrm{p}} = q_{\perp \textrm {p}} + q_{\parallel \textrm {p}} \end{equation}

is the proton heat flux and

(4.12) \begin{equation} F_{\mathrm{w}} = \left ( \frac {3 }{2}U + v_{\mathrm{A}} \right ){\mathcal E}_{\mathrm{w}} \end{equation}

is the AW enthalpy flux.

Figure 7. The individual energy fluxes defined in (4.7) through (4.12) expressed as fractions of the total energy flux $F_{\mathrm{tot}}$ defined in (4.6). Each panel corresponds to a different heliolatitude $\Theta$ , asymptotic wind speed $U_{\infty }$ and energy per proton $E$ .

As mentioned in §2.3, outside the low solar corona, radiation can be neglected, and hence, (2.19) implies that

(4.13) \begin{equation} a(r)F_{\mathrm{total}}(r) = \mbox{ constant} \end{equation}

in steady state. As the solar wind flows from the corona out into the distant interplanetary medium, power $a(r) F_{\mathrm{total}}(r)$ is converted from one form to another without loss. At large $r$ , $F_{\mathrm{tot}}$ is dominated by $F_U$ , but in the corona $F_U$ is negligible, and thus, $F_U/F_{\mathrm{tot}}$ increases with $r$ . Likewise, $F_{\mathrm{g}}$ is large and negative in the corona but negligible at large $r$ , and thus, $F_{\mathrm{g}}/F_{\mathrm{tot}}$ increases with $r$ . As $(F_U+F_{\mathrm{g}})/F_{\mathrm{tot}}$ increases with $r$ , other fractional energy fluxes, such as $F_{\mathrm{w}}/F_{\mathrm{tot}}$ , must decrease to maintain total energy conservation, and the fractional energy flux that decreases the most corresponds to the mechanism that is doing the most to accelerate the wind and lift it out of the Sun’s gravitational potential well.

Figure 7(a) illustrates a slow-wind stream with asymptotic wind speed $U_{\infty } = 306 \mbox{ km} \mbox{ s}^{-1}$ , where $U_\infty$ is defined in (4.3). This value of $U_\infty$ corresponds to an energy per proton

(4.14) \begin{equation} E = \frac {1}{2} m_{\mathrm{p}} U_\infty ^2 \end{equation}

of 487 eV, as indicated in the figure legend. As can be seen in the figure, the primary mechanism that powers the wind is the AW enthalpy flux at $r\lesssim 1.5 R_{\odot }$ , the electron heat flux between $\sim 1.5R_{\odot }$ and $\sim 4 R_{\odot }$ , and the electron enthalpy flux at $r\gtrsim 4 R_{\odot }$ .

Figure 7(b,c,d) illustrates wind streams with larger values of $U_\infty$ and $E$ . In these plots the AW enthalpy flux is the dominant driver of the wind out to $r\sim 20{-}30 R_{\odot }$ , beyond which some form of enthalpy flux becomes dominant: either a roughly equal mix of $F_{\mathrm{enth,p}}$ and $F_{\mathrm{enth,\!e}}$ in figure 7(b) or mostly $F_{\mathrm{p}}$ in figures 7(c) and 7(d).

Halekas et al. (Reference Halekas2023) constructed similar plots using PSP data, but included the electric-potential-energy flux $F_{\Phi _E} = e\Phi _E n U$ , which does not appear in our total energy equation. The significance of $F_{\Phi _E}$ can be understood by considering single-species versions of the total energy equation, as in equation (23) of Lemaire & Scherer (Reference Lemaire and Scherer1973),

(4.15) \begin{align} & a(r)\left \{ q_{\mathrm{e}} + nU\left [ \frac {1}{2} m_{\mathrm{e}} U^2 + k_{\mathrm{B}} \left (\frac {3}{2} T_{\parallel \textrm e} + T_{\perp \textrm e}\right ) + m_{\mathrm{e}} \Phi _{\mathrm{g}} - e \Phi _E\right ]\right \}\nonumber\\&\quad{} = a(r) F_{\mathrm{tot, e}} = \mbox{ constant}, \end{align}

(4.16) \begin{align} & a(r)\left \{ q_{\mathrm{p}} + nU\left [ \frac {1}{2} m_{\mathrm{p}} U^2 + k_{\mathrm{B}} \left (\frac {3}{2} T_{\parallel \textrm {p}} + T_{\perp \textrm {p}}\right ) + m_{\mathrm{p}} \Phi _{\mathrm{g}} + e \Phi _E\right ]\right \}\nonumber\\&\quad{} =a(r) F_{\mathrm{tot, p}} = \mbox{constant}, \\[6pt] \nonumber \end{align}

where $T_{\parallel \textrm e}$ and $T_{\perp \textrm e}$ are the parallel and perpendicular electron temperatures, and $\Phi _g = - GM_{\odot }/r$ is the gravitational potential. Equations (4.15) and (4.16) do not account for AWs, a point that will come up again following (4.17). When (4.15) and (4.16) are added, the $\Phi _E$ terms cancel, which is why $\Phi _E$ does not appear in our (4.6) for the total energy flux. Fundamentally, this is because the electric field does no work on a neutral proton--electron plasma, although it does do work on the protons and equal and opposite work on the electrons. Far from the Sun, almost all of the solar-wind energy flux is in the form of the bulk-flow kinetic energy of the protons, which implies that $F_{\mathrm{tot,e}} \ll F_{\mathrm{tot,p}}$ at all $r$ . As the proton bulk-flow kinetic energy flux and gravitational potential-energy flux are at most comparable to $F_{\mathrm{tot,p}}$ , the electron bulk-flow kinetic energy flux and gravitational potential-energy flux are everywhere at least a factor of ${\simeq} m_{\mathrm{e}}/m_{\mathrm{p}}$ smaller than $F_{\mathrm{tot, p}}$ , and hence, negligible. After dividing (4.15) by $a(r)$ , we are left with three terms that are everywhere ${\ll} F_{\mathrm{tot,p}}$ , namely $F_{\mathrm{tot,e}}$ , $n m_{\mathrm{e}} U^3/2$ and $n m_{\mathrm{e}} U \Phi _g$ , and three terms that can, at least at some locations, be of the same order of magnitude as $F_{\mathrm{tot,p}}$ : $q_{\mathrm{e}}$ , $n U e \Phi _E$ and $5n k_{\mathrm{B}} T_{\mathrm{e}}/2$ , where we have set $T_{\parallel \textrm e} = T_{\perp \textrm e} = T_{\mathrm{e}}$ . After dropping the three negligible terms, we can rewrite (4.15) as

(4.17) \begin{equation} F_{\Phi _E} \equiv n U e \Phi _E = q_{\mathrm{e}} + \frac {5}{2} n U k_{\mathrm{B}} T_{\mathrm{e}}. \end{equation}

In other words, the electric-potential-energy flux $F_{\Phi _E}$ is approximately the sum of the electron heat flux and electron enthalpy flux. The AW enthalpy flux can lead to violations of (4.17) if the AWs dissipate and heat the electrons. However, we expect (4.17) to be accurate when $F_{\Phi _E}$ is large compared with $F_{\mathrm{w}}$ and of the same order of magnitude as $F_{\mathrm{tot,p}}$ .

The PSP measurements reported in figure 3 of Halekas et al. (Reference Halekas2023) show that in the fastest wind streams observed in the study, in which $2000 \mbox{ eV} \lt E \lt 2500 \mbox{ eV}$ , AWs were the dominant energisation mechanism over the comparatively narrow radial interval $13 R_{\odot } \lt r \lt 23 R_{\odot }$ . On the other hand, in the slowest wind streams with $E\lt 500 \mbox{ eV}$ , $F_{\Phi _E}$ was the dominant energisation mechanism at $13 R_{\odot } \lt r \lt 48 R_{\odot }$ , with $q_{\mathrm{e}}$ significantly smaller than $F_{\Phi _E}$ . This latter point implies via (4.17) that $F_{\mathrm{enth,\!e}}$ was the dominant energisation mechanism in these low-speed streams at these radii. Taken on their own, these observations are consistent with a long-standing conjecture that the fast solar wind is driven primarily by an AW energy flux, while the slow solar wind is driven by some other mechanism.

However, our results, which provide a model for what is happening interior to the orbit of PSP, paint a different picture. Although the energy flux ratios shown in figure 7 agree with the results from Halekas et al. (Reference Halekas2023) described in the previous paragraph, in our model $F_{\mathrm{w}}$ is the dominant outward (i.e. positive) energy flux close to the Sun for all values of $E$ . As mentioned previously, for the wind stream modelled in figure 7(a), the dominant energy flux powering the wind is $F_{\mathrm{w}}$ at $r\lt 1.5R_{\odot }$ , $q_{\mathrm{e}}$ between $1.5 R_{\odot }$ and $4 R_{\odot }$ , and $F_{\mathrm{enth,\!e}}$ at $r\gtrsim 4 R_{\odot }$ . Our results show that the observation that $F_{\mathrm{enth,\!e}}$ or $q_{\mathrm{e}}$ dominates at large $r$ does not mean that one of these mechanisms dominates near the Sun. Indeed, it would not make sense for $F_{\mathrm{enth,\!e}}$ to be the dominant mechanism powering the solar wind between the coronal base and $r\gtrsim 10 R_{\odot }$ , for if it were then the divergence of the electron enthalpy flux would be large compared with both the turbulent heating rate and the divergence of the electron heat flux. This would cause the electrons to cool approximately adiabatically over this range of radii, which is ruled out by observational constraints on $T_{\mathrm{e}}$ . Nor would it make sense for $q_{\mathrm{e}}$ to be the dominant energy flux powering the solar wind at $r\lesssim 1.5 R_{\odot }$ , because in this region $\lambda _{\mathrm{mfp}} T_{\mathrm{e}}^{-1} \mathrm{d} T_{\mathrm{e}}/\mathrm{d}r \ll 1$ and $\mathrm{d}T_{\mathrm{e}}/\mathrm{d} r \gt 0$ (see figures 4 and 8), implying that the electron heat flux is well approximated by the Spitzer--Harm formula (2.22) and directed in towards the Sun. In other words, in the low corona, $q_{\mathrm{e}}$ is an energy sink, not an energy source, and therefore, cannot be the dominant energy flux that powers the solar wind. Instead, some non-thermal energy source, such as AWs, is required to power the corona and solar wind (cf. Parker Reference Parker1965).

Figure 8. Our model solution that reaches heliolatitude $\Theta =1.1^\circ$ at large $r$ . The various quantities plotted are defined in the captions to figures 4 and 5 and the text following (4.5).

Figure 8 provides more detailed information about the solar-wind stream illustrated in figure 7(a), which reaches a heliolatitude of $\Theta = 1.1^\circ$ at large $r$ . As illustrated in figure 2(b), in this wind stream $|B|$ drops rapidly between $r=1.2 R_{\odot }$ and $r=1.6 R_{\odot }$ , causing $v_{\mathrm{A}}$ to drop by a factor of almost 10 over this same range of radii, as illustrated in figure 8(b). The resulting large Alfvén-speed gradient leads to strong AW reflection and intense AW dissipation and turbulent heating via (3.19), most of which goes into electron heating because $\delta v_{\mathrm{rms}}/v_{\mathrm{A}}$ is still small enough at these radii (figure 8 b) that stochastic ion heating is weaker than electron heating, as shown in figure 8(d). At $1.2 R_{\odot } \lt r \lesssim 2 R_{\odot }$ , this strong electron heating causes $T_{\mathrm{e}}$ in this slow-wind stream to exceed the value of $T_{\mathrm{e}}$ in the moderately fast solar-wind stream illustrated in figure 4. As the Spitzer–Harm heat flux in (2.22) is $\propto T_{\mathrm{e}}^{7/2}$ , the increase in $T_{\mathrm{e}}$ leads to a large increase in $q_{\mathrm{e}}$ , which causes $q_{\mathrm{e}}$ to become the dominant outward energy flux at $1.5 R_{\odot } \lesssim r \lesssim 4 R_{\odot }$ , as shown in figure 7. The reason that the intense heating inside $2R_{\odot }$ decreases $U_{\infty }$ relative to the model solution shown in figure 4 is explained following (4.3). Finally, we note that radiative losses are somewhat larger in figure 8(f) than in figure 5(b).