An approximate analytic solution to the coupled problems of coronal heating and solar-wind acceleration

Abstract

Between the base of the solar corona at $r=r_\textrm {b}$ and the Alfvén critical point at $r=r_\textrm {A}$, where $r$ is heliocentric distance, the solar-wind density decreases by a factor $ \mathop > \limits_\sim 10^5$, but the plasma temperature varies by a factor of only a few. In this paper, I show that such quasi-isothermal evolution out to $r=r_\textrm {A}$ is a generic property of outflows powered by reflection-driven Alfvén-wave (AW) turbulence, in which outward-propagating AWs partially reflect, and counter-propagating AWs interact to produce a cascade of fluctuation energy to small scales, which leads to turbulent heating. Approximating the sub-Alfvénic region as isothermal, I first present a brief, simplified calculation showing that in a solar or stellar wind powered by AW turbulence with minimal conductive losses, $\dot {M} \simeq P_\textrm {AW}(r_\textrm {b})/v_\textrm {esc}^2$, $U_{\infty } \simeq v_\textrm {esc}$, and $T\simeq m_\textrm {p} v_\textrm {esc}^2/[8 k_\textrm {B} \ln (v_\textrm {esc}/\delta v_\textrm {b})]$, where $\dot {M}$ is the mass outflow rate, $U_{\infty }$ is the asymptotic wind speed, $T$ is the coronal temperature, $v_\textrm {esc}$ is the escape velocity of the Sun, $\delta v_\textrm {b}$ is the fluctuating velocity at $r_\textrm {b}$, $P_\textrm {AW}$ is the power carried by outward-propagating AWs, $k_\textrm {B}$ is the Boltzmann constant, and $m_\textrm {p}$ is the proton mass. I then develop a more detailed model of the transition region, corona, and solar wind that accounts for the heat flux $q_\textrm {b}$ from the coronal base into the transition region and momentum deposition by AWs. I solve analytically for $q_\textrm {b}$ by balancing conductive heating against internal-energy losses from radiation, $p\,\textrm {d} V$ work, and advection within the transition region. The density at $r_\textrm {b}$ is determined by balancing turbulent heating and radiative cooling at $r_\textrm {b}$. I solve the equations of the model analytically in two different parameter regimes. In one of these regimes, the leading-order analytic solution reproduces the results of the aforementioned simplified calculation of $\dot {M}$, $U_\infty$, and $T$. Analytic and numerical solutions to the model equations match a number of observations.

1. Introduction

Pioneering works by Parker (1958, 1965), Hartle & Sturrock (1968), and Durney (1972) modelled the solar wind as a steady-state, spherical outflow powered by the outward conduction of heat from the base of the corona. These models succeeded in producing a supersonic wind, but were unable to explain the large outflow velocities measured in fast-solar-wind streams near Earth. They also had little predictive power for the mass outflow rate from the Sun, $\dot {M}$, because they specified the temperature of the coronal base as a boundary condition, and $\dot {M}$ is highly sensitive to the coronal temperature (Hansteen & Leer 1995).

A possible solution to these problems was proposed almost as soon as the difficulties became apparent, namely that the solar wind is powered by an Alfvén-wave (AW) energy flux (Parker 1965, p. 686; Hollweg 1973, 1978). This idea received strong support from the discovery of large-amplitude AWs in the interplanetary medium that propagate away from the Sun in the local plasma frame (Belcher & Davis 1971) as well as the remote observation of AW-like motions in the low solar atmosphere carrying an energy flux sufficient to power the solar wind (De Pontieu et al. 2007). A leading paradigm for how AWs energise the solar wind is based on reflection-driven AW turbulence. As AWs propagate away from the Sun, they undergo partial reflection because of the radial variation in the Alfvén speed (Heinemann & Olbert 1980; Velli 1993). Counter-propagating AWs then interact nonlinearly (Iroshnikov 1963; Kraichnan 1965), causing wave energy to cascade from large wavelengths to small wavelengths and dissipate, thereby heating the ambient plasma (Velli, Grappin & Mangeney 1989; Matthaeus et al. 1999; Cranmer & van Ballegooijen 2005; Verdini & Velli 2007; Perez & Chandran 2013; van Ballegooijen & Asgari-Targhi 2016, 2017). The presence of turbulence in the interplanetary medium is confirmed by spacecraft measurements (Coleman 1968; Matthaeus & Goldstein 1982; Tu & Marsch 1995; Smith et al. 2001; Bruno & Carbone 2005; Horbury, Forman & Oughton 2008; Wicks et al. 2010; Chen et al. 2020), and numerical solar-wind models based on reflection-driven AW turbulence have proven successful at explaining a number of observations of the solar wind and corona (Cranmer, van Ballegooijen & Edgar 2007; Verdini et al. 2010; Chandran et al. 2011; van der Holst et al. 2014; Usmanov, Goldstein & Matthaeus 2014). Three-dimensional compressible magnetohydrodynamic (MHD) simulations have also shown that the solar wind can be self-consistently generated by an AW energy flux (Shoda et al. 2019).

Although the aforementioned numerical models reproduce many of the observed properties of the solar wind (see also Riley et al. 2011; Gressl et al. 2014), analytic formulae that determine $\dot {M}$ and the outflow speed far from the Sun ($U_\infty$) remain elusive. A number of studies have obtained a single equation that constrains the two unknowns $\dot {M}$ and $U_{\infty }$. For example, Sandbaek, Leer & Hansteen (1994) pointed out that if the energy flux far from the Sun is mostly in the form of bulk-flow kinetic energy, and if the energy flux at the coronal base is dominated by the flux of gravitational potential energy, heat, and some additional form of mechanical energy (from, e.g. AWs), then energy conservation implies that

(1.1)\begin{equation} \dot{M} = \frac{\dot{E}_\textrm{m0} + \dot{E}_{q0}}{\frac{1}{2}(v_\textrm{esc}^2 + U_\infty^2)}, \end{equation}

where $\dot {E}_\textrm {m0}$ is the mechanical-energy input into the solar wind at the corona base, $\dot {E}_{q0}$ is the energy input at the coronal base from thermal conduction (which is negative in models that include the lower solar atmosphere),

(1.2)\begin{equation} v_\textrm{esc} = \left(\frac{2 G M_{{\odot}}}{R_{{\odot}}}\right)^{1/2} = 617.7 \mbox{ km} \mbox{ s}^{{-}1} \end{equation}

is the escape velocity of the Sun, $G$ is the gravitational constant, $M_\odot$ is the solar mass, and $R_{\odot }$ is the solar radius. Schwadron & McComas (2003) derived a variant of (1.1) that explicitly relates $\dot {E}_{q0}$ to the altitude of the coronal-temperature maximum. Hansteen & Leer (1995) and Hansteen, Leer & Holzer (1997) showed that $\dot {E}_{q0} \ll \dot {E}_\textrm {m0}$, and Hansteen & Velli (2012) made use of this finding to further refine (1.1), obtaining

(1.3)\begin{equation} \dot{M} = \frac{\dot{E}_\textrm{m0}}{\frac{1}{2}(v_\textrm{esc}^2 + U_\infty^2)}. \end{equation}

Another important result was obtained by Leer & Holzer (1980), who found that heating inside the sonic critical point enhances $\dot {M}$ but has little effect on $U_{\infty }$, whereas heating beyond the sonic critical point increases $U_{\infty }$ but has little effect on $\dot {M}$.

The main goal of the present paper is to obtain approximate analytic solutions for $\dot {M}$, $U_{\infty }$, the temperature of the corona, the heat flux from the coronal base into the lower solar atmosphere, and the plasma density at the coronal base under the assumption that AW turbulence is the primary energisation mechanism of the solar wind. Section 2 takes a first step towards this goal by presenting a simplified approximate calculation of $\dot {M}$, $U_{\infty }$, and the coronal temperature. Section 3 develops a more detailed solar-wind model that accounts for physical processes neglected in § 2. Approximate analytic solutions to the equations of this model are presented in § 3, and numerical solutions are presented in § 4. Section 5 discusses and summarises the main results of the paper.

2. Heuristic calculation of $\dot {M}$, $U_{\infty }$, and the coronal temperature in AW-driven winds with minimal conductive losses

Approximate expressions for $\dot {M}$, $U_\infty$, and the coronal temperature can be quickly obtained by modelling the solar wind as a spherically symmetric, steady-state outflow and assuming that: (1) AW turbulence is the dominant heating mechanism; (2) solar rotation can be neglected, so that the magnetic field $\boldsymbol {B}$ and flow velocity are aligned; (3) $\boldsymbol {B}\propto r^{-2}\boldsymbol {\hat {r}}$, where $\boldsymbol {\hat {r}}$ is the radial unit vector (i.e. a split monopole, with $B_r>0$ in one hemisphere and $B_r<0$ in the other); (4) momentum deposition by AWs can be neglected between the coronal base and sonic critical point; and (5) $p\,\textrm {d} V$ work is the dominant sink of internal energy in the sub-Alfvénic region of the solar wind, in which the solar-wind outflow velocity $U$ is smaller than the Alfvén speed

(2.1)\begin{equation} v_\textrm{A} = \frac{B}{\sqrt{4{\rm \pi}\rho}}, \end{equation}

where $\rho$ is the mass density. Assumptions (3)–(5) are relaxed in the next section.

In steady state, given assumption (2) above, mass and flux conservation imply (see § 3.1) that

(2.2)\begin{equation} v_\textrm{Ab} = y_\textrm{b} U_\textrm{b}, \end{equation}

where

(2.3)\begin{equation} y_\textrm{b} \equiv \left[\frac{\rho_\textrm{b}}{\rho(r_\textrm{A})}\right]^{1/2}, \end{equation}

$\rho _\textrm {b} = \rho (r_\textrm {b})$, $r_\textrm {b}$ is the radius of the coronal base, which in this section (but not the next) is simply set equal to $R_{\odot }$, $r_\textrm {A}$ is the Alfvén critical point at which $U=v_\textrm {A}$, $v_\textrm {Ab} = v_\textrm {A}(r_\textrm {b})$, and $U_\textrm {b} = U(r_\textrm {b})$. The mass outflow rate can thus be written in the form

(2.4)\begin{equation} \dot{M} = 4{\rm \pi} R_{{\odot}}^2 \rho_\textrm{b} U_\textrm{b} = 4{\rm \pi} R_{{\odot}}^2 \rho_\textrm{b} v_\textrm{Ab} y_\textrm{b}^{{-}1}. \end{equation}

As shown in § 3.2, heating by reflection-driven AW turbulence causes the sub-Alfvénic part of the solar wind at $r< r_\textrm {A}$ to become quasi-isothermal, in the sense that

(2.5)\begin{equation} \frac{1}{c_\textrm{s}^2} \left|\frac{\mathrm{d} c_\textrm{s}^2}{\mathrm{d}r}\right| \ll \frac{1}{\rho} \left| \frac{\mathrm{d}\rho}{\mathrm{d}r} \right|, \end{equation}

where

(2.6)\begin{equation} c_\textrm{s}^2 \equiv \frac{p}{\rho} = \frac{2 k_\textrm{B} T}{m_\textrm{p}} \end{equation}

is the square of the isothermal sound speed, $k_\textrm {B}$ is the Boltzmann constant, $T$ is the temperature, and $m_\textrm {p}$ is the proton mass. This argument is consistent with observations and models of the solar wind, which suggest that the temperature varies by a factor of only a few between $r_\textrm {b}$ and $r_\textrm {A}$, whereas $\rho$ varies by a factor of approximately $10^5$ (see, e.g. Cranmer et al. 2007). When the sub-Alfvénic solar wind is approximated as an isothermal plasma, the momentum equation at $r<r_\textrm {A}$ can be written in the form

(2.7)\begin{equation} U \frac{\mathrm{d}U}{\mathrm{d}r} ={-}\frac{c_\textrm{s}^2}{\rho} \frac{\mathrm{d}\rho}{\mathrm{d}r} - \frac{v_\textrm{esc}^2 R_{{\odot}}}{2 r^2}. \end{equation}

As $\dot {M} = 4{\rm \pi} r^2 U \rho$ is independent of $r$, $U^{-1} \,\textrm {d} U/\textrm {d} r = - \rho ^{-1} \,\textrm {d} \rho / \textrm {d} r - 2/r$. Upon substituting this relation into (2.7) and rearranging terms, one obtains

(2.8)\begin{equation} \frac{(c_\textrm{s}^2 - U^2)}{\rho}\frac{\mathrm{d}\rho}{\mathrm{d}r} = \frac{2 U^2}{r} - \frac{v_\textrm{esc}^2 R_{{\odot}}}{2r^2}. \end{equation}

In order for (2.8) to have a smooth transonic solution, the right-hand side of (2.8) must vanish at the radius $r_\textrm {c}$ at which $U=c_\textrm {s}$, so that $\textrm {d} \rho /\textrm {d} r$ remains finite. This leads to the two critical-point conditions

(2.9)\begin{equation} U_\textrm{c} = c_\textrm{s} \qquad \frac{r_\textrm{c}}{R_{{\odot}}} = \frac{v_\textrm{esc}^2}{4 c_\textrm{s}^2}, \end{equation}

where $U_\textrm {c} = U(r_\textrm {c})$. Integrating (2.7), one obtains the Bernoulli integral

(2.10)\begin{equation} \frac{1}{2} U^2 + c_\textrm{s}^2 \ln\left(\frac{\rho}{\rho_\textrm{b}}\right) - \frac{v_\textrm{esc}^2 R_{{\odot}}}{2 r} = \mbox{constant} ={-} \frac{v_\textrm{esc}^2}{2}, \end{equation}

where the second equality in (2.10) results from evaluating the left-hand side of (2.10) at $r=r_\textrm {b}$ and dropping the $U_\textrm {b}^2/2$ term, which is $\ll v_\textrm {esc}^2/2$. After evaluating the left-hand side of (2.10) at $r=r_\textrm {c}$ and using (2.9) to rewrite $U_\textrm {c}$ and $r_\textrm {c}$ in terms of $c_\textrm {s}^2$, one obtains

(2.11)\begin{equation} \ln\left( \frac{ \rho_\textrm{c}}{\rho_\textrm{b}}\right) ={-} \frac{v_\textrm{esc}^2}{2 c_\textrm{s}^2} + \frac{3}{2}, \end{equation}

where $\rho _\textrm {c} = \rho (r_\textrm {c})$. Upon setting $\dot {M} = 4{\rm \pi} r_\textrm {c}^2 \rho _\textrm {c} U_\textrm {c}$ and using (2.9) and (2.11) to rewrite $r_\textrm {c}$, $\rho _\textrm {c}$, and $U_\textrm {c}$ in terms of $c_\textrm {s}$, one obtains (Hansteen & Velli 2012)

(2.12)\begin{equation} \dot{M} = \frac{{\rm \pi} R_{{\odot}}^2 v_\textrm{esc}^4 \rho_\textrm{b}}{4 c_\textrm{s}^3} \exp\left( - \frac{v_\textrm{esc}^2}{2 c_\textrm{s}^2} + \frac{3}{2}\right). \end{equation}

The exponential appearing on the right-hand side of (2.12) reflects the fact that, at $r< r_\textrm {c}$, the flow is subsonic and the density drops off approximately as in a static atmosphere (Hansteen & Velli 2012). Equating the right-hand sides of (2.4) and (2.12), one finds that

(2.13)\begin{equation} \ln y_\textrm{b} = \frac{v_\textrm{esc}^2}{2 c_\textrm{s}^2} - \ln c_1 \simeq \frac{v_\textrm{esc}^2}{2 c_\textrm{s}^2}, \end{equation}

where $c_1 \equiv e^{3/2} v_\textrm {esc}^4/(16 c_\textrm {s}^3 v_\textrm {Ab})$.

When radiative cooling and thermal conduction are neglected, the internal-energy equation becomes

(2.14)\begin{equation} -c_\textrm{s}^2 U \frac{\mathrm{d}\rho}{\mathrm{d}r} + \frac{\rho U}{\gamma-1} \frac{\mathrm{d}c_\textrm{s}^2}{\mathrm{d}r} = Q, \end{equation}

where $Q$ is the turbulent heating rate, and $\gamma$ is the ratio of specific heats. In the quasi-isothermal approximation, the second term on the left-hand side of (2.14) can be neglected. Multiplying (2.14) by $4{\rm \pi} r^2$ and integrating over the quasi-isothermal sub-Alfvénic region, one obtains

(2.15)\begin{equation} -c_\textrm{s}^2 \dot{M} \int_{r_\textrm{b}}^{r_\textrm{A}} \frac{1}{\rho} \frac{\mathrm{d}\rho}{\mathrm{d}r} \textrm{d} r = 4{\rm \pi} \int_{r_{\textrm{b}}}^{r_\textrm{A}} r^2 Q(r) \, \textrm{d} r \simeq P_\textrm{AW}(r_\textrm{b}), \end{equation}

where the approximate equality assumes that the volume-integrated turbulent heating rate between $r_\textrm {b}$ and $r_\textrm {A}$ is comparable to $P_\textrm {AW}(r_\textrm {b})$, the power carried by outward-propagating AWs at the coronal base, consistent with direct numerical simulations (Perez et al. 2021). As $\ln (\rho _\textrm {b}/\rho (r_\textrm {A})) = 2 \ln y_\textrm {b}$, (2.15) yields

(2.16)\begin{equation} \dot{M} \simeq \frac{P_\textrm{AW}(r_\textrm{b})}{2 c_\textrm{s}^2 \ln y_\textrm{b}} \simeq \frac{P_\textrm{AW}(r_\textrm{b})}{v_\textrm{esc}^2}, \end{equation}

where the second relation in (2.16) follows from (2.13).

Equations (2.15) and (2.16) can be understood as follows. Within the quasi-isothermal sub-Alfvénic region, each time the density of a parcel of plasma decreases by a factor of $e$, its thermal energy must be replaced via heating to offset internal-energy losses from $p\,\textrm {d} V$ work. The quantity $c_\textrm {s}^2 \ln (\rho _\textrm {b}/\rho (r_\textrm {A})) = 2 c_\textrm {s}^2 \ln y_\textrm {b} \simeq v_\textrm {esc}^2$ is thus the heating cost per unit mass for plasma to transit the sub-Alfvénic region. The reason that the product $c_\textrm {s}^2 \ln (\rho _\textrm {b}/\rho (r_\textrm {A}))$ is approximately constant is that increasing $c_\textrm {s}^2$ leads to an exponential increase in $\dot {M}$ and the solar-wind density and an exponential reduction of $\rho _\textrm {b}/\rho (r_\textrm {A})$, leaving $c_\textrm {s}^2 \ln (\rho _\textrm {b}/\rho (r_\textrm {A}))$ approximately unchanged. Equations (2.15) and (2.16) state that $\dot {M}$ is the net heating power within the sub-Alfvénic region (which is taken to be $\simeq P_\textrm {AW}(r_\textrm {b})$) divided by the heating cost per unit mass.

Assuming that the AW-energy flux and gravitational-potential-energy flux are the dominant mechanical-energy fluxes at the coronal base and that the kinetic-energy flux is the dominant mechanical-energy flux at large $r$, and equating the mechanical luminosities at $r=r_\textrm {b}$ and at large $r$, one obtains

(2.17)\begin{equation} P_\textrm{AW}(r_\textrm{b}) - \tfrac{1}{2}\dot{M} v_\textrm{esc}^2 = \tfrac{1}{2} \dot{M} U_{\infty}^2. \end{equation}

Substituting (2.16) into (2.17) yields

(2.18)\begin{equation} U_{\infty} \simeq v_\textrm{esc}. \end{equation}

Reflection-driven AW turbulence thus changes the energy per unit mass from $\simeq - v_\textrm {esc}^2/2$ at the coronal base to $\simeq v_\textrm {esc}^2/2$ far from the Sun.

In the absence of super-radial expansion, $P_\textrm {AW}(r_\textrm {b}) = 4 {\rm \pi}R_{\odot }^2 v_\textrm {Ab} \rho _\textrm {b} (\delta v_\textrm {b})^2$, where $\delta v_\textrm {b}$ is the root-mean-square (r.m.s.) amplitude of the fluctuating velocity at the coronal base. This relation, in conjunction with (2.4) and (2.16), implies that $y_\textrm {b} \simeq v_\textrm {esc}^2 / (\delta v_\textrm {b})^2$, which leads via (2.6) and (2.13) to the approximate value of the coronal temperature,

(2.19)\begin{equation} T \simeq \frac{m_\textrm{p} v_\textrm{esc}^2}{8 k_\textrm{B} \ln(v_\textrm{esc}/\delta v_\textrm{b})}. \end{equation}

A number of factors can cause $\dot {M}$, $U_{\infty }$, and the coronal temperature to deviate from the estimates in (2.16) (2.18), and (2.19). For example, some of the AW power survives out to $r_\textrm {A}$, which reduces the total turbulent heating in the sub-Alfvénic region appearing on the right-hand side of (2.15), which, in turn, reduces $\dot {M}$. The heat that is conducted from the corona into the transition region adds a negative sink term to the right-hand sides of (2.14) and (2.15), which likewise acts to reduce $\dot {M}$. On the other hand, the AW pressure force helps drive plasma away from the Sun at the critical point, which acts to increase $\dot {M}$. These effects, as well as super-radial expansion of the magnetic field, are included in the more detailed solar-wind model developed in the next section.

3. Steady-state model of the transition region, corona, and solar wind

The lowest layer of the solar atmosphere is the chromosphere, which extends 2000–3000 km above the photosphere with a temperature $T$ ranging from several thousand $K$ to around $10^4$ K. Bounding the chromosphere from above is the transition region, a narrow layer approximately $100$ km thick, in which the density $\rho$ drops (and $T$ increases) by approximately two orders of magnitude. Above the transition region lies the corona, which extends out a few solar radii from the Sun and has a temperature of around $10^6$ K. The corona contains both closed magnetic loops, which connect back to the Sun at both ends, and open magnetic-field lines that connect the solar surface to the distant interplanetary medium. Regions of the corona permeated by open magnetic-field lines have lower densities than closed-field-line regions and are referred to as coronal holes. In the analysis to follow, $r_\textrm {b}$ denotes the radius of the coronal base just above the transition region in Sun-centred spherical coordinates $(r, \theta , \phi )$, which is taken to be

(3.1)\begin{equation} r_\textrm{b} = 1.005 R_{{\odot}}. \end{equation}

The steady-state model developed in this section describes the outflowing plasma within an open magnetic flux tube from the transition region to beyond the Alfvén critical point $r_\textrm {A}$, at which the plasma outflow velocity $U$ equals $v_\textrm {A}$. Figure 1 provides a schematic overview of the model, which determines five unknowns – the density at the coronal base $\rho _\textrm {b}$, the temperature of the quasi-isothermal sub-Alfvénic region, the mass outflow rate $\dot {M}$, the asymptotic outflow velocity $U_\infty$, and the flux of heat from the coronal base into the transition region $q_\textrm {b}$ – through the following five steps:

Figure 1. Schematic overview of model.

1. (i) balancing turbulent heating at $r=r_\textrm {b}$ against radiative cooling at $r=r_\textrm {b}$;

2. (ii) balancing the total turbulent heating between $r_\textrm {b}$ and $r_\textrm {A}$ against the two primary sinks of internal energy in this region, $p\,\textrm {d} V$ work and the flux of heat into the transition region;

3. (iii) balancing, within the transition region, conductive heating against internal-energy losses from $p\,\textrm {d} V$ work, advection, and radiation;

4. (iv) equating the mass outflow rate at $r=r_\textrm {b}$ with the mass outflow rate at the wave-modified sonic critical point $r=r_\textrm {c}$; and

5. (v) equating the wave-modified Bernoulli integral at $r=r_\textrm {b}$ and $r=r_\textrm {c}$.

These five steps are detailed in §§ 3.4–3.7. Section 3.1 reviews some identities that follow from the conservation of mass and magnetic flux. Sections 3.2 and 3.3 review previous results on reflection-driven AW turbulence in the solar wind and present the simplified model of reflection-driven AW turbulence that is used in this paper. A glossary is given in table 1.

Table 1. Glossary.

3.1. Flux and mass conservation

To simplify the analysis, solar rotation is neglected, and the magnetic field is taken to be radial, except in the corona, where open magnetic-field lines fan out to fill the space above closed magnetic loops at lower altitudes. Mathematically,

(3.2)\begin{equation} B(r) = \frac{\bar{B} \eta(r) R_{{\odot}}^2}{r^2}, \end{equation}

where $\eta (r)$ is the local super-radial expansion factor, which approaches 1 when $r/R_{\odot } \gg 1$, and $\bar {B}$ is the magnetic-field strength that would arise at the photosphere in the absence of super-radial expansion (i.e. if $\eta (r)$ were unity everywhere.) Given (3.2), the magnetic-field strength at the coronal base is

(3.3)\begin{equation} B_\textrm{b} \equiv B(r_\textrm{b}) = \bar{B} \eta_\textrm{b} \psi, \end{equation}

where here and in the following a ‘b’ subscript indicates that the subscripted quantity is evaluated at $r=r_\textrm {b}$, and

(3.4)\begin{equation} \psi \equiv \frac{R_{{\odot}}^2}{r_\textrm{b}^2} = 0.9901. \end{equation}

Because rotation is neglected, the steady-state solar-wind outflow velocity is aligned with the background magnetic field (Mestel 1961):

(3.5)\begin{equation} \boldsymbol{v} = U \frac{\boldsymbol{B}}{B}. \end{equation}

The density and velocity satisfy the steady-state continuity equation,

(3.6)\begin{equation} {\boldsymbol{\nabla}\,\,\boldsymbol{\cdot}}\,\,(\rho \boldsymbol{v}) = 0. \end{equation}

It follows from (3.5), (3.6), and ${\boldsymbol{\nabla}\,\boldsymbol{\cdot}}\, \boldsymbol {B}=0$ that $\boldsymbol {B} \cdot \boldsymbol {\nabla }( \rho U / B) = 0$, and, hence,

(3.7)\begin{equation} \frac{\rho U}{B} = \mbox{constant}. \end{equation}

Equation (3.7) is equivalent to

(3.8)\begin{equation} \dot{M} \equiv A(r) \rho(r) U(r) = \mbox{constant}, \end{equation}

where

(3.9)\begin{equation} A(r) = \frac{4 {\rm \pi}r^2}{\eta(r)} \end{equation}

is the cross-sectional area of the flow, which satisfies $A(r) = \mbox {constant}/ B(r)$, as required by flux conservation. Equations (3.8) and (3.9) follow the convention of expressing the mass outflow rate as the total solar mass-loss rate that would arise if the local mass flux at some $r$ characterised the entire outflow at that $r$, even though the actual solar wind is comprised of different wind streams with different properties. All of the results of this paper can be applied to an individual flux tube accounting for some fraction $\nu$ of the total flow area by multiplying the right-hand side of (3.9) by $\nu$.

It follows from (2.1) and (3.7) that $\rho ^{1/2} U/v_\textrm {A}$ is a constant. Because $U(r_\textrm {A}) = v_\textrm {A}(r_\textrm {A})$, this constant must be $\rho _\textrm {A}^{1/2}$, where

(3.10)\begin{equation} \rho_\textrm{A} \equiv \rho(r_\textrm{A}). \end{equation}

Thus,

(3.11)\begin{equation} v_\textrm{A} = y U, \end{equation}

where

(3.12)\begin{equation} y \equiv \left(\frac{\rho}{\rho_\textrm{A}}\right)^{1/2}. \end{equation}

With the aid of (3.3), (3.4), (3.9) and (3.11), (3.8) can be rewritten as

(3.13)\begin{equation} \dot{M} = A_\textrm{b} \rho_\textrm{b} U_\textrm{b} = \frac{4 {\rm \pi}R_{{\odot}}^2}{\psi \eta_\textrm{b}} \times \rho_\textrm{b} \times \frac{v_\textrm{Ab}}{y_\textrm{b}} = \frac{R_{{\odot}}^2 \bar{B} (4 {\rm \pi}\rho_\textrm{b})^{1/2}}{y_\textrm{b}} = \bar{B} R_{{\odot}}^2 (4 {\rm \pi}\rho_\textrm{A})^{1/2}. \end{equation}

Thus, $\dot {M}$ is determined uniquely by the value of $\rho _\textrm {A}$ and the single-hemisphere open magnetic flux, $2 {\rm \pi}R_{\odot }^2 \bar {B}$.

3.2. Reflection-driven AW turbulence

Dmitruk et al. (2002) (hereafter D02) developed an analytic model of reflection-driven AW turbulence in the solar corona valid in the limit of small $L_\perp$, where $L_\perp$ is the correlation length of the AW fluctuations measured in the plane perpendicular to the background magnetic field. Chandran & Hollweg (2009) (hereafter CH09) generalised this model by accounting for the solar-wind outflow velocity. Section 3.2.1 summarises the main results of the CH09 model, and § 3.2.2 uses the CH09 model to show that heating by reflection-driven AW turbulence causes the sub-Alfvénic region of the solar wind (at $r<r_\textrm {A}$, where $U< v_\textrm {A}$) to become approximately isothermal. Section 3.2.3 presents a modified version of the CH09 model that is easier to work with analytically, which is used to incorporate AW turbulence into the solar-wind model developed in this section.

3.2.1. The Chandran & Hollweg (2009) model of reflection-driven AW turbulence

In classical mechanics, a simple harmonic oscillator with frequency $\omega$ and energy $E$ possesses an adiabatic invariant $E/\omega$. If the parameters of the oscillator vary on a time scale $t_0$ satisfying $t_0 \gg \omega ^{-1}$ (e.g. if the length of a pendulum is slowly varied), then $E/\omega$ is almost exactly conserved. As $(\omega t_0)^{-1} \rightarrow 0$, changes in $E/\omega$ vanish faster than any power of $(\omega t_0)^{-1}$ (Landau & Lifshitz 1960).

An AW is like a space-filling collection of harmonic oscillators, and the wave action is analogous to the harmonic oscillator's adiabatic invariant. The wave action per unit volume per unit $\omega$ is ${\mathcal {E}}_\omega / \omega ^\prime$, where $\omega$ is the AW frequency in an inertial frame centred on the Sun, $\omega ^\prime$ is the AW frequency in the local plasma frame, and ${\mathcal {E}}_\omega$ is the AW energy per unit volume per unit $\omega$. In the Wentzel-Kramers-Brillouin (WKB) limit, in which the wave period is much shorter than the time scale on which the plasma parameters vary appreciably and the wave length is much shorter than the length scales over which the background plasma varies appreciably, the wave action satisfies the conservation law

(3.14)\begin{equation} \frac{\partial }{\partial t} \left(\frac{{\mathcal{E}}_\omega}{\omega^\prime}\right) +{\boldsymbol{\nabla}\,\boldsymbol{\cdot}}\, \left(\frac{\boldsymbol{c} {\mathcal{E}}_\omega}{\omega^\prime}\right) = 0, \end{equation}

where $\boldsymbol {c}$ is the group velocity of the waves (Bretherton & Garrett 1968; Dewar 1970). For outward-propagating AWs in the solar wind, $\omega = k_r(U+v_\textrm {A})$, $\omega ^\prime = k_r v_\textrm {A}$, and $\boldsymbol {c} = (U+v_\textrm {A})\boldsymbol {\hat {r}}$, where $\boldsymbol {\hat {r}}$ is the radial unit vector and $k_r$ is the radial component of the wave vector. In a steady-state solar wind, $\omega$ depends on neither position nor time. Upon multiplying (3.14) by $\omega$, integrating over $\omega$, and assuming a steady state, one obtains

(3.15)\begin{equation} {\boldsymbol{\nabla}\,\boldsymbol{\cdot}}\, \left[ \frac{\boldsymbol{\hat{r}} (U+v_\textrm{A})^2 {\mathcal{E}}_\textrm{tot}}{v_\textrm{A}}\right] = 0, \end{equation}

where

(3.16)\begin{equation} {\mathcal{E}}_\textrm{tot} = \int {\mathcal{E}}_\omega \,\textrm{d} \omega = \tfrac{1}{4} \rho z_+^2 \end{equation}

is the total AW energy density,

(3.17)\begin{equation} z_\pm\ {\equiv}\ \langle |\boldsymbol{z}^\pm|^2 \rangle^{1/2}, \end{equation}

$\langle \dots \rangle$ indicates a time average,

(3.18)\begin{equation} \boldsymbol{z}^\pm\ {=}\ \delta \boldsymbol{v} \mp \frac{\delta \boldsymbol{B}}{\sqrt{4{\rm \pi} \rho}} \end{equation}

are the Elsasser variables, $\delta \boldsymbol {v}$ and $\delta \boldsymbol {B}$ are the fluctuating velocity and magnetic field, and $\boldsymbol {z}^+$ ($\boldsymbol {z}^-$) corresponds to AW fluctuations propagating away from (toward) the Sun.¹ In (3.17) and the following, a $\pm$ sign is used as a subscript (as opposed to a superscript) when the subscripted quantity is an r.m.s. value. As ${\boldsymbol{\nabla}\,\boldsymbol{\cdot}}\, (\rho U \boldsymbol {\hat {r}}) = 0$, (3.15) can be rewritten as

(3.19)\begin{equation} \frac{\mathrm{d}}{\mathrm{d}r}g^2 = 0, \end{equation}

where

(3.20)\begin{equation} g^2 = \frac{(U+v_\textrm{A})^2 z_+^2}{Uv_\textrm{A}} = \frac{(1+y)^2 z_+^2}{y} \end{equation}

is the wave-action flux per unit mass flux per unit $\omega$ times $4 \omega$ integrated over $\omega$, or, for brevity, the ‘wave action flux per unit mass flux’.

When the finite radial wavelength of the AWs is taken into account, the radial gradient in $v_\textrm {A}$ causes partial non-WKB reflection of $\boldsymbol {z}^+$ fluctuations, leading to the production of $\boldsymbol {z}^-$ fluctuations. Counter-propagating AWs then interact, causing fluctuation energy to cascade to small scales and dissipate. In the CH09 model, this loss of fluctuation energy causes $g^2$ to decrease with radius according to the equation²

(3.21)\begin{equation} (U+v_\textrm{A}) \frac{\mathrm{d}}{\mathrm{d}r}g^2 ={-} \frac{z_-}{L_\perp} g^2. \end{equation}

Equation (3.21) states that in a reference frame that follows an outward-propagating AW, the wave action flux per unit mass flux decays on the eddy turnover time scale $L_\perp / z_-$. The reason that only $z_-$ (and not $z_+$) appears in this eddy turnover time scale is that ${\boldsymbol z}^+$ fluctuations are not sheared or distorted by other $\boldsymbol {z}^+$ fluctuations, but they are sheared and distorted by counter-propagating $\boldsymbol {z}^-$ fluctuations (Iroshnikov 1963; Kraichnan 1965). The radial decay of $g^2$ is accompanied by turbulent heating at the rate

(3.22)\begin{equation} Q = \frac{\rho z_+^2 z_-}{4 L_\perp}, \end{equation}

which is the energy density of the outward-propagating AWs $\rho z_+^2 / 4$ divided by their eddy turnover time scale $L_\perp / z_-$. Equation (3.22) drops a term $\rho z_-^2 z_+ / (4 L_\perp )$ that is normally included in the turbulent heating rate because of CH09's assumption that

(3.23)\begin{equation} z_-\ {\ll}\ z_+ , \end{equation}

an inequality that holds in the small-$L_\perp$ limit, as can be seen later in (3.25).

Following D02, CH09 determined $z_-$ by balancing the rate at which $z_-$ is produced through reflections against the rate at which $z_-$ cascades to small scales through nonlinear interactions in the small-$L_\perp$ limit, obtaining

(3.24)\begin{equation} \frac{(U+v_\textrm{A})}{v_\textrm{A}}\left| \frac{\mathrm{d}v_\textrm{A}}{\mathrm{d}r}\right| z_+\ {=}\ \frac{z_+}{L_\perp} z_-, \end{equation}

or, equivalently,

(3.25)\begin{equation} \frac{z_-}{L_\perp} = \frac{(U+v_\textrm{A})}{v_\textrm{A}}\left| \frac{\mathrm{d}v_\textrm{A}}{\mathrm{d}r}\right|. \end{equation}

Upon substituting (3.25) into (3.21), assuming that $v_\textrm {A}(r)$ has a single maximum at $r_\textrm {m} > r_\textrm {b}$, and solving for $g(r)$, CH09 found that

(3.26)\begin{equation} g^2 = g_\textrm{b}^2 h(r) , \end{equation}

where

(3.27)\begin{equation} h(r) = \left\{ \begin{array}{@{}ll} v_\textrm{Ab}/v_\textrm{A}(r) & \mbox{ for } r_\textrm{b} < r < r_\textrm{m} \\ v_\textrm{Ab} v_\textrm{A}(r) / v_\textrm{Am}^2 & \mbox{ if } r > r_\textrm{m} \end{array}\right., \end{equation}

where $v_\textrm {Am} = v_\textrm {A}(r_\textrm {m})$ is the maximum value of $v_\textrm {A}$. Conceptually, (3.26) and (3.27) state that an appreciable fraction of the local AW action flux per unit mass flux dissipates within each Alfvén-speed scale height, which causes $z_+^2(r)$ to drop below the WKB value that would be predicted from (3.20) with constant $g^2$.

Upon substituting (3.26) into (3.22) and making use (3.20), CH09 obtained

(3.28)\begin{equation} Q = \rho (\delta v_\textrm{b})^2 \left(\frac{U+v_\textrm{A}}{v_\textrm{A}}\right) \left|\frac{\mathrm{d}v_\textrm{A}}{\mathrm{d}r}\right| \left(\frac{y}{y_\textrm{b}}\right)\left(\frac{1 + y_\textrm{b}}{1+y}\right)^2 h(r), \end{equation}

where

(3.29)\begin{equation} (\delta v)^2 \equiv \langle | \delta \boldsymbol{ v}|^2 \rangle = \frac{z_+^2}{4}. \end{equation}

The second equality in (3.29) follows from (3.23). As

(3.30)\begin{equation} y_\textrm{b} \gg 1, \end{equation}

equation (3.28) can be rewritten to a good approximation in the following simplified form with the aid of (3.11):

(3.31)\begin{equation} Q = \rho (\delta v_\textrm{b})^2 y_\textrm{b}\left(\frac{U }{U+v_\textrm{A} }\right) \left| \frac{\mathrm{d}v_\textrm{A}}{\mathrm{d}r}\right| h(r). \end{equation}

3.2.2. Approximate isothermality between the transition region and Alfvén critical point

In steady state, the plasma internal-energy equation takes the form

(3.32)\begin{equation} {\boldsymbol{\nabla}\,\boldsymbol{\cdot}}\,\left[\boldsymbol{v} \left(\frac{p}{\gamma-1}\right)\right] ={-}p {\boldsymbol{\nabla}\,\boldsymbol{\cdot}}\, \boldsymbol{v} -{\boldsymbol{\nabla}\,\boldsymbol{\cdot}}\,\boldsymbol{q} + Q - R, \end{equation}

where $p$ is the pressure, $p/(\gamma -1)$ is the internal-energy density, $-p {\boldsymbol{\nabla}\,\boldsymbol{\cdot}}\, \boldsymbol {v}$ is the rate at which $p\,\textrm {d} V$ work is done on the plasma per unit volume, $Q$ is the rate of turbulent heating per unit volume, $\boldsymbol {q}$ is the heat flux, which is written in the form³

(3.33)\begin{equation} \boldsymbol{q} = q_r \frac{\boldsymbol{B}}{B}, \end{equation}

and $R$ is the rate of radiative cooling per unit volume. In the corona and solar wind, the density is sufficiently small that radiative cooling can be neglected, and, to a good approximation, (3.32) can be rewritten with the aid of (3.6) as

(3.34)\begin{equation} - c_\textrm{s}^2 U \frac{\textrm{d} \rho}{\textrm{d} r} + \frac{\rho U}{\gamma-1} \frac{\mathrm{d}}{\mathrm{d}r} c_\textrm{s}^2 ={-} \frac{1}{A} \frac{\mathrm{d}}{\mathrm{d}r}\left(A q_r\right) + Q . \end{equation}

A generic consequence of heating by reflection-driven AW turbulence is that, when other forms of heating (including conduction) are subdominant, the flow becomes quasi-isothermal at $r_\textrm {b} < r<r_\textrm {A}$, meaning that

(3.35)\begin{equation} \left| \frac{1}{c_\textrm{s}^2} \frac{\mathrm{d}c_\textrm{s}^2}{\mathrm{d}r}\right| \ll \left| \frac{1}{\rho} \frac{\mathrm{d}\rho}{\mathrm{d}r}\right|, \end{equation}

whereas (3.35) is not satisfied at $r>r_\textrm {A}$. This can be demonstrated by (1) neglecting conductive heating in (3.34), (2) assuming that (3.35) is satisfied so that the second term on the left-hand side of (3.34) can be neglected, and (3) solving (3.34) for $c_\textrm {s}^2$. If the resulting expression for $c_\textrm {s}(r)^2$ satisfies (3.35), then the neglect of the second term on the left-hand side of (3.34) is self-consistent, and this expression for $c_\textrm {s}(r)^2$ is a reasonable approximation for the full solution of (3.34). On the other hand, if the resulting value of $c_\textrm {s}^2(r)$ does not satisfy (3.35), then (3.34) does not possess a quasi-isothermal solution.

Carrying out this procedure and equating the first term on the left-hand side of (3.34) with the turbulent heating term on the right-hand side (which is given by (3.31)), one obtains

(3.36)\begin{equation} c_\textrm{s}^2 = y_\textrm{b} (\delta v_\textrm{b})^2 \left(\frac{ L_\rho}{L_{v_\textrm{A}}}\right) \left(\frac{v_\textrm{A}}{U+v_\textrm{A}}\right)h(r), \end{equation}

where

(3.37)\begin{equation} L_\rho = \rho\left|\frac{\mathrm{d}\rho}{\mathrm{d}r}\right|^{{-}1} \qquad L_{v_\textrm{A}} = v_\textrm{A}\left|\frac{\mathrm{d}v_\textrm{A}}{\mathrm{d}r}\right|^{{-}1} \end{equation}

are the density and Alfvén-speed scale heights. These scale heights are generally some constant of order unity times $r$, with $L_\rho /L_{v_{\textrm {A}}}$ increasing by a factor of a few between the low corona and $r_\textrm {A}$. On the other hand, $h(r)$ decreases by a factor of a few between the low corona and $r_\textrm {A}$, so that the product $(L_\rho/L_{v_{\rm A}})h(r)$ varies quite weakly with $r$.⁴ The $v_\textrm {A}/(U+v_\textrm {A})$ term in (3.36) likewise exhibits little variation in the sub-Alfvénic region, ranging from $\simeq 1$ at $r=r_\textrm {b}$ to $0.5$ at $r=r_\textrm {A}$. On the other hand, $\rho$ varies by a factor of $\sim 10^5$ between the coronal base and $r_\textrm {A}$ (see, e.g. Cranmer et al. 2007). Thus, the approximate solution for $c_\textrm {s}^2$ in (3.36) satisfies (3.35) at $r_\textrm {b} < r < r_\textrm {A}$, and AW heating indeed causes the sub-Alfvénic region to become quasi-isothermal. In contrast, at $r>r_\textrm {A}$, $U$ asymptotes toward a constant value, $h \propto v_\textrm {A} \propto \rho ^{1/2} U \propto \rho ^{1/2}$, $v_\textrm {A}/(U+v_\textrm {A}) \propto \rho ^{1/2}$, and the right-hand side of (3.36) becomes proportional to $\rho$, contradicting (3.35).

3.2.3. A modified version of the Chandran & Hollweg (2009) model

There are two difficulties with incorporating the CH09 model into an analytic solar-wind model that includes the region immediately above the transition region. First, (3.25) is consistent with (3.23) only if

(3.38)\begin{equation} \frac{v_\textrm{A} L_\perp}{z_+} \ll \frac{v_\textrm{A}}{|\textrm{d} v_\textrm{A}/ \textrm{d} r|}. \end{equation}

The left-hand side of (3.38) is the characteristic distance a $\boldsymbol {z}^-$ fluctuation at scale $L_\perp$ (measured perpendicular to the background magnetic field) propagates along the magnetic field before cascading and dissipating, and the right-hand side is the Alfvén-speed scale height. Just above the transition region, the magnetic-field strength has a scale height of order $10^{-2} R_{\odot }$ (Hackenberg, Marsch & Mann 2000; Cranmer et al. 2007), the $v_\textrm {A}$ scale height has a similar value, and (3.38) is not satisfied (see, e.g. equation (3.13-a) and table 3 of Chandran & Perez 2019). On the other hand, beyond the low corona, the $v_\textrm {A}$ scale height grows to values $\sim r$ that are $\gg L_\perp$. Thus, the CH09 model should be applied only beyond some minimum heliocentric distance, and some other approximation is needed to treat the region immediately above the transition region.

The second difficulty with incorporating the CH09 model into an analytic solar-wind model is that the function $h(r)$ in (3.27) depends on the number of local extrema in the Alfvén-speed profile and the locations of these extrema. In models that account for the rapid increase in $B(r)$ as $r$ drops below $\simeq 1.01 R_{\odot }$, $v_\textrm {A}$ typically has two local extrema within the corona: a local minimum just above the coronal base and a local maximum a few tenths of a solar radius above the coronal base. (See, e.g. figure 3 of Cranmer & van Ballegooijen (2005) or figure 9 of Cranmer et al. (2007).) In the context of the analytic solar-wind model developed in the following, each local extremum introduces two new unknowns: the location of the extremum, and either the density or flow speed at the extremum.

In this paper, the first of the two aforementioned difficulties is handled by replacing the CH09 result for $Q(r)$ at $r=r_\textrm {b}$ with the expression

(3.39)\begin{equation} Q_\textrm{b} =\frac{\rho_\textrm{b} (\delta v_\textrm{b})^2 v_\textrm{Ab} }{l_\textrm{b}}, \end{equation}

where $l_\textrm {b}$ is a free parameter, which is set equal to $0.3 R_{\odot }$ in the numerical examples presented in § 4. The second difficulty is resolved by replacing (3.25) with the expression

(3.40)\begin{equation} \frac{z_-}{L_\perp} ={-} \sigma (U + v_\textrm{A}) \frac{\mathrm{d}}{\mathrm{d}r} \ln(1+y) \end{equation}

at $r_\textrm {b} < r < r_\textrm {A}$, where $\sigma$ is a free parameter, whose value is approximately 0.1–0.5 in the numerical examples in § 4. Whereas $z_- / L_\perp = (U+v_\textrm {A})/L_{v_\textrm {A}}$ in (3.25), (3.40) takes $z_- / L_\perp$ to be $\simeq (U+v_\textrm {A})/L_{\rho }$ times a free parameter.⁵ Because $\rho$ (and, hence, $y$) is a monotonically decreasing function of $r$, the minus sign on the right-hand side of (3.40) ensures that $z_- > 0$. The right-hand side of (3.40) is taken to be proportional to $\ln (1+y)$ instead of $\ln y$ simply to make some of the expressions encountered later on easier to integrate analytically. Another motivation for using (3.40) instead of (3.25) is that the parameter-free CH09 result in (3.25) overestimates the rate at which AWs cascade and dissipate in the solar wind (Chandran & Hollweg 2009), and the introduction of a free parameter in (3.40) in principle makes it possible to correct for this.

Upon substituting (3.40) into (3.21), solving for $g^2(r)$, and rewriting $g^2$ in terms of $z_+^2$ using (3.20), one finds that

(3.41)\begin{equation} z_+^2 = z_{+\textrm{b}}^2 \frac{y}{y_\textrm{b}} \left(\frac{1+y_\textrm{b}}{1+y} \right)^{2-\sigma} . \end{equation}

Equations (3.22), (3.29), (3.40), and (3.41) can then be used to show that

(3.42)\begin{equation} H = \chi_\textrm{H} P_\textrm{AW}(r_\textrm{b}), \end{equation}

where

(3.43)\begin{equation} H = \int_{r_\textrm{b}}^{r_\textrm{A}} Q(r) A(r) \,\textrm{d} r \end{equation}

is the total turbulent heating power between $r_\textrm {b}$ and $r_\textrm {A}$,

(3.44)\begin{equation} P_\textrm{AW}(r) = \rho (\delta v)^2 (U + v_\textrm{A}) A \end{equation}

is the power (energy flux times area) of outward-propagating AWs at radius $r$, and

(3.45)\begin{equation} \chi_\textrm{H} = 1 - 2^{\sigma -1} \left(\frac{2-\sigma}{1-\sigma}\right) \frac{(1+y_\textrm{b})^{1-\sigma}}{y_\textrm{b}} + \frac{1}{y_\textrm{b}(1-\sigma)} \end{equation}

is the fraction of $P_\textrm {AW}(r_\textrm {b})$ that is dissipated between $r_\textrm {b}$ and $r_\textrm {A}$.

3.3. Relating the AW amplitudes at the coronal base and photosphere

In the absence of wave reflection and dissipation, $P_\textrm {AW}$ (defined in (3.44)) would have almost exactly the same value at $r_\textrm {b}$ and $R_{\odot }$.⁶ However, the steep Alfvén-speed gradient in the transition region leads to strong AW reflection, and a vigorous energy cascade in the chromosphere leads to substantial AW dissipation (van Ballegooijen et al. 2011; Chandran & Perez 2019). This reduces $P_\textrm {AW}(r_\textrm {b}) / P_\textrm {AW }(R_{\odot })$ to some value $f_\textrm {chr} < 1$, where $f_\textrm {chr}$ is an effective AW transmission coefficient for the chromosphere and transition region. Equivalently, because $U\ll v_\textrm {A}$ at $r\leqslant r_\textrm {b}$,

(3.46)\begin{equation} P_\textrm{AW}(r_\textrm{b}) = \rho_\textrm{b} \delta v_\textrm{b}^2 v_\textrm{Ab} A_\textrm{b} = f_\textrm{chr} \rho_{{\odot}} \delta v_{{\odot}}^2 v_{\textrm{A}{\odot}} A_{{\odot}}, \end{equation}

where a $\odot$ subscript indicates that the subscripted quantity is evaluated at the photosphere. For the numerical calculations in § 4, I set

(3.47)\begin{equation} \rho_\odot\ {=}\ 10^{17} m_\textrm{p} \mbox{ cm}^{{-}3}. \end{equation}

As $B A = \mbox {constant}$ and $v_\textrm {A} = B /(4{\rm \pi} \rho )^{1/2}$, (3.46) can be rewritten as

(3.48)\begin{equation} \delta v_\textrm{b}^2 = \delta v_{{\odot}\textrm{eff}}^2 \left(\frac{\rho_{{\odot}}}{\rho_\textrm{b}}\right)^{1/2} , \end{equation}

where

(3.49)\begin{equation} \delta v_{{\odot}\textrm{eff}} = f_\textrm{chr}^{1/2} \delta v_{{\odot}}. \end{equation}

The difference between $\delta v_{\odot }$ and $\delta v_{\odot \textrm {eff}}$ is that $\delta v_{\odot }$ gives rise to the fluctuating velocity $\delta v_\textrm {b}$ at $r=r_\textrm {b}$ when reflection and dissipation are accounted for, whereas $\delta v_{\odot \textrm {eff}}$ would give rise to the same value of $\delta v_\textrm {b}$ via WKB propagation (i.e. without reflection or dissipation).

The value of $\delta v_{\odot \textrm {eff}}$ can be constrained in different ways. For example, observations fix $\delta v_{\odot }$ at a value of $\simeq 1 \mbox { km} \mbox { s}^{-1}$ (Richardson & Schwarzschild 1950), and numerical simulations suggest that $f_\textrm {chr}$ is $\simeq 0.04 \text {--} 0.08$ (van Ballegooijen et al. 2011; Chandran & Perez 2019). Alternatively, if the solar wind is assumed to be powered primarily by an AW energy flux, then $\delta v_{\odot \textrm {eff}}$ can be inferred directly from measurements of the mass flux and energy flux far from the Sun. This latter method is used to determine $\delta v_{\odot \textrm {eff}}$ in some of the numerical solutions presented in § 4 and is described in more detail in § 4.2. It is worth noting that in contrast to $\delta v_{\odot }$, $f_\textrm {chr}$, and $\delta v_{\odot \textrm {eff}}$, which are plausibly independent of the properties of the coronal plasma and coronal magnetic field, $\delta v_\textrm {b}$ depends upon $\rho _\textrm {b}$, which varies between different flux tubes with different super-radial expansion factors $\eta _\textrm {b}$, as shown later in (3.55).

3.4. Balancing turbulent heating and radiative cooling at the coronal base

Within the corona and transition region, the plasma is optically thin, and the rate of radiative cooling is given by

(3.50)\begin{equation} R = \left(\frac{\rho}{m_\textrm{p}}\right)^2 \varLambda(T), \end{equation}

where $m_\textrm {p}$ is the proton mass, and $\varLambda (T)$ is the optically thin radiative loss function. Figure 2 shows three different approximations to $\varLambda (T)$. The dashed lines correspond to equation (A1) of Rosner, Tucker & Vaiana (1978). The solid line is a plot of

(3.51)\begin{equation} \varLambda(T) = c_\textrm{R} T^{{-}1/2}, \end{equation}

where

(3.52)\begin{equation} c_\textrm{R} = 1.549 \times 10^{{-}19} \mbox{ erg} \mbox{ cm}^3 \mbox{ s}^{{-}1} \mbox{ K}^{1/2}. \end{equation}

Equation (3.51) is the temperature derivative of equation (A3) of Rosner et al. (1978). The dotted line in figure 2 is a piecewise-continuous linear approximation to the low-temperature range of $\varLambda (T)$ in figure 1 of Cranmer et al. (2007), which is included to illustrate that optically thin radiative cooling becomes extremely weak at $T\lesssim 10^4$ K.

Figure 2. Three approximations to the optically thin radiative loss function $\varLambda (T)$; see the text for further discussion.

Throughout most of the corona, $\rho$ is sufficiently small that radiative cooling is negligible. However, as $r$ decreases within the corona, $\rho$ increases by orders of magnitude, which causes $R/Q$ to increase, because $R\propto \rho ^2$. There is thus some radius $r_\textrm {b}$ at which

(3.53)\begin{equation} R(r_\textrm{b}) = Q(r_\textrm{b}), \end{equation}

which marks the transition between the corona, in which radiative cooling is negligible, and the low solar atmosphere, in which radiation is thermodynamically important. In this paper, the radius $r_\textrm {b}$ is identified as the base of the corona, as already noted in (3.1). As $r$ decreases below $r_\textrm {b}$, $\rho (r)$ increases above $\rho _\textrm {b}$, $R/Q$ increases to values $\gg 1$, and the temperature gradient length scale must decrease so that conductive heating can balance radiative cooling (as well as internal-energy losses from advection and $p\,\textrm {d} V$ work). This shortening of the temperature gradient length scale gives rise to the transition region, which is discussed further in § 3.6. To facilitate an analytic solution, I neglect turbulent heating at $r<r_\textrm {b}$ and radiative cooling at $r>r_\textrm {b}$.

With the aid of (3.39), (3.46), and (3.50), (3.53) can be written in the form

(3.54)\begin{equation} \frac{\rho_\textrm{b}^2}{m_\textrm{p}^2} \varLambda(T_\textrm{b}) = \frac{\rho_\textrm{b} (\delta v_\textrm{b})^2 v_\textrm{Ab} }{l_\textrm{b}} = \frac{\rho_{{\odot}} (\delta v_{{\odot}\textrm{eff}})^2 v_{\textrm{A} \odot} B_\textrm{b} }{l_\textrm{b} B_{{\odot}}}. \end{equation}

Solving for $\rho _\textrm {b}$, one finds via (3.51) that

(3.55)\begin{equation} \rho_\textrm{b} = \left(\frac{\delta v_{{\odot}\textrm{eff}}^2 \psi \bar{B} \eta_\textrm{b} m_\textrm{p}^{5/2} \overline{c_\textrm{s}}}{c_\textrm{R} l_\textrm{b}}\right)^{1/2} \left(\frac{\rho_{{\odot}}}{8 {\rm \pi}k_\textrm{B}}\right)^{1/4}, \end{equation}

where

(3.56)\begin{equation} \overline{ c_\textrm{s}} = \mbox{constant} \end{equation}

is the sound speed (2.6) in the sub-Alfvénic region ($r_\textrm {b} < r < r_\textrm {A}$), which is approximated as a constant for the reasons discussed in § 3.2.2. Equation (3.55) can be rewritten as

(3.57)\begin{equation} \rho_\textrm{b} = \frac{B_\textrm{ref}^2}{4 {\rm \pi}v_\textrm{esc}^2} \tilde{\rho}_{{\odot}}^{1/4} B_\ast^2 \xi^2 \psi^{1/2} x^{1/4}, \end{equation}

where

(3.58)\begin{equation} x = \frac{\overline{ c_\textrm{s}}^2}{v_\textrm{esc}^2} \end{equation}

is the dimensionless temperature of the sub-Alfvénic region,

(3.59)\begin{equation} B_\textrm{ref} = \left(\frac{4 {\rm \pi}m_\textrm{p}^{5/2} v_\textrm{esc}^6}{c_\textrm{R} R_{{\odot}} \sqrt{2 k_\textrm{B }}}\right)^{1/2} = 118.8 \mbox{ G} \end{equation}

is the magnetic-field strength for which $v_\textrm {A} = v_\textrm {esc}$ when the radiative cooling time $\rho \overline { c_\textrm {s}}^2 / R$, free-fall time $R_{\odot }/v_\textrm {esc}$, and sound-crossing time $R_{\odot }/\overline {c_\textrm {s}}$ are all equal,

(3.60)\begin{gather} \epsilon_{{\odot}} = \frac{(\delta v_{{\odot}\textrm{eff}})^2}{v_\textrm{esc}^2}, \end{gather}

(3.61)\begin{gather} \tilde{\rho}_{{\odot}} = \frac{4 {\rm \pi}v_\textrm{esc}^2 \rho_{{\odot}} }{B_\textrm{ref}^2}, \end{gather}

(for reference, $\tilde {\rho }_\odot = 5.7 \times 10^5$ given (3.47)), and

(3.62)\begin{equation} \xi = \left(\frac{\epsilon_{{\odot}} \eta_\textrm{b}}{B_\ast^3 \tilde{l}_\textrm{b}}\right)^{1/4} \qquad B_\ast\ {=}\ \frac{\bar{B}}{B_\textrm{ref}} \qquad \tilde{l}_\textrm{b} = \frac{l_\textrm{b}}{R_{{\odot}}}. \end{equation}

Equation (3.57) can used to express $v_\textrm {Ab}$ in the form

(3.63)\begin{equation} v_\textrm{Ab} = v_\textrm{esc} \psi^{3/4} \xi^{{-}1} \eta_\textrm{b} (x\tilde{\rho}_\odot)^{{-}1/8}, \end{equation}

and (3.48), (3.57), and (3.60) can be used to write

(3.64)\begin{equation} \epsilon \equiv \frac{(\delta v_\textrm{b})^2}{v_\textrm{esc}^2} = \epsilon_{{\odot}} \tilde{\rho}_\odot^{3/8} B_\ast^{{-}1} \xi^{{-}1} \psi^{{-}1/4} x^{{-}1/8}. \end{equation}

3.5. Internal-energy equilibrium within the sub-Alfvénic region

In the quasi-isothermal approximation (3.35), the second term on the left-hand side of (3.34) is negligible, and the volume integral of (3.34) between $r=r_\textrm {b}$ and $r=r_\textrm {A}$ yields

(3.65)\begin{equation} -\!\dot{M} \overline{ c_\textrm{s}}^2 \int_{r_\textrm{b}}^{r_\textrm{A}} \frac{1}{\rho}\frac{\mathrm{d}\rho}{\mathrm{d}r} \textrm{d} r ={-} \left(A q_{r}\right) |^{r_\textrm{A}}_{r_\textrm{b}} + \chi_\textrm{H} \dot{M} \delta v_\textrm{b}^2 (1 + y_\textrm{b}), \end{equation}

where (3.42)–(3.44) have been used to rewrite the volume integral of the turbulent heating rate, (3.8) has been used to pull a constant factor of $\dot {M} = \rho A U$ out of the integral on the left-hand side, and (3.11) has been used to write $U_\textrm {b} + v_\textrm {Ab} = U_\textrm {b} (1 + y_\textrm {b})$.

The heat flux at $r=r_\textrm {A}$ is a small fraction of the total energy flux in AW-driven solar-wind models and is thus neglected in (3.65). Upon evaluating the integral on the left-hand side of (3.65), noting that $q_{r\textrm {b}} = - |q_{r \textrm {b}}| \equiv - q_\textrm {b}$, making use of (3.12), and rearranging terms, one obtains

(3.66)\begin{equation} \chi_\textrm{H} \dot{M}\delta v_\textrm{b}^2 (1 + y_\textrm{b}) = 2 \dot{M}\overline{ c_\textrm{s}}^2 \ln y_\textrm{b} + A_\textrm{b} q_\textrm{b}. \end{equation}

The left-hand side of (3.66) is the total turbulent heating rate within the sub-Alfvénic region, which represents the source of internal energy between $r_\textrm {b}$ and $r_\textrm {A}$. The two terms on the right-hand side of (3.66) are the dominant sinks of internal energy in the sub-Alfvénic region: $p\,\textrm {d} V$ work and thermal conduction into the transition region. Dividing (3.66) by $\dot {M} v_\textrm {esc}^2$ leads to

(3.67)\begin{equation} \epsilon \chi_\textrm{H} (1 + y_\textrm{b}) = 2 x \ln y_\textrm{b} + \frac{ q_\textrm{b}}{\rho_\textrm{b} U_\textrm{b} v_\textrm{esc}^2}. \end{equation}

3.6. The flux of heat from the corona into the transition region

The temperature structure within the transition region can be determined using a method similar to the methods of Rosner et al. (1978) and Schwadron & McComas (2003). The Knudsen number $N_\textrm {K}$ (the electron Coulomb mean free path $\lambda _\textrm {mfp}$ divided by the temperature gradient scale length $l_T$) is approximately $10^{-3}$ in the low corona, approximately $10^{-3}$ at the upper end of the transition region, and approximately $10^{-6}$ at the lower end of the transition region.⁷ Because $N_\textrm {K} \ll 1$, the radial component of the heat flux in the transition region is well approximated by the Spitzer & Härm (1953) formula,

(3.68)\begin{equation} q_r ={-}\alpha T^{5/2} \frac{\mathrm{d}T}{\mathrm{d}r}, \end{equation}

where

(3.69)\begin{equation} \alpha =\frac{1.84 \times 10^{{-}5} \mbox{ erg} \mbox{ cm}^{{-}1} \mbox{ s}^{{-}1} \mbox{ K}^{{-}7/2}}{\ln \varLambda_\textrm{Coul}}, \end{equation}

and $\ln \varLambda _\textrm {Coul}$ is the Coulomb logarithm. In the numerical examples of § 4,

(3.70)\begin{equation} \ln \varLambda_\textrm{Coul} = 18.1, \end{equation}

the value for electron-electron collisions in a proton-electron plasma with $\rho /m_\textrm {p} = 10^{9} \mbox { cm}^{-3}$ and $T = 10^6 \mbox { K}$ (Huba 2013). The magnetic field near $r=r_\textrm {b}$ is taken to be approximately radial, and the width of the transition region ($\sim 10^2 \mbox { km}$) is so narrow that $p$, $B$, and $A$ are treated as constants within the transition region, with

(3.71)\begin{equation} {\boldsymbol{\nabla}\,\boldsymbol{\cdot}}\, \boldsymbol{q} = \frac{\mathrm{d}q_r}{\mathrm{d}r}. \end{equation}

As mentioned previously, turbulent heating is neglected at $r< r_\textrm {b}$. Within the transition region, (3.32) thus becomes

(3.72)\begin{equation} \left(\frac{\gamma}{\gamma-1}\right) p {\boldsymbol{\nabla}\,\boldsymbol{\cdot}}\, \boldsymbol{v} ={-} \frac{\mathrm{d}q_\textrm{r}}{\mathrm{d}r} - \frac{c_\textrm{R} p^2}{4 k_\textrm{B}^2 T^{5/2}}, \end{equation}

where $\rho$ has been expressed in terms of $p$ and $T$ using (2.6). The velocity divergence in (3.72) can be expressed in terms of the heat flux via

(3.73)\begin{equation} {\boldsymbol{\nabla}\,\boldsymbol{\cdot}}\, \boldsymbol{v} ={-} \frac{\boldsymbol{v}}{\rho} \boldsymbol{\cdot} \boldsymbol{\nabla} \rho = \frac{\boldsymbol{v}}{T} \boldsymbol{\cdot} \boldsymbol{\nabla} T = \frac{U}{T} \left(-\frac{q_r}{\alpha T^{5/2}}\right). \end{equation}

where the first equality follows from (3.6). With the aid of (3.73), and using (3.68) to write

(3.74)\begin{equation} \frac{\mathrm{d}q_r}{\mathrm{d}r} ={-} \frac{\mathrm{d}q_r}{\mathrm{d}T} \frac{q_r}{\alpha T^{5/2}} , \end{equation}

one can rewrite (3.72) in the form

(3.75)\begin{equation} -\!a_2 q_r = q_r \frac{\mathrm{d}q_r}{\mathrm{d}T} - a_1, \end{equation}

where

(3.76)\begin{equation} a_1 = \frac{\alpha p^2 c_\textrm{R}}{4 k_\textrm{B}^2} \qquad a_2 = \left(\frac{\gamma}{\gamma-1}\right) \frac{p U}{T} \end{equation}

are both constants. Upon defining

(3.77)\begin{equation} w ={-} \frac{a_2}{a_1} q_\textrm{r}, \end{equation}

which is positive, one can rewrite (3.75) as

(3.78)\begin{equation} \frac{a_2^2}{a_1} \textrm{d} T = \frac{w }{1 + w} \textrm{d} w, \end{equation}

which can be integrated from the chromospheric values of $T$ and $w$, denoted $T_\textrm {chr}$ and $w_\textrm {chr}$, to the values of $T$ and $w$ at $r=r_\textrm {b}$, denoted $T_\textrm {b}$ and $w_\textrm {b}$. As $T_\textrm {chr} \ll T_\textrm {b}$ and the chromospheric heat flux is negligible ($q$ scaling as $T^{7/2}$ divided by the temperature-gradient scale length), the chromospheric terms are dropped, and the integral becomes

(3.79)\begin{equation} w_\textrm{b} - \ln(1+w_\textrm{b}) = a_3, \end{equation}

where

(3.80)\begin{gather} a_3 = \frac{a_2^2 T_\textrm{b}}{a_1} = \frac{2}{I_1^2} \left(\frac{\gamma}{\gamma-1}\right)^2 \left(\frac{U_\textrm{b}}{\overline{ c_\textrm{s}}}\right)^2, \end{gather}

(3.81)\begin{gather} I_1 = \left(\frac{\alpha c_\textrm{R} m_\textrm{p}}{4 k_\textrm{B}^3}\right)^{1/2} = 0.1581, \end{gather}

and the numerical value on the right-hand side of (3.81) is calculated using (3.70). Equation (3.79) can be solved in terms of the lower branch of the Lambert $W$ function, $W_{-1}$. This solution, in conjunction with (3.77), yields

(3.82)\begin{equation} q_\textrm{b} ={-}q_{r{\rm b}} ={-} \frac{a_1}{a_2}[1 + W_{{-}1}({-}e^{-(1+a_3)})], \end{equation}

which is positive since $\displaystyle W_{-1}(- e^{-(1+a_3)})<-1$.

If $U_\textrm {b}/ \overline { c_\textrm {s}}$ is sufficiently small, then $a_3 \ll 1$. In this low-Mach-number limit, (3.82) becomes, to leading order in $a_3$,

(3.83)\begin{equation} q_\textrm{b} = I_1 \rho_\textrm{b} \overline{ c_\textrm{s}}^3. \end{equation}

Equation (3.83) is equivalent to equation (3.15) of Rosner et al. (1978) when thermal conduction is the only source of heating in the transition region, i.e. when $f_\textrm {H}$ is set equal to zero in their equation (3.15). Equation (3.82) follows Schwadron & McComas (2003) in generalising the work of Rosner et al. (1978) to include internal-energy losses from $p\,\textrm {d} V$ work and advection. However, the analysis leading to (3.82) differs from that of Schwadron & McComas (2003) in that (3.82) completely neglects turbulent heating at $r<r_\textrm {b}$ and does not assume that $\int q(T) \,\textrm {d} T \propto T$.⁸

3.7. Constraints associated with the wave-modified sonic critical point

In the presence of mostly outward-propagating AWs, a radial background magnetic field, and a radial outflow, the momentum equation in the approximately isothermal sub-Alfvénic region takes the form

(3.84)\begin{equation} \rho U \frac{\mathrm{d}U}{\mathrm{d}r} ={-} \overline{ c_\textrm{s}}^2\frac{\mathrm{d}\rho}{\mathrm{d}r} - \frac{\mathrm{d}}{\mathrm{d}r}\left(\frac{\rho z_+^2}{8}\right) - \frac{G M_{{\odot}} \rho}{r^2}, \end{equation}

where $\rho z_+^2 / 8$ is the AW pressure, which is one-half the AW energy density (Dewar 1970). Equation (3.7) implies that

(3.85)\begin{equation} \frac{1}{\rho}\frac{\mathrm{d}\rho}{\mathrm{d}r} + \frac{1}{U} \frac{\mathrm{d}U}{\mathrm{d}r} - \frac{1}{B}\frac{\mathrm{d}B}{\mathrm{d}r} = 0. \end{equation}

Equations (3.41) and (3.85) can be used to rewrite (3.84) in the form

(3.86)\begin{equation} \frac{1}{\rho}\frac{\mathrm{d}\rho}{\mathrm{d}r} \left\{{-}U^2 + \overline{ c_\textrm{s}}^2 + \frac{(\delta v_\textrm{b})^2 (1+y_\textrm{b})^{2-\sigma}[y^2(1+\sigma) + 3y]}{4 y_\textrm{b}(1+y)^{3-\sigma}}\right\} = \frac{2}{r}\left(\gamma_B U^2 - \frac{v_\textrm{esc}^2 R_{{\odot}}}{4 r}\right), \end{equation}

where

(3.87)\begin{equation} \gamma_B \equiv{-} \frac{r}{2B} \frac{\mathrm{d}B}{\mathrm{d}r} = 1 - \frac{r}{2\eta} \frac{\mathrm{d}\eta}{\mathrm{d}r}, \end{equation}

and $\eta$ is defined in (3.2). In order for (3.86) to possess a transonic-wind solution for $U$, the quantity in braces on the left-hand side must be positive near the Sun and negative far from the Sun; i.e. it must pass through zero at some radius $r_\textrm {c}$ (the wave-modified sonic critical point). In order for $\textrm {d} \rho / \textrm {d} r$ to remain finite at $r=r_\textrm {c}$, the right-hand side of (3.86) must also vanish at $r_\textrm {c}$. Together, these conditions yield the constraints

(3.88)\begin{equation} r_\textrm{c} = \frac{v_\textrm{esc}^2 R_{{\odot}}}{4 \gamma_{B{\rm c}} U_\textrm{c}^2} \end{equation}

and

(3.89)\begin{equation} U_\textrm{c}^2 = \overline{ c_\textrm{s}}^2 + \frac{(\delta v_\textrm{b})^2 (1+y_\textrm{b})^{2-\sigma}[y_\textrm{c}^2(1+\sigma) + 3y_\textrm{c}]}{4 y_\textrm{b}(1+y_\textrm{c})^{3-\sigma}} , \end{equation}

where, here and in the following, a ‘c’ subscript indicates that the subscripted quantity is evaluated at $r=r_\textrm {c}$. An implicit assumption underlying (3.88) and (3.89) is that

(3.90)\begin{equation} r_\textrm{c} < r_\textrm{A}, \end{equation}

so that the quasi-isothermal approximation applies at $r_\textrm {c}$.

Mass and flux conservation (i.e. (3.7)) imply that

(3.91)\begin{equation} \frac{\rho_\textrm{c} U_\textrm{c}}{B_\textrm{c}} = \frac{\rho_\textrm{b} U_\textrm{b}}{B_\textrm{b}} = \frac{\rho_\textrm{b} v_\textrm{Ab}}{B_\textrm{b} y_\textrm{b}}, \end{equation}

where the second equality in (3.91) follows from (3.11). Equation (3.2) implies that

(3.92)\begin{equation} \frac{B_\textrm{b}}{B_\textrm{c}} = \frac{r_\textrm{c}^2}{R_{{\odot}}^2} \psi \frac{\eta_\textrm{b}}{\eta_\textrm{c}} = \frac{v_\textrm{esc}^4}{16 \gamma_{B{\rm c}}^2 U_\textrm{c}^4} \psi \frac{\eta_\textrm{b}}{\eta_\textrm{c}}, \end{equation}

where the second equality in (3.92) follows from (3.88). Upon substituting (3.92) into (3.91) and evaluating $v_\textrm {Ab}$ using (3.63), one obtains

(3.93)\begin{equation} U_\textrm{c}^2 = v_\textrm{esc}^2 \left[\frac{y_\textrm{c}^2 \psi^{1/4} \xi (x \tilde{\rho}_{{\odot}})^{1/8}}{16 \gamma_{B{\rm c}}^2 \eta_\textrm{c} y_\textrm{b}}\right]^{2/3}. \end{equation}

Substituting (3.93) into (3.89) yields

(3.94)\begin{equation} \left[\frac{y_\textrm{c}^2 \psi^{1/4} \xi (x \tilde{\rho}_{{\odot}})^{1/8}}{16 \gamma_{B{\rm c}}^2 \eta_\textrm{c} y_\textrm{b}}\right]^{2/3} - x - \frac{\epsilon (1+y_\textrm{b})^{2-\sigma}[y_\textrm{c}^2(1+\sigma) + 3 y_\textrm{c}]}{4 y_\textrm{b}(1+y_\textrm{c})^{3-\sigma}} = 0. \end{equation}

The integral over $r$ of $\rho ^{-1}$ times the momentum equation (3.84) yields the Bernoulli integral,

(3.95)\begin{align} \frac{U^2}{2} + \overline{ c_\textrm{s}}^2 \ln\left(\frac{\rho}{\rho_\textrm{b}}\right) - \frac{v_\textrm{esc}^2 R_{{\odot}}}{2 r} - \frac{(\delta v_\textrm{b})^2 (1+y_\textrm{b})^{2-\sigma}}{2y_\textrm{b}} \left[ \left(\frac{1+\sigma}{1-\sigma}\right) (1+y)^{\sigma - 1} + (1+y)^{\sigma -2}\right] = \varGamma, \end{align}

where $\varGamma$ is independent of $r$. Evaluating (3.95) at $r=r_\textrm {b}$ leads to the equality

(3.96)\begin{equation} \varGamma = \frac{U_\textrm{b}^2}{2} - \frac{v_\textrm{esc}^2 \psi^{1/2}}{2} - \frac{(\delta v_\textrm{b})^2}{2}\left[ \frac{1+\sigma}{1-\sigma} + \frac{2}{y_\textrm{b}(1-\sigma)}\right]. \end{equation}

Evaluating (3.95) at $r=r_\textrm {c}$, rewriting $r_\textrm {c}$ and $U_\textrm {c}$ using (3.88) and (3.89), rewriting $\varGamma$ using (3.96), and multiplying the resulting equation by $2/v_\textrm {esc}^2$ yields

(3.97)\begin{align} &\psi^{1/2} + x\left[4 \ln\left( \frac{y_\textrm{c}}{y_\textrm{b}}\right) + 1 - 4 \gamma_{B{\rm c}}\right] + \frac{\epsilon (1 + y_\textrm{b})^{2-\sigma} \varPhi}{4 (1-\sigma)y_\textrm{b} (1+y_\textrm{c})^{3-\sigma}}\nonumber\\ &\quad - \frac{\eta_\textrm{b}^2 \psi^{3/2}}{y_\textrm{b}^2 \xi^2 (x \tilde{\rho}_{{\odot}})^{1/4}} + \epsilon \left[\frac{1+\sigma}{1-\sigma} + \frac{2}{y_\textrm{b}(1-\sigma)}\right] = 0 , \end{align}

where

(3.98)\begin{equation} \varPhi \equiv y_\textrm{c}^2[(1-4\gamma_{B\textrm{c}})(1-\sigma^2) - 4(1+\sigma)] + y_\textrm{c}[3(1-4\gamma_{B\textrm{c}})(1-\sigma) - 12 - 4\sigma] - 8. \end{equation}

3.8. Mathematical structure of the model and approximate analytic solutions

The various quantities appearing in the model equations can be divided into five groups: (1) quantities that are determined observationally ($R_{\odot }$, $M_{\odot }$, $v_\textrm {esc}$, $\rho _{\odot }$, $\tilde {\rho }_{\odot }$, $\psi$, $\delta v_{\odot \textrm {eff}}$, $\bar {B}$); (2) free parameters ($\sigma$, $l_\textrm {b}$); (3) the super-radial expansion factor $\eta (r)$, which takes on different values in different magnetic flux tubes and in different models for the solar magnetic field; (4) the three principal unknowns, $y_\textrm {b}$, $y_\textrm {c}$, and $x$; and (5) additional unknowns that can be determined once $y_\textrm {b}$, $y_\textrm {c}$, and $x$ are found ($\dot {M}$, $\chi _\textrm {H}$, $U(r)$, $\rho _\textrm {b}$, $q_\textrm {b}$, $v_\textrm {Ab}$, $r_\textrm {c}$, $r_\textrm {A}$). The three principal unknowns $y_\textrm {b}$, $y_\textrm {c}$, and $x$ are determined by solving the three simultaneous equations (3.67), (3.94), and (3.97), where it must be remembered that $\epsilon$ is itself a function of $x$ via (3.64). The additional unknowns $\dot {M}$, $\chi _\textrm {H}$, $\rho _\textrm {b}$, $v_\textrm {Ab}$, $q_\textrm {b}$, and $r_\textrm {c}$ then follow immediately from (3.13), (3.45), (3.57), (3.63), (3.82), and (3.88), respectively. For example,

(3.99)\begin{equation} \dot{M} = \frac{R_\odot^2 \bar{B}^2}{v_\textrm{esc}} y_\textrm{b}^{{-}1} (x \tilde{\rho}_\odot)^{1/8} \xi \psi^{1/4} . \end{equation}

The procedures for determining $r_\textrm {A}$ and $U(r)$ involve a few more steps, which are described in Appendix A.

Two approximate analytic solutions to (3.67), (3.94) and (3.97), valid in two different parameter regimes, are derived in Appendix B. Both solutions rely on the approximations

(3.100)\begin{equation} y_\textrm{b} \gg 1 \qquad \psi = 1 \qquad \epsilon \ll 1 \qquad \eta_\textrm{c} = \gamma_{B \textrm{c}} = 1. \end{equation}

The last equality in (3.100) amounts to taking the magnetic-field lines to be purely radial at $r\geqslant r_\textrm {c}$. The first of the two approximate analytic solutions is valid in the conduction-dominated regime, in which the dominant sink of internal energy in the sub-Alfvénic region is the flux of heat from the corona into the transition region, rather than $p\,\textrm {d} V$ work. As discussed further in the following, this regime arises only for values of $\delta v_{\odot \textrm {eff}}$ much smaller than the solar value. This limit is thus not directly relevant to solar-wind observations. To leading order in $\epsilon _\odot$, the mass outflow rate $\dot {M}^\textrm {(cond)}$ and asymptotic flow velocity $U_{\infty }^\textrm {(cond)}$ in this parameter regime are given by

(3.101)\begin{equation} \dot{M}^\textrm{(cond)}=\frac{R_{{\odot}}^2 B_\textrm{ref}^2}{v_\textrm{esc}} \frac{1}{I_1^{1/14} I_2} [ \epsilon_\odot^{14-4\sigma} (\eta_\textrm{b} B_\ast)^{{-}4\sigma} \tilde{l}_\textrm{b}^{3\sigma}\tilde{\rho}_\odot^{7-2\sigma}]^{1/(7-7\sigma)}, \end{equation}

where $I_2$ is a numerical constant given in (B 12), and

(3.102)\begin{align} U_\infty^\textrm{(cond)} &= v_\textrm{esc}\Bigg[ \frac{2^{9-4\sigma} \zeta} {(1+\sigma)^2 I_1^{\sigma/14}} \left( \frac{1-\sigma}{7-3\sigma}\right)^{(7-3\sigma)/2}\nonumber\\ &\quad\vphantom{ \frac{2^{9-4\sigma} \zeta} {(1+\sigma)^2 I_1^{\sigma/14}}}\epsilon_\odot^{-(7-2\sigma)/7} B_\ast^{(7-5\sigma)/7} \eta_\textrm{b}^{2\sigma/7} \tilde{\rho}_\odot^{({-}7 + 2 \sigma)/14} \tilde{l}_\textrm{b}^{{-}3\sigma/14}\Bigg]^{1/2}. \end{align}

The coronal temperature in the conduction-dominated limit follows directly from (2.6) and (B 8).

The second approximate analytic solution is valid in the expansion-dominated regime, in which the $p\,\textrm {d} V$ work resulting from expansion is the dominant sink of internal energy in the sub-Alfvénic region. To leading order, the mass outflow rate in this parameter regime is given by

(3.103)\begin{equation} \dot{M}^\textrm{(exp)}_0 = \epsilon_\odot \bar{B} R_\odot^2 \sqrt{4{\rm \pi} \rho_\odot} = \frac{P_\textrm{AW}(r_\textrm{b})}{v_\textrm{esc}^2}, \end{equation}

where $P_\textrm {AW}(r_\textrm {b})$ is the AW power (3.44) evaluated at the coronal base. The asymptotic wind speed in the expansion-dominated regime is to leading order given by

(3.104)\begin{equation} U_{\infty, 0}^\textrm{(exp)} = v_\textrm{esc}, \end{equation}

and the coronal temperature in the expansion-dominated regime is to leading order

(3.105)\begin{equation} T = \frac{m_\textrm{p} v_\textrm{esc}^2}{k_\textrm{B} \ln(\epsilon_\odot^{{-}3} \tilde{\rho}_{{\odot}}^{{-}3/2} \tilde{l}_\textrm{b}^{{-}1} \eta_\textrm{b} B_\ast)}. \end{equation}

Equations (3.103) and (3.104) reproduce the approximate scalings of the simplified calculation presented in § 2. Equation (3.105) matches the right-hand side of (2.19) to within five percent for Sun-like parameters, the difference arising because in the model developed in this section, $\delta v_\textrm {b}$ has a weak dependence on the coronal temperature via (3.48) and (3.57) that is not accounted for in § 2. Appendix B presents higher-order corrections to (3.103), (3.104), and (3.105) that account for conductive losses, wave momentum deposition inside the wave-modified sonic critical point, and the fact that only part of $P_\textrm {AW}(r_\textrm {b})$ is dissipated within the sub-Alfvénic region. Second- and fourth-order approximations to $\dot {M}$ and $U_\infty$ are shown in figure 5 below.

Analytic estimates for the ranges of $\epsilon _{\odot }$ values corresponding to the conduction-dominated and expansion-dominated limits are given in Appendix B. There are three constraints on $\epsilon _\odot$ in the conduction-dominated limit, the most stringent of which is that the wave-energy term dominates over the internal-energy term in the Bernoulli integral at the wave-modified sonic critical point. The resulting range of $\epsilon _\odot$ values is much smaller than the solar value, as illustrated in figure 4. The expansion-dominated limit corresponds to a finite range of $\epsilon _\odot$ values that is sufficiently large that $p\,\textrm {d} V$ work dominates over conduction as the primary internal-energy sink within the sub-Alfvénic region, and sufficiently small that the sound speed makes the dominant contribution to the outflow velocity at the wave-modified sonic point in (3.89). This range of $\epsilon _\odot$ values is relevant to the solar case, as shown in the next section.

4. Numerical examples

This section presents several numerical solutions and approximate analytic solutions to the equations of the model developed in § 3. The numerical solutions are obtained by solving (3.67), (3.94), and (3.97) for $y_\textrm {b}$, $y_\textrm {c}$, and $x$ using Newton's method.⁹ Once $y_\textrm {b}$, $y_\textrm {c}$, and $x$ are determined, $\dot {M}$, $\chi _\textrm {H}$, $\rho _\textrm {b}$, $v_\textrm {Ab}$, $q_\textrm {b}$, $r_\textrm {c}$, $r_\textrm {A}$, and $U_{\infty }$ are computed from (3.13), (3.45), (3.57), (3.63), (3.82), (3.88), (A 3), and (A 13), respectively. The approximate analytic solutions are derived in Appendix B.

4.1. Magnetic-field model

The equations in § 3 are compatible with any model for the radial profile of the magnetic-field strength, or, equivalently, any choice of $\eta (r)$. On the other hand, the approximate analytic solutions derived in Appendix B assume that

(4.1)\begin{equation} \eta_\textrm{c} = \gamma_{B\textrm{c}} = 1. \end{equation}

To maintain consistency between the numerical and analytic solutions, and to avoid introducing additional complexity and free parameters, I take (4.1) to hold when computing the numerical solutions presented in this section. In essence, (4.1) amounts to assuming that all of the super-radial expansion of the magnetic field occurs inside the wave-modified sonic critical point. In measurements from the FIELDS experiment on the Parker Solar Probe (PSP) (Bale et al. 2016) during PSP's first few orbits, $|B_r| \simeq 2.2 \mbox { nT} (1 \mbox { a.u.}/r)^2$ (Sam Badman, private communication), where a.u. is the abbreviation for astronomical unit. In order to match this $B_r$ profile, I set

(4.2)\begin{equation} B_\ast\ {=}\ 0.00856 \hspace{0.3cm} \longleftrightarrow \hspace{0.3cm} B_r(\mbox{1 a.u.}) = 2.2 \mbox{ nT} \end{equation}

when computing the results shown in figures 3 through 6.

Figure 3. The outflow velocity $U$, density $\rho$, r.m.s. amplitude of the velocity fluctuation $\delta v$, and temperature $T$ as functions of heliocentric distance $r$ in a numerical solution to the model equations with $\eta _\textrm {b} = 100$, $\sigma = 0.5$, and the parameter values in (4.1), (4.2), (4.6), and (4.7). The dotted lines in the top-left panel are described in the text. In the top panels and the lower-right panel, the circles are measurements from a 1 h interval containing PSP's first perihelion on 6 November 2018, from figure 1 of Kasper et al. (2019). The error bars around the PSP data points in these three panels indicate the approximate range of values with a relative occurrence rate of at least 50 % of the peak occurrence rate within that 1 h interval. In the top-right panel, the error bars lie within the data point. In the lower-left panel, the PSP data point is from Chen et al. (2020), and the error bars around that data point show the approximate range of measured values near $r=35.7 R_{\odot }$ in figure 7 of Chen et al. (2020). The triangle in the lower-left panel is the value obtained by De Pontieu et al. (2007) from an analysis of the motion of filamentary structures in the low solar atmosphere based on observations from the Solar Optical Telescope on the Hinode satellite.

Figure 4. Mass outflow rate as a function of the effective fluctuating velocity at the photosphere $\delta v_{\odot \textrm {eff}} = f_\textrm {chr}^{1/2} \delta v_\odot$ (see (3.49)) for the parameter values in (4.1), (4.2), (4.7), $\eta _\textrm {b} = 30$, and $\sigma =0.18$, where $f_\textrm {chr}$ is the chromospheric/transition-region AW transmission coefficient in (3.46). The corresponding r.m.s. photospheric velocity $\delta v_\odot$ is shown at the top of the plot for the case in which $f_\textrm {chr} = 0.05.$ The solid line plots (3.99) using the numerical solution to (3.67), (3.94), and (3.97). The dotted and short-dashed lines plot the approximate analytic results from (3.101) and (3.103), respectively. The vertical dash-dot line corresponds to $\epsilon _{\odot \textrm {cond}}$ in (B 14). The left and right edges of the shaded region correspond to $\epsilon _{\odot \textrm {exp}, \textrm {min}}$ and $\epsilon _{\odot \textrm {exp}, \textrm {max}}$ in (B 35) and (B 36), respectively. To the right of the vertical long-dashed line, $r_\textrm {c} > r_\textrm {A}$, which violates (3.90) and the assumptions of the model.

Figure 5. The dependence of various flow properties on the Wang-Sheeley super-radial expansion factor $f_{\textrm {max}, \textrm {WS}}$ (4.8) for the parameter values in (4.1), (4.2), and (4.6)–(4.9). All quantities are evaluated using numerical solutions to the model equations, with the exception of $\dot {M}_2^\textrm {(exp)}$, $\dot {M}_4^\textrm {(exp)}$, $U_{\infty , 2}^\textrm {(exp)}$, and $U_{\infty , 4}^\textrm {(exp)}$, which are defined in (B 32) and (B 33), and the data points labeled ‘WS empirical’, which are taken from table 2 of Wang & Sheeley (1990). The horizontal error bars on these data points convey the half widths of the $f_{\textrm {max}, \textrm {WS}}$ data bins in that table. The vertical error bars correspond to one-half of the $100 \mbox { km} \mbox { s}^{-1}$ increment between the discretised $U_{\infty }$ values that define four of the five data bins. The quantity $I_1 \rho _\textrm {b} \overline {c_\textrm {s}}^3$ in the left panel of the second row is the low-Mach-number approximation to $q_\textrm {b}$ given in (3.83).

Figure 6. The dependence of various flow properties on the Wang-Sheeley super-radial expansion factor $f_{\textrm {max},\textrm {WS}}$ (4.8) and the effective fluctuating velocity at the photosphere $\delta v_{\odot \textrm {eff}}$ defined in (3.49).

4.2. Values of the AW power at the coronal base and $\delta v_{\odot {\rm eff}}$

Equations (3.46) and (3.49) and the flux-conservation relation $A_\textrm {b} B_\textrm {b} = A_{\odot } B_{\odot }$ imply that

(4.3)\begin{equation} \delta v_{{\odot} \textrm{eff}} = \left(\frac{4{\rm \pi}}{\rho_{{\odot}}}\right)^{1/4} \left( \frac{P_\textrm{AWb}}{B_\textrm{b} A_\textrm{b}} \right)^{1/2}. \end{equation}

When the solar wind is powered primarily by AWs, $P_\textrm {AWb}$ is approximately equal to the value of $\dot {E}_\textrm {m0}$ in (1.3), i.e.,

(4.4)\begin{equation} P_\textrm{AWb} \simeq \tfrac{1}{2} \dot{M} ( v_\textrm{esc}^2 + U_{\infty}^2). \end{equation}

Upon evaluating the right-hand side of (4.4) using the data in table 1 of Schwadron & McComas (2008) for Ulysses’ third northern polar pass (3NPP), one obtains

(4.5)\begin{equation} P_\textrm{AWb} = 3.6 \times 10^{27} \mbox{ erg} \mbox{ s}^{{-}1}, \end{equation}

which corresponds to $P_\textrm {AWb} / (4 {\rm \pi}R_{\odot }^2) = 0.59 \times 10^5 \mbox { erg} \mbox { cm}^{-2} \mbox { s}^{-1}$. I use the 3NPP data because this is the part of Ulysses's first three orbits during which the average scaled radial magnetic field $|\langle B_r \rangle | \cdot (r/ 1 \mbox { a.u.})^2$ was most consistent with (4.2) (in the 3NPP data, $|\langle B_r \rangle | \cdot (r/ 1 \mbox { a.u.})^2 = 2.1 \pm 0.08 \mbox { nT}$), and because $\dot {M}$ is strongly correlated with the average scaled radial magnetic field (Schwadron & McComas 2008). Equations (3.47), (4.2), (4.3), and (4.5) and the relation $B_\textrm {b} A_\textrm {b} = |B_r(1 \mbox { a.u.})| 4 {\rm \pi}(1 \mbox { a.u.})^2$ imply that

(4.6)\begin{equation} \delta v_{{\odot} \textrm{eff}} = 0.22 \mbox{ km} \mbox{ s}^{{-}1}. \end{equation}

This value of $\delta v_{\odot \textrm {eff}}$ is used in figures 3, 5, and 7.

Figure 7. The dependence of various flow properties on the Wang-Sheeley super-radial expansion factor $f_{\textrm {max},\textrm {WS}}$ (4.8) and the strength of the radial magnetic field at $r= 1 \mbox { a.u.}$

4.3. Free parameters

As discussed in § 3.8, there are two free parameters in the model: $l_\textrm {b}$ (the AW dissipation length scale at $r_\textrm {b}$) and $\sigma$ (the dimensionless coefficient in the turbulent heating rate). The solutions to (3.67), (3.94), and (3.97) are not very sensitive to the value of $l_\textrm {b}$. Throughout the rest of this section, $l_\textrm {b}$ is thus simply fixed at a value that seems physically reasonable:

(4.7)\begin{equation} l_\textrm{b} = 0.3 R_{{\odot}}. \end{equation}

On the other hand, the solutions depend sensitively on the value of $\sigma$. Larger values of $\sigma$ cause a larger fraction of the AW power at the coronal base $P_\textrm {AWb}$ to be dissipated in the quasi-isothermal sub-Alfvénic region, which increases $\dot {M}$ (see, e.g. (2.15) and (2.16)). For fixed $P_\textrm {AWb}$, increasing $\dot {M}$ decreases $U_{\infty }$. In some of the examples to follow, $\sigma$ is varied to optimise agreement between the model and observations, as described further in §§ 4.4 and 4.6. The super-radial expansion factor at the coronal base, $\eta _\textrm {b}$, also varies in the model, but this quantity is in principle observable for individual flux tubes. The value of $\eta _\textrm {b}$ is thus treated as an input into the model rather than an adjustable parameter.

4.4. Fiducial solution matching measurements from Parker Solar Probe's first perihelion

Figure 3 illustrates the $r$ dependence of the outflow velocity, density, fluctuating velocity, and temperature in a numerical solution to the model equations that is designed to agree with measurements from PSP's first perihelion encounter (E1) on 6 November 2018. The procedure for computing $U(r)$ and $T(r)$ is described in Appendix A. In order to solve for $U$ at some radius $r$, the value of $\eta$ must be specified at that $r$. For figure 3, I set $\eta (r) = 1 + (\eta _\textrm {b} - 1) \exp (-(r-r_\textrm {b})^2/(0.02 R_{\odot } )^2)$, which is effectively unity at $r= r_\textrm {c}$, consistent with (4.1). This choice of $\eta (r)$ is not realistic for the low corona (see, e.g. Cranmer et al. 2007), but the flow profile in the low corona is not a focus of this work. It should be noted that while the value of $\delta v$ at $r=r_\textrm {b}$ (i.e. $\delta v_\textrm {b}$) shown in the lower-left panel of figure 3 depends on $\eta _\textrm {b}$ through (3.48), (3.57), and (3.62), $\delta v_\textrm {b}$ does not on the way in which $\eta (r)$ decreases from $\eta _\textrm {b }$ to 1.

The solution in figure 3 is based upon the somewhat arbitrary assumption that $\eta _\textrm {b} = 100$ in the magnetic flux tube encountered by PSP at the time of its first perihelion. The value of $\sigma$ is set equal to 0.5 to optimise the agreement between the model and the data. Fits of comparable (in some cases superior) quality can be obtained for different values of $\eta _\textrm {b}$. For example, the parameter combinations $(\eta _\textrm {b}, \sigma ) = (300, 0.47)$ and $(\eta _\textrm {b}, \sigma ) = (30,0.55)$ lead to similar agreement with the data. Modelling of the solar magnetic field during PSP E1 suggests that the solar-wind stream encountered by PSP at the time of PSP's first perihelion originated in a small equatorial coronal hole (Bale et al. 2019). The fact that all four quantities plotted in figure 3 can be matched by varying the single parameter $\sigma$ suggests that the model is reasonably successful at capturing the physical processes that control the heating and acceleration of coronal-hole outflows.

The dotted lines in the top-left panel of figure 3 are solutions of the Bernoulli equation (3.95) for different values of the Bernoulli constant $\varGamma$. One of the dotted lines intersects the solid line, and that intersection occurs at the wave-modified sonic critical point, $r=r_\textrm {c}$. That dotted line is an accretion-like solution of (3.95) with the same value of $\varGamma$ as in the model solution but with $\textrm {d} U / \textrm {d} r < 0$ at $r=r_\textrm {c}$.

4.5. Illustration of the conduction-dominated and expansion-dominated regimes

Figure 4 plots $\dot {M}$ as a function of $\delta v_{\odot \textrm {eff}}$ when $\eta _\textrm {b} = 30$ and $\sigma = 0.18$. This figure includes the numerical solution to (3.67), (3.94), and (3.97) as well as the approximate analytic results $\dot {M}^\textrm {(cond)}$ and $\dot {M}^\textrm {(exp)}_0$ from (3.101) and (3.103). The vertical dash-dot line in this figure corresponds to $\epsilon _{\odot \textrm {cond}}$ in (B 14). The conduction-dominated approximation that gives rise to $\dot {M}^\textrm {(cond)}$ in (3.101) assumes that $\epsilon _\odot \ll \epsilon _{\odot \textrm {cond}}$.¹⁰ The shaded region corresponds to $\epsilon _{\odot \textrm {exp}, \textrm {min}}< \epsilon _\odot < \epsilon _{\odot \textrm {exp}, \textrm {max}}$, where $\epsilon _{\odot \textrm {exp}, \textrm {min}}$ and $\epsilon _{\odot \textrm {exp}, \textrm {max}}$ are defined in (B 35) and (B 36). The expansion-dominated approximation in § B.2 assumes that $\epsilon _{\odot \textrm {exp}, \textrm {min}}\ll \epsilon _\odot \ll \epsilon _{\odot \textrm {exp}, \textrm {max}}$. The vertical long-dashed line corresponds to the condition $r_\textrm {c} = r_\textrm {A}$. To the right of this line, $\dot {M}$ and the solar-wind density become so large that $r_\textrm {A} < r_\textrm {c}$, violating (3.90) and the assumptions underlying the model. In conjunction with (4.6), figure 4 illustrates that the expansion-dominated regime is relevant to the solar wind and that the conduction-dominated regime corresponds to values of $\delta v_{\odot \textrm {eff}}$ much smaller than the solar value.

4.6. Anti-correlation between the coronal super-radial expansion factor and $U_{\infty }$

Wang & Sheeley (1990) showed that the outflow velocity in a solar-wind stream at $r=1 {\rm a.u.}$ is anti-correlated with the super-radial expansion factor in the coronal magnetic flux tube from which the solar-wind stream originated. To compare their results with the model developed in § 3, I follow Cranmer et al. (2007) by defining the Wang & Sheeley (1990) coronal super-radial expansion factor $f_{\textrm {max}, \textrm {WS}}$ to be $\eta (1.04 R_{\odot })/ \eta (2.5 R_{\odot })$. I then compute $\eta (1.04 R_{\odot })$ and $\eta (2.5 R_{\odot })$ using a magnetic-field model similar to that used by Cranmer et al. (2007). In particular, I employ the global solar-magnetic-field model of Banaszkiewicz, Axford & McKenzie (1998) with the parameter values listed below their equation (2) ($K=1$, $M=1.789$, $a_1^\textrm {(B98)} = 1.538$, and $Q=1.5$). I supplement this global magnetic field with the low-solar-atmosphere magnetic-field model of Hackenberg et al. (2000) using the same parameter values listed in their figure 1 ($L = 30 \mbox { Mm}$, $d = 0.34 \mbox { Mm}$, $B_0 = 11.8 \mbox { G}$, and $B_\textrm {max} = 1.5 \mbox { kG}$). I then integrate out along a magnetic-field line at the edge of the polar corona hole in the Banaszkiewicz et al. (1998) model, compute $\eta (r)$, and obtain

(4.8)\begin{equation} f_{\textrm{max}, \textrm{WS}} = \frac{\eta_\textrm{b}}{5.75}. \end{equation}

I take (4.8) to hold even as $\eta _\textrm {b}$ is varied in figures 5, 6, and 7. This is equivalent to assuming that if local magnetic structures on the Sun cause $\eta _\textrm {b}$ to increase or decrease by some factor relative to the Hackenberg/Banaszkiewicz model just described, then $\eta (1.04 R_{\odot })$ increases or decreases by the same factor, but $\eta (2.5 R_{\odot })$ does not change.

If $\sigma$ is fixed while $\eta _\textrm {b}$ is varied, the model of § 3 does not agree with results of Wang & Sheeley (1990) shown in the top-right panel of figure 5. On the other hand, the model becomes consistent with those results if $\sigma$ is taken to have a power-law dependence on $\eta _\textrm {b}$ of the form

(4.9)\begin{equation} \sigma = 5.7\times 10^{{-}2} \, \eta_\textrm{b}^{0.34}. \end{equation}

Equation (4.9) is used to determine $\sigma$ in figures 5 through 7. The implication of (4.9) that $\sigma$ is an increasing function of $\eta _\textrm {b}$ is plausible because increasing $\eta _\textrm {b}$ increases $v_\textrm {Ab}$, as can be seen from (3.63). This, in turn, increases the number of Alfvén-speed scale heights between the coronal base and the Alfvén critical point ($\int _{r_\textrm {b}}^{r_\textrm {A}} (1/v_\textrm {A}) | \,\textrm {d} v_\textrm {A}/ \textrm {d} r|\, \textrm {d} r$), which increases the fraction of $P_\textrm {AWb}$ that is dissipated at $r<r_\textrm {A}$ (Chandran & Hollweg 2009). Further work, however, is needed to clarify how the structure of the coronal magnetic field influences the rate of turbulent dissipation along different magnetic flux tubes with different values of $\eta _\textrm {b}$ and to test the extent to which (3.41) and (4.9) are consistent with more rigorous treatments of AW turbulence in the sub-Alfvénic region of the solar wind.

The top-left panel of figure 5 shows that as $f_{\textrm {max}, \textrm {WS}}$ ranges from $\simeq 3$ to $\simeq 30$, $\dot {M}$ varies from $10^{-14} M_{\odot } \mbox { yr}^{-1}$ to $2\times 10^{-14} M_{\odot } \mbox { yr}^{-1}$, values that are similar to the solar-mass-loss rates inferred from Ulysses and PSP measurements (McComas et al. 2000; Schwadron & McComas 2008; Kasper et al. 2019). The fourth-order analytic approximation $\dot {M}_4^\textrm {(exp)}$ from (B 32) reproduces the numerical solution to the model equations reasonably well. The second-order analytic approximation $\dot {M}_2^\textrm {(exp)}$ is also reasonably accurate at $f_{\textrm {max}, \textrm {WS}} >rsim 3$, but deviates markedly from the numerical solution at $f_{\textrm {max}, \textrm {WS}} \lesssim 2$. A similar comment applies to the top-right panel of figure 5, which also shows that the model agrees fairly well with observational constraints on $U_{\infty }(f_{\textrm {max}, \textrm {WS}})$ when (4.9) holds. The left panel of the second row of figure 5 shows that the low-Mach-number approximation to $q_\textrm {b}$ given in (3.83) is quite accurate at small $f_{\textrm {max}, \textrm {WS}}$, where $U_\textrm {b}/ \overline {c_\textrm {s}} \ll 1$ (see the right panel in the second row of this figure). However, as $f_{\textrm {max}, \textrm {WS}}$ increases, $U_\textrm {b}$ increases, because the outflow is concentrated into a smaller cross-sectional area at the coronal base. This increase in $U_\textrm {b}$ leads to larger losses of internal energy within the transition region from $p\,\textrm {d} V$ work and advection, which, in turn, causes $q_\textrm {b}$ to increase above the low-Mach-number scaling so that conductive heating within the transition region can balance the additional non-radiative cooling. This same panel also shows that $q_\textrm {b}$ is significantly smaller than the AW energy flux at the coronal base ($\simeq \rho _\textrm {b} (\delta v_\textrm {b})^2 v_\textrm {Ab} = \rho _{\odot }^{1/2} (\delta v_{\odot \textrm {eff}})^2 B_\textrm {b} /\sqrt {4{\rm \pi} }$), which increases with $f_{\textrm {max}, \textrm {WS}}$ because increasing $f_{\textrm {max}, \textrm {WS}}$ increases $\eta _\textrm {b}$ and hence $B_\textrm {b}$. The lower-left panel of figure 5 shows that $\rho _\textrm {b}$ increases by a factor of $\simeq 10$ as $f_{\textrm {max}, \textrm {WS}}$ increases from 1 to 100, consistent with the $\rho _\textrm {b} \propto \eta _\textrm {b}^{1/2} x^{1/4}$ scaling in (3.57) and (3.62), given that $x = \overline {c_\textrm {s}}^2/v_\textrm {esc}^2$ is approximately constant over this range of $f_\textrm {max}$. The lower-right panel shows that $r_\textrm {c}$ is typically several solar radii, whereas $r_\textrm {A}$ ranges from $12 R_{\odot }$ to $15 R_{\odot }$ for this set of parameters.

4.7. Dependence of the flow properties on the AW power, super-radial expansion factor, and strength of the interplanetary magnetic field

The contour plots in figure 6 show how several quantities vary as functions of $f_{\textrm {max}, \textrm {WS}}$ and $\delta v_{\odot \textrm {eff}}$ in numerical solutions to the model equations based on the parameter values in (4.1), (4.2), (4.7), and (4.9). The top panels of figure 6 show that $\dot {M}$ is an increasing function of both $f_{\textrm {max}, \textrm {WS}}$ and $\delta v_{\odot \textrm {eff}}$, whereas $U_\infty$ is a decreasing function of both $f_{\textrm {max}, \textrm {WS}}$ and $\delta v_{\odot \textrm {eff}}$. The second row shows that $r_\textrm {A}$ is a strongly decreasing function of $\delta v_{\odot \textrm {eff}}$ and only weakly dependent on $f_{\textrm {max}, \textrm {WS}}$ for values of $\delta v_{\odot \textrm {eff}}$ comparable to the Sun-like value in (4.6). The right panel of this row shows that $r_\textrm {c}$ is an increasing function of $f_{\textrm {max}, \textrm {WS}}$ and, for the most part, a decreasing function of $\delta v_{\odot \textrm {eff}}$. The lower-left panel of figure 6 shows that $\overline {c_\textrm {s}}$ varies only weakly with $f_{\textrm {max}, \textrm {WS}}$ and $\delta v_{\odot \textrm {eff}}$. As shown in the lower-right panel, $\chi _\textrm {H}$ varies by a factor of $\simeq 2$ over the parameter range shown. If $\delta v_{\odot \textrm {eff}}$ is set equal to the value in (4.6) and $f_\textrm {max}$ is restricted to the interval $(2, 8)$ so that $U_{\infty }$ in the upper-right panel takes on fast-wind-like values of 600–800 $\mbox { km} \mbox { s}^{-1}$, then $\chi _\textrm {H} \simeq 0.5 \text {--}0.7$, consistent with direct numerical simulations of reflection-driven AW turbulence (Perez et al. 2021).

Figure 7 displays contour plots of the same quantities as in figure 6, but this time as functions of $f_{\textrm {max}, \textrm {WS}}$ and $B_\ast$, or, equivalently, $B_r(\mbox {1 a.u.})$, using the parameters in (4.1), (4.6), (4.7), and (4.9). The top panels show that $\dot {M}$ increases approximately linearly with $B_r(\mbox {1 a.u.})$ at fixed $f_{\textrm {max}, \textrm {WS}}$, and that $U_\infty$ varies only weakly with $B_r(\mbox {1 a.u.})$, consistent with measurements from the Ulysses spacecraft (Schwadron & McComas 2008; Riley et al. 2010). The left panel of the middle row shows that $r_\textrm {A}$ depends more strongly on $B_r(\mbox {1 a.u.})$ than on $f_{\textrm {max}, \textrm {WS}}$, whereas the reverse is true for $r_\textrm {c}$. As in figure 6, $\overline {c_\textrm {s}}$ varies very weakly across the entire parameter range observed. The lower-right panel shows that, when (4.9) holds, $\chi _\textrm {H}$ depends more strongly on $f_{\textrm {max}, \textrm {WS}}$ than on $B_r(\mbox {1 a.u.})$.

5. Discussion and conclusion

The main goal of this paper is to obtain an approximate analytic solution to the coupled problems of coronal heating and solar-wind acceleration under the assumption that the solar wind is powered primarily by an AW energy flux. Section 2 presents a first step toward this goal, namely, a simplified calculation of $\dot {M}$, $U_{\infty }$, and the coronal temperature in a spherically symmetric, steady-state solar wind in the absence of solar rotation. This calculation is based upon: (1) the assumption that $p\,\textrm {d} V$ work rather than thermal conduction is the primary sink of internal energy in the sub-Alfvénic region as a whole; (2) the result of § 3.2.2 that a solar or stellar wind heated primarily by AW turbulence becomes approximately isothermal between the coronal base $r_\textrm {b}$ and Alfvén critical point $r_\textrm {A}$; and (3) the finding in recent direct numerical simulations that most of the AW power at the coronal base $P_\textrm {AW}(r_\textrm {b})$ is dissipated between $r=r_\textrm {b}$ and $r=r_\textrm {A}$ (Perez et al. 2021). The calculation of § 2 shows that

(5.1)\begin{equation} \dot{M} \simeq \frac{P_\textrm{AW}(r_\textrm{b})}{v_\textrm{esc}^2} \qquad U_\infty \simeq v_\textrm{esc} \qquad T \simeq \frac{m_\textrm{p} v_\textrm{esc}^2}{8 k_\textrm{B} \ln(v_\textrm{esc}/\delta v_\textrm{b})} . \end{equation}

Deviations from (5.1) can be caused by a number of factors, including conductive losses into the transition region, wave momentum deposition inside the wave-modified sonic-critical point, and the fact that part of the AW power reaches the super-Alfvénic region, where it enhances $U_{\infty }$ without contributing to the heating that helps drive the outflow of mass past the wave-modified sonic critical point. These factors are accounted for in the more detailed solar-wind model developed in § 3. It is worth noting that in this model:

1. (i) the plasma density at the coronal base $\rho (r_\textrm {b}) = \rho _\textrm {b}$ is determined by equating the turbulent heating rate $Q(r_\textrm {b})$ and radiative cooling rate $R(r_\textrm {b})$; and

2. (ii) an analytic solution for the heat flux from the corona into the transition region $q_\textrm {b}$ is obtained by balancing, within the transition region, conductive heating against internal-energy losses from $p\,\textrm {d} V$ work, advection, and radiative cooling.

The expression (3.57) for $\rho _\textrm {b}$ that results from setting $Q(r_\textrm {b}) = R(r_\textrm {b})$ contains a factor of $x^{1/4}$, where $x = \overline { c_\textrm {s}}^2/v_\textrm {esc}^2$ is the dimensionless temperature. Thus, $x$ must be found before the exact value of $\rho _\textrm {b}$ can be determined. However, because $\dot {M}$ is so sensitive to the coronal temperature, $x^{1/4}$ can to a good approximation be treated as a constant (see, e.g. figures 6 and 7), and (3.57) with $x^{1/4} = 0.44$ can be used to obtain the approximate value of $\rho _\textrm {b}$ without solving the full model equations.

The equations of the solar-wind model developed in § 3 are solved analytically in two different parameter regimes. One of these is the conduction-dominated limit, in which heat conduction into the transition region is the dominant mechanism for draining internal energy from the sub-Alfvénic region, wave pressure makes the dominant contribution to the critical-point velocity $U_\textrm {c}$ in (3.89), and the wave-energy term dominates over the plasma-internal-energy term in the Bernoulli equation (3.95). This limit corresponds to photospheric velocities much smaller than those of the Sun. The second parameter regime is the expansion-dominated limit, in which $p\,\textrm {d} V$ work is the dominant sink of internal energy in the sub-Alfvénic region, the sound speed makes the dominant contribution to the critical-point velocity $U_\textrm {c}$ in (3.89), and the plasma-internal-energy term in the Bernoulli integral (3.95) dominates over the wave-energy term. As illustrated in figure 4, the expansion-dominated regime is relevant to the solar wind. The leading-order solution in the expansion-dominated regime reproduces (5.1), with a small difference in the coronal temperature arising from the fact that $\delta v_\textrm {b}$ has a weak dependence on the coronal temperature in the model of § 3. Numerical solutions to the model equations approach the approximate analytic solutions in the appropriate parameter regimes, match a range of solar-wind observations, and illustrate how the properties of the solar wind depend upon the r.m.s. photospheric velocity, super-radial expansion factor, and interplanetary magnetic-field strength.

5.1. Top-down causality for determining $\dot {M}$, $q_{\rm b}$, and the pressure, temperature range, and altitude of the transition region within the solar atmosphere

In this paper, as in Parker's original model (Parker 1958, 1965), the average rate at which mass flows out through the lower solar atmosphere is determined in large part by the outflow condition at the wave-modified sonic critical point $r_\textrm {c}$, several $R_\odot$ out from the Sun, at which the plasma transitions from being gravitationally bound to gravitationally unbound. Since the gravitational force weakens with increasing $r$, there is no physical mechanism that can prevent plasma at $r_\textrm {c}$ from flowing outward at approximately the wave-enhanced effective sound speed, $c_\textrm {s, eff} \equiv \displaystyle [(p + p_\textrm {wave})/\rho ]^{1/2} = [ c_\textrm {s}^2 + 0.5 (\delta v)^2]^{1/2}$, evaluated at $r=r_\textrm {c}$, which is comparable to the square root of the right-hand side of (3.89). This is why $\dot {M} \sim A(r_\textrm {c}) \rho _\textrm {c} c_{\textrm {s}, \textrm {eff}}(r_\textrm {c})$.

On the other hand, localised motions near the transition region at speeds $\ll v_\textrm {esc}$ are gravitationally bound, and the mass flux that they carry is not determinative of $\dot {M}$. For example, if at some time, such motions led to an overall mass outflow rate at $r_\textrm {b}$ exceeding the rate $\sim A(r_\textrm {c}) \rho _\textrm {c} c_{\textrm {s}, \textrm {eff}}(r_\textrm {c})$ at which mass flows past the critical point $r_\textrm {c}$, then plasma would build up in the corona. This would, in turn, weaken the pressure gradient relative to the gravitational force per unit volume in the vicinity of the transition region, thereby reducing the amount of plasma flowing up from the chromosphere.

Nevertheless, the transition region and chromosphere do affect $\dot {M}$ in two ways. First, $q_\textrm {b}$ reduces the net heating power within the quasi-isothermal sub-Alfvénic region, $P_\textrm {net}$. As discussed following (2.16), the heating cost per unit mass for plasma to transit the quasi-isothermal sub-Alfvénic region is $\overline { c_\textrm {s}}^2 \ln (\rho _\textrm {b}/\rho (r_\textrm {A})) \simeq v_\textrm {esc}^2$, and $\dot {M} \simeq P_\textrm {net}/v_\textrm {esc}^2$. By reducing $P_\textrm {net}$, the conduction of heat from the corona into the transition region reduces $\dot {M}$. The second way that the lower solar atmosphere affects $\dot {M}$ is via the chromospheric/transition-region transmission coefficient, $f_\textrm {chr} = P_\textrm {AW}(r_\textrm {b})/P_\textrm {AW}(R_{\odot })$, whose value again influences the value of $P_\textrm {net}$. To summarise this paragraph and the preceding paragraph, the chromosphere and transition region influence $\dot {M}$ thermodynamically, but not dynamically.

The regulation of the mass flux at $r_\textrm {b}$ by the critical-point condition at $r_\textrm {c}$ is an example of ‘top-down’ causality, in which physical processes at larger $r$ control the plasma properties at smaller $r$. In the model of this paper, top-down causality also characterises the determination of $q_\textrm {b}$, the pressure within the transition region, the altitude of the transition region in the solar atmosphere, and the temperature jump across the transition region. As mentioned previously, $\rho _\textrm {b}$ is determined by the condition that $Q(r_\textrm {b}) = R(r_\textrm {b})$, without reference to conditions in the chromosphere or the value of $q_\textrm {b}$. The sound speed at $r=r_\textrm {b}$, which is $\simeq \overline {c_\textrm {s}}$, is approximately determined by balancing turbulent heating against internal-energy losses from $p\,\textrm {d} V$ work within the corona and sub-Alfvénic solar wind. The outflow velocity at $r=r_\textrm {b}$, $U_\textrm {b}$, follows from the value of $\dot {M}$, which is controlled by the critical-point condition, as described above. The values of $\rho _\textrm {b}$, $\overline { c_\textrm {s}}$, and $U_\textrm {b}$ jointly determine $q_\textrm {b}$ via (3.82), which embodies the requirement that conductive heating offset internal-energy losses (from radiation, $p\,\textrm {d} V$ work, and advection) within the transition region. The values of $\rho _\textrm {b}$ and $\overline { c_\textrm {s}}$ are sufficient to determine the approximate transition-region pressure, $p_\textrm {tr} \simeq \rho _\textrm {b} \overline { c_\textrm {s}}^2$. The pressure within the upper chromosphere, $p_\textrm {chr}(r)$, is an approximately exponentially decreasing function of altitude. The altitude of the transition region is determined by setting $p_\textrm {chr}(r)= p_\textrm {tr}$. The shape of the radiative loss function plotted in figure 2 constrains the temperature at the bottom of the transition region to be approximately $10^4$ K, so that radiative cooling within the comparatively dense upper chromosphere can be balanced by local heating mechanisms in the absence of strong conductive heating. Given this constraint, the factor by which the temperature changes across the transition region is determined by the value of $T(r_\textrm {b})$, which is controlled by the balance between heating and $p\,\textrm {d} V$ work in the corona, as discussed previously.

Although the condition $R(r_\textrm {b}) = Q(r_\textrm {b})$ determines the density $\rho _\textrm {b}$ at the upper boundary of the transition region within the corona (subject to the caveats at the end of the paragraph following (5.1)), setting $R(r)=Q(r)$ within the chromosphere does not determine the density at the lower edge of the transition region, because small changes in $T$ within the upper chromosphere lead to dramatic changes in the radiative loss function $\varLambda (T)$, as shown in figure 2. The density at the bottom of the transition region is instead largely determined by the values of $\rho _\textrm {b}$ and $\overline {c_\textrm {s}}$, the near constancy of the pressure across the transition region, and the above-mentioned constraint (arising from the shape of the radiative loss function) that $T\sim 10^4$ K in the upper chromosphere.

5.2. Limitations and future work

The model of § 3 has a number of shortcomings. First, the sub-Alfvénic region is not truly isothermal, and, hence, the quasi-isothermal approximation in (2.5) leads to some error. Second, for the solutions shown in figure 5, the temperature $m_\textrm {p} \overline {c_\textrm {s}}^2/(2k_\textrm {B})$ of the sub-Alfvénic region increases from $\simeq 7\times 10^5 \mbox { K}$ when $U_\infty \simeq 800 \mbox { km} \mbox { s}^{-1}$ to $\simeq 8.4 \times 10^5 \mbox { K}$ when $U_\infty \simeq 400 \mbox { km} \mbox { s}^{-1}$. In contrast, in measurements from the Ulysses spacecraft, the coronal freeze-in temperature increases from $\simeq 9 \times 10^5 \mbox { K}$ when $U_{\infty } \simeq 800 \mbox { km} \mbox { s}^{-1}$ to $\simeq 1.35 \times 10^6 \mbox { K}$ when $U_\infty \simeq 400 \mbox { km} \mbox { s}^{-1}$ (McComas et al. 2002; Schwadron & McComas 2003). Although the numerical solutions in figure 5 reproduce the observed anti-correlation between coronal temperature and asymptotic wind speed, the isothermal-sub-Alfvénic-region approximation of the present paper is too crude to be able match the measured freeze-in temperatures in detail. Another shortcoming is that the dimensionless coefficient $\sigma$ that appears in the turbulent heating rate has a large impact on the solution to the model equations, but is an adjustable free parameter. Further work is needed to provide a physical basis for determining $\sigma$ and how it varies from one flux tube to another.

In a future study, the model developed in § 3 could be used in conjunction with studies that map the magnetic-field line traversed by PSP back to a source region on the Sun to rapidly predict flow properties at PSP's location based on the observed super-radial expansion factor within the source region (see, e.g. Bale et al. 2019; Kasper et al. 2019). In addition, the modelling results and approaches developed in this paper could be applied to outflows from other astrophysical objects, such as stars with differing masses and winds from accretion disks around compact astrophysical objects.