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

Between the base of the solar corona and the Alfven critical point, the solar-wind density decreases by approximately five orders of magnitude, but the temperature varies by a factor of only a few. In this paper, I show that such quasi-isothermal evolution out to the Alfven critical point is a generic property of outflows powered by reflection-driven Alfven-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-Alfvenic region as isothermal, I first present a simplified calculation of the mass outflow rate, asymptotic wind speed, and coronal temperature that neglects conductive losses and the wave pressure force. I then develop a more detailed model of the transition region, corona, and solar wind that accounts for the heat flux from the coronal base into the transition region and momentum deposition by AWs. I solve analytically for this heat flux by balancing, within the transition region, conductive heating against internal-energy losses from radiation, pdV work, and advection. The density at the coronal base is determined by locally balancing turbulent heating and radiative cooling. I solve the equations of the model analytically in two different parameter regimes. Analytic and numerical solutions to the model equations agree with a number of observations.


Introduction
Pioneering works by Parker (1958Parker ( , 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,Ṁ , because they specified the temperature of the coronal base as a boundary condition, andṀ 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 1973Hollweg , 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 et al. 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 et al. 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 et al. 2007;Verdini et al. 2010;Chandran et al. 2011;van der Holst et al. 2014;Usmanov et al. 2014). Threedimensional 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Ṁ and the outflow speed far from the Sun (U ∞ ) remain elusive. A number of studies have obtained a single equation that constrains the two unknownsṀ and U ∞ . For example, Sandbaek et al. (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 thaṫ , (1.1) whereĖ m0 is the mechanical-energy input into the solar wind at the corona base,Ė q0 is the energy input at the coronal base from thermal conduction (which is negative in models that include the lower solar atmosphere), is the escape velocity of the Sun, G is the gravitational constant, M ⊙ is the solar mass, and R ⊙ is the solar radius. Schwadron & McComas (2003) derived a variant of (1.1) that explicitly relatesĖ q0 to the altitude of the coronal-temperature maximum. Hansteen & Leer (1995) and Hansteen et al. (1997) showed thatĖ q0 ≪Ė m0 , and Hansteen & Velli (2012) made use of this finding to further refine (1.1), obtaininġ . (1.3) Another important result was obtained by Leer & Holzer (1980), who found that heating inside the sonic critical point enhancesṀ but has little effect on U ∞ , whereas heating beyond the sonic critical point increases U ∞ but has little effect onṀ .
The main goal of the present paper is to obtain approximate analytic solutions forṀ , U ∞ , 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Ṁ , U ∞ , 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Ṁ , U ∞ , and the coronal temperature in AW-driven winds with minimal conductive losses Approximate expressions forṀ , U ∞ , 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 B and flow velocity are aligned; (3) B ∝ r −2r , wherer 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 dV 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 where ρ 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 , r b is the radius of the coronal base, which in this section (but not the next) is simply set equal to R ⊙ , r A is the Alfvén critical point at which U = v A , v Ab = v A (r b ), and U b = U (r b ). The mass outflow rate can thus be written in the forṁ As shown in § 3.2, heating by reflection-driven AW turbulence causes the sub-Alfvénic part of the solar wind at r < r A to become quasi-isothermal, in the sense that is the square of the isothermal sound speed, k B is the Boltzmann constant, T is the temperature, and m 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 b and r A , whereas ρ 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 A can be written in the form (2.7) AsṀ = 4πr 2 U ρ is independent of r, U −1 dU/dr = −ρ −1 dρ/dr − 2/r. Upon substituting this relation into (2.7) and rearranging terms, one obtains (2.8) In order for (2.8) to have a smooth transonic solution, the right-hand side of (2.8) must vanish at the radius r c at which U = c s , so that dρ/dr remains finite. This leads to the two critical-point conditions where U c = U (r c ). Integrating (2.7), one obtains the Bernoulli integral where the second equality in (2.10) results from evaluating the left-hand side of (2.10) at r = r b and dropping the U 2 b /2 term, which is ≪ v 2 esc /2. After evaluating the left-hand side of (2.10) at r = r c and using (2.9) to rewrite U c and r c in terms of c 2 s , one obtains where ρ c = ρ(r c ). Upon settingṀ = 4πr 2 c ρ c U c and using (2.9) and (2.11) to rewrite r c , ρ c , and U c in terms of c s , one obtains (Hansteen & Velli 2012 (2.12) The exponential appearing on the right-hand side of (2.12) reflects the fact that, at r < r 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 where c 1 ≡ e 3/2 v 4 esc /(16c 3 s v Ab ). When radiative cooling and thermal conduction are neglected, the internal-energy equation becomes where Q is the turbulent heating rate, and γ 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πr 2 and integrating over the quasi-isothermal sub-Alfvénic region, one obtains where the approximate equality assumes that the volume-integrated turbulent heating rate between r b and r A is comparable to P AW (r b ), the power carried by outward-propagating AWs at the coronal base, consistent with direct numerical simulations (Perez et al. in press). As ln(ρ b /ρ(r A )) = 2 ln y b , equation (2.15) yieldṡ 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 dV work. The quantity c 2 s ln(ρ b /ρ(r A )) = 2c 2 s ln y b ≃ v 2 esc is thus the heating cost per unit mass for plasma to transit the sub-Alfvénic region. The reason that the product c 2 s ln(ρ b /ρ(r A )) is approximately constant is that increasing c 2 s leads to an exponential increase inṀ and the solar-wind density and an exponential reduction of ρ b /ρ(r A ), leaving c 2 s ln(ρ b /ρ(r A )) approximately unchanged. Equations (2.15) and (2.16) state thatṀ is the net heating power within the sub-Alfvénic region (which is taken to be ≃ P AW (r 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 b and at large r, one obtains (2.17) Substituting (2.16) into (2.17) yields Reflection-driven AW turbulence thus changes the energy per unit mass from ≃ −v 2 esc /2 at the coronal base to ≃ v 2 esc /2 far from the Sun. In the absence of super-radial expansion, where δv 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 b ≃ v 2 esc /(δv b ) 2 , which leads via (2.6) and (2.13) to the approximate value of the coronal temperature, A number of factors can causeṀ , U ∞ , 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 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Ṁ . 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Ṁ . On the other hand, the AW pressure force helps drive plasma away from the Sun at the critical point, which acts to increaseṀ . 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.

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 ρ 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 b denotes the radius of the coronal base just above the transition region in Sun-centred spherical coordinates (r, θ, φ), which is taken to be (3.1) 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 A , at which the plasma outflow velocity U equals v A . Figure 1 provides a schematic overview of the model, which determines five unknowns -the density at the coronal base ρ b , the temperature of the quasi-isothermal sub-Alfvénic region, the mass outflow rateṀ , the asymptotic outflow velocity U ∞ , and the flux of heat from the coronal base into the transition region q b -through the following five steps: 1. balancing turbulent heating at r = r b against radiative cooling at r = r b ; 2. balancing the total turbulent heating between r b and r A against the two primary sinks of internal energy in this region, p dV work and the flux of heat into the transition region; 3. balancing, within the transition region, conductive heating against internalenergy losses from p dV work, advection, and radiation; 4. equating the mass outflow rate at r = r b with the mass outflow rate at the wave-modified sonic critical point r = r c ; and 5. equating the wave-modified Bernoulli integral at r = r b and r = r 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.

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, where η(r) is the local super-radial expansion factor, which approaches 1 when r/R ⊙ ≫ 1, andB is the magnetic-field strength that would arise at the photosphere in the absence  of super-radial expansion (i.e., if η(r) were unity everywhere.) Given (3.2), the magneticfield strength at the coronal base is where here and in the following a 'b' subscript indicates that the subscripted quantity is evaluated at r = r b , and Because rotation is neglected, the steady-state solar-wind outflow velocity is aligned with the background magnetic field (Mestel 1961): (3.5) The density and velocity satisfy the steady-state continuity equation, ∇ · (ρv) = 0.
(3.6) It follows from (3.5), (3.6), and ∇ · B = 0 that B · ∇(ρU/B) = 0, and, hence, is the cross-sectional area of the flow, which satisfies A(r) = 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 ν of the total flow area by multiplying the right-hand side of (3.9) by ν. It follows from (2.1) and (3.7) that this constant must be ρ With the aid of (3.3), (3.4), (3.9) and (3.11), equation (3.8) can be rewritten aṡ Thus,Ṁ is determined uniquely by the value of ρ A and the single-hemisphere open magnetic flux, 2πR 2 ⊙B .

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 ⊥ , where L ⊥ 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 A , where U < v 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.

The Chandran & Hollweg (2009) model of reflection-driven AW turbulence
In classical mechanics, a simple harmonic oscillator with frequency ω and energy E possesses an adiabatic invariant E/ω. If the parameters of the oscillator vary on a time scale t 0 satisfying t 0 ≫ ω −1 (e.g., if the length of a pendulum is slowly varied), then E/ω is almost exactly conserved. As (ωt 0 ) −1 → 0, changes in E/ω vanish faster than any power of (ω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 ω is E ω /ω ′ , where ω is the AW frequency in an inertial frame centred on the Sun, ω ′ is the AW frequency in the local plasma frame, and E ω is the AW energy per unit volume per unit ω. 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 where c is the group velocity of the waves (Bretherton & Garrett 1968;Dewar 1970). For outward-propagating AWs in the solar wind, ω = k r (U + v A ), ω ′ = k r v A , and c = (U + v A )r, wherer is the radial unit vector and k r is the radial component of the wave vector. In a steady-state solar wind, ω depends on neither position nor time. Upon multiplying (3.14) by ω, integrating over ω, and assuming a steady state, one obtains is the total AW energy density, . . . indicates a time average, are the Elsasser variables, δv and δB are the fluctuating velocity and magnetic field, and z + (z − ) corresponds to AW fluctuations propagating away from (toward) the Sun. † In (3.17) and the following, a ± sign is used as a subscript (as opposed to a superscript) when the subscripted quantity is an r.m.s. value. As ∇·(ρUr) = 0, (3.15) can be rewritten as d dr is the wave-action flux per unit mass flux per unit ω times 4ω integrated over ω, 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 A causes partial non-WKB reflection of z + fluctuations, leading to the production of 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) 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 ⊥ /z − . The reason that only z − (and not z + ) appears in this eddy turnover time scale is that z + fluctuations are not sheared or distorted by other z + fluctuations, but they are sheared and distorted by counter-propagating z − fluctuations (Iroshnikov 1963;Kraichnan 1965). The radial decay of g 2 is accompanied by turbulent heating at the rate (3.22) † The use of ∓ on the right-hand side of (3.18) instead of ± implies that z + fluctuations propagate parallel to the background magnetic field, and z − fluctuations propagate anti-parallel to the background magnetic field. The identification of z + with outward-propagating AWs thus corresponds to the case Br > 0. If Br < 0, the same analysis goes through by replacing the ∓ on the right-hand side of (3.18) with ±. ‡ CH09 derived (3.21) starting from the MHD equations; the alternative derivation presented here is given for brevity.
which is the energy density of the outward-propagating AWs ρz 2 + /4 divided by their eddy turnover time scale L ⊥ /z − . Equation (3.22) drops a term ρz 2 − z + /(4L ⊥ ) that is normally included in the turbulent heating rate because of CH09's assumption that (3.23) an inequality that holds in the small-L ⊥ 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 ⊥ limit, obtaining Upon substituting (3.25) into (3.21), assuming that v A (r) has a single maximum at r m > r b , and solving for g(r), CH09 found that 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 The second equality in (3.29) follows from (3.23). As can be rewritten to a good approximation in the following simplified form with the aid of (3.11): (3.31)

Approximate isothermality between the transition region and Alfvén critical point
In steady state, the plasma internal-energy equation takes the form where p is the pressure, p/(γ − 1) is the internal-energy density, −p∇ · v is the rate at which p dV work is done on the plasma per unit volume, Q is the rate of turbulent heating per unit volume, q is the heat flux, which is written in the form † 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 A generic consequence of heating by reflection-driven AW turbulence is that, when other forms of heating (including conduction) are subdominant, the flow becomes quasiisothermal at r b < r < r A , meaning that whereas (3.35) is not satisfied at r > r 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 2 s . If the resulting expression for c 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 s (r) 2 is a reasonable approximation for the full solution of (3.34). On the other hand, if the resulting value of c 2 s (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 are the density and Alfvén-speed scale heights. These scale heights are generally some constant of order unity times r, with L ρ /L vA increasing by a factor of a few between the low corona and r A . On the other hand, h(r) decreases by a factor of a few between the low corona and r A , so that the product (L ρ /L vA )h(r) varies quite weakly with r. ‡ The v A /(U + v A ) term in (3.36) likewise exhibits little variation in the sub-Alfvénic region, ranging from ≃ 1 at r = r b to 0.5 at r = r A . On the other hand, ρ varies by a factor of ∼ 10 5 between the coronal base and r A (see, e.g., Cranmer et al. 2007). Thus, the approximate solution for c 2 s in (3.36) satisfies (3.35) at r b < r < r A , and AW heating indeed causes the sub-Alfvénic region to become quasi-isothermal. In contrast, at r > r A , , and the right-hand side of (3.36) becomes proportional to ρ, contradicting (3.35).

A modified version of the Chandran & Hollweg (2009) model
There are two difficulties with incorporating the CH09 model into an analytic solarwind model that includes the region immediately above the transition region. First, The left-hand side of (3.38) is the characteristic distance a z − fluctuation at scale L ⊥ (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 ⊙ (Hackenberg et al. 2000;Cranmer et al. 2007), the v A scale height has a similar value, and (3.38) is not satisfied (see, e.g., equation (3.13-a) and 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 ≃ 1.01R ⊙ , v 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 b with the expression where l b is a free parameter, which is set equal to 0.3R ⊙ in the numerical examples presented in § 4. The second difficulty is resolved by replacing (3.25) with the expression times a free parameter. † Because ρ (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) Equations (3.22), (3.29), (3.40), and (3.41) can then be used to show that is the total turbulent heating power between r b and r A , is the power (energy flux times area) of outward-propagating AWs at radius r, and is the fraction of P AW (r b ) that is dissipated between r b and r A .
3.3. Relating the AW amplitudes at the coronal base and photosphere In the absence of wave reflection and dissipation, P AW (defined in (3.44)) would have almost exactly the same value at r b and R ⊙ . † 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 AW (r b )/P AW (R ⊙ ) to some value f chr < 1, where f chr is an effective AW transmission coefficient for the chromosphere and transition region. Equivalently, because U ≪ v A at r r b , where a ⊙ subscript indicates that the subscripted quantity is evaluated at the photosphere. For the numerical calculations in § 4, I set ( 3.47) As BA = constant and v A = B/(4πρ) 1/2 , equation (3.46) can be rewritten as where (3.49) The difference between δv ⊙ and δv ⊙eff is that δv ⊙ gives rise to the fluctuating velocity δv b at r = r b when reflection and dissipation are accounted for, whereas δv ⊙eff would give rise to the same value of δv b via WKB propagation (i.e., without reflection or dissipation). The value of δv ⊙eff can be constrained in different ways. For example, observations fix δv ⊙ at a value of ≃ 1 km s −1 (Richardson & Schwarzschild 1950), and numerical simulations suggest that f chr is ≃ 0.04−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 δv ⊙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 δv ⊙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 δv ⊙ , f chr , and δv ⊙eff , which are plausibly independent of the properties of the coronal plasma and coronal magnetic field, δv b depends upon ρ b , which varies between different flux tubes with different super-radial expansion factors η b , as shown later in (3.55).

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 where m p is the proton mass, and Λ(T ) is the optically thin radiative loss function. Figure 2 shows three different approximations to Λ(T ). The dashed lines correspond to equation A1 of Rosner et al. (1978). The solid line is a plot of Rosner et al. (1978). The dotted line in figure 2 is a piecewise-continuous linear approximation to the lowtemperature range of Λ(T ) in figure 1 of Cranmer et al. (2007), which is included to illustrate that optically thin radiative cooling becomes extremely weak at T 10 4 K. Throughout most of the corona, ρ is sufficiently small that radiative cooling is negligible. However, as r decreases within the corona, ρ increases by orders of magnitude, which causes R/Q to increase, because R ∝ ρ 2 . There is thus some radius r b at which 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 b is identified as the base of the corona, as already noted in (3.1). As r decreases below r b , ρ(r) increases above ρ b , R/Q increases to values ≫ 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 dV 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 b and radiative cooling at r > r b .
With the aid of (3.39), (3.46), and (3.50), equation (3.53) can be written in the form (3.54) Solving for ρ b , one finds via (3.51) that is the sound speed (2.6) in the sub-Alfvénic region (r b < r < r A ), which is approximated as a constant for the reasons discussed in § 3.2.2. Equation (3.55) can be rewritten as is the dimensionless temperature of the sub-Alfvénic region, is the magnetic-field strength for which v A = v esc when the radiative cooling time ρc s 2 /R, free-fall time R ⊙ /v esc , and sound-crossing time R ⊙ /c s are all equal, (3.61) (for reference,ρ ⊙ = 5.7 × 10 5 given (3.47)), and and (3.48), (3.57), and (3.60) can be used to write 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 b and r = r A yields where (3.42)-(3.44) have been used to rewrite the volume integral of the turbulent heating rate, equation (3.8) has been used to pull a constant factor ofṀ = ρAU out of the integral on the left-hand side, and (3.11) has been used to write The heat flux at r = r A is a small fraction of the total energy flux in AW-driven solarwind models and is thus neglected in (3.65). Upon evaluating the integral on the left-hand side of (3.65), noting that q rb = −|q rb | ≡ −q b , making use of (3.12), and rearranging terms, one obtains 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 b and r 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 dV work and thermal conduction into the transition region. Dividing (3.66) byṀ v 2 esc leads to ( 3.67) 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 K (the electron Coulomb mean free path λ 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 K ≪ 1, the radial component of the heat flux in the transition region is well approximated by the Spitzer & Härm (1953) formula, the value for electron-electron collisions in a proton-electron plasma with ρ/m p = 10 9 cm −3 and T = 10 6 K (Huba 2013). The magnetic field near r = r b is taken to be approximately radial, and the width of the transition region (∼ 10 2 km) is so narrow that p, B, and A are treated as constants within the transition region, with As mentioned previously, turbulent heating is neglected at r < r b . Within the transition region, equation (3.32) thus becomes where ρ 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) where the first equality follows from (3.6). With the aid of (3.73), and using (3.68) to write which can be integrated from the chromospheric values of T and w, denoted T chr and w chr , to the values of T and w at r = r b , denoted T b and w b . As T chr ≪ T b and the chromospheric heat flux is negligible (q scaling as T 7/2 divided by the temperaturegradient scale length), the chromospheric terms are dropped, and the integral becomes where (3.80) 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 .
If U b /c s is sufficiently small, then a 3 ≪ 1. In this low-Mach-number limit, (3.82) becomes, to leading order in a 3 , (3.83) 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 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 dV 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 b and does not assume that q(T )dT ∝ 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 where ρz 2 + /8 is the AW pressure, which is one-half the AW energy density (Dewar 1970 and η 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 c (the wave-modified sonic critical point). In order for dρ/dr to remain finite at r = r c , the right-hand side of (3.86) must also vanish at r c . Together, these conditions yield the constraints † Schwadron & McComas (2003) integrate the internal-energy equation from the upper chromosphere all the way out to the coronal temperature maximum, where the heat flux vanishes. The value of q b in (3.82) corresponds to r = r b , at which turbulent heating and radiative cooling balance, which, in the low-Mach-number limit, corresponds to the maximum of the heat flux. This heat-flux maximum is lower down in the solar atmosphere than the temperature maximum. and where, here and in the following, a 'c' subscript indicates that the subscripted quantity is evaluated at r = r c . An implicit assumption underlying (3.88) and (3.89) is that so that the quasi-isothermal approximation applies at r c . Mass and flux conservation (i.e., (3.7)) imply that where the second equality in (3.91) follows from (3.11). Equation (3.2) implies that where the second equality in (3.92) follows from (3.88). Upon substituting (3.92) into (3.91) and evaluating v Ab using (3.63), one obtains Substituting (3.93) into (3.89) yields The integral over r of ρ −1 times the momentum equation (3.84) yields the Bernoulli integral, (3.95) where Γ is independent of r. Evaluating (3.95) at r = r b leads to the equality . (3.96) Evaluating (3.95) at r = r c , rewriting r c and U c using (3.88) and (3.89), rewriting Γ using (3.96), and multiplying the resulting equation by 2/v 2 esc yields 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 ⊙ , M ⊙ , v esc , ρ ⊙ ,ρ ⊙ , ψ, δv ⊙eff ,B); (2) free parameters (σ, l b ); (3) the super-radial expansion factor η(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 b , y c , and x; and (5) additional unknowns that can be determined once y b , y c , and x are found (Ṁ , χ H , U (r), ρ b , q b , v Ab , r c , r A ). The three principal unknowns y b , y c , and x are determined by solving the three simultaneous equations (3.67), (3.94), and (3.97), where it must be remembered that ǫ is itself a function of x via (3.64). The additional unknownsṀ , χ H , ρ b , v Ab , q b , and r c then follow immediately from (3.13), (3.45), (3.57), (3.63), (3.82), and (3.88), respectively. For example,Ṁ The procedures for determining r 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 The last equality in (3.100) amounts to taking the magnetic-field lines to be purely radial at r r c . The first of the two approximate analytic solutions is valid in the conductiondominated 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 dV work. As discussed further in the following, this regime arises only for values of δv ⊙eff much smaller than the solar value. This limit is thus not directly relevant to solar-wind observations. To leading order in ǫ ⊙ , the mass outflow rateṀ (cond) and asymptotic flow velocity U (cond) ∞ in this parameter regime are given bẏ where I 2 is a numerical constant given in (B 12), and 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 dV 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Ṁ where P AW (r 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 and the coronal temperature in the expansion-dominated regime is to leading order . (3.105) 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, δv 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 higherorder 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 AW (r b ) is dissipated within the sub-Alfvénic region. Second-and fourth-order approximations toṀ and U ∞ are shown in figure 5 below. Analytic estimates for the ranges of ǫ ⊙ values corresponding to the conductiondominated and expansion-dominated limits are given in Appendix B. There are three constraints on ǫ ⊙ 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 ǫ ⊙ values is much smaller than the solar value, as illustrated in figure 4. The expansion-dominated limit corresponds to a finite range of ǫ ⊙ values that is sufficiently large that p dV 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 ǫ ⊙ values is relevant to the solar case, as shown in the next section.

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 b , y c , and x using Newton's method. † Once y b , y c , and x are determined,Ṁ , χ H , ρ b , v Ab , q b , r c , r A , and U ∞ are computed from (3.13),

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 η(r). On the other hand, the approximate analytic solutions derived in Appendix B assume that η c = γ Bc = 1. (4.1) 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, equation (4.1) amounts to assuming that all of the super-radial expansion of the magnetic field occurs inside the wave-modified sonic critical point.
4.3. Free parameters As discussed in § 3.8, there are two free parameters in the model: l b (the AW dissipation length scale at r b ) and σ (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 b . Throughout the rest of this section, l b is thus simply fixed at a value that seems physically reasonable: (4.7) On the other hand, the solutions depend sensitively on the value of σ. Larger values of σ cause a larger fraction of the AW power at the coronal base P AWb to be dissipated in the quasi-isothermal sub-Alfvénic region, which increasesṀ (see, e.g., (2.15) and (2.16)). For fixed P AWb , increasingṀ decreases U ∞ . In some of the examples to follow, σ 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, η b , also varies in the model, but this quantity is in principle observable for individual flux tubes. The value of η b is thus treated as an input into the model rather than an adjustable parameter. 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 η must be specified at that r. For figure 3, I set η(r) = 1 + (η b − 1) exp(−(r − r b ) 2 /(0.02R ⊙ ) 2 ), which is effectively unity at r = r c , consistent with (4.1). This choice of η(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 δv at r = r b (i.e., δv b ) shown in the lower-left panel of Figure 3 depends on η b through (3.48), (3.57), and (3.62), δv b does not on the way in which η(r) decreases from η b to 1. The solution in Figure 3 is based upon the somewhat arbitrary assumption that η b = 100 in the magnetic flux tube encountered by PSP at the time of its first perihelion. The value of σ is set equal to 0.5 to optimise the agreement between the model and the data. 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 . The fact that all four quantities plotted in Figure 3 can be matched by varying the single parameter σ suggests that the model is reasonably successful at capturing the physical processes that control the heating and acceleration of coronal-hole outflows.

Fiducial solution matching measurements from Parker Solar Probe's first perihelion
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 Γ . One of the dotted lines intersects the solid line, and that intersection occurs at the wave-modified sonic critical point, r = r c . That dotted line is an accretion-like solution of (3.95) with the same value of Γ as in the model solution but with dU/dr < 0 at r = r c .  in (B 14). The conduction-dominated approximation that gives rise toṀ (cond) in (3.101) assumes that ǫ ⊙ ≪ ǫ ⊙cond . † The shaded region corresponds to ǫ ⊙exp,min < ǫ ⊙ < ǫ ⊙exp,max , where ǫ ⊙exp,min and ǫ ⊙exp,max are defined in (B 35) and (B 36). The expansion-dominated approximation in § B.2 assumes that ǫ ⊙exp,min ≪ ǫ ⊙ ≪ ǫ ⊙exp,max . The vertical long-dashed line corresponds to the condition r c = r A . To the right of this line,Ṁ and the solar-wind density become so large † It should be noted, however, that the solution for U (cond) ∞ in (3.102) exceeds the speed of light when δv ⊙eff 2 × 10 −3 km s −1 , and a relativistic treatment would be needed to model the outflow correctly in this limit. that r A < r 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 δv ⊙eff much smaller than the solar value. 4.6. Anti-correlation between the coronal super-radial expansion factor and U ∞ Wang & Sheeley (1990) showed that the outflow velocity in a solar-wind stream at r = 1 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 max,WS to be η(1.04R ⊙ )/η(2.5R ⊙ ). I then compute η(1.04R ⊙ ) and η(2.5R ⊙ ) 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 et al. (1998)  I take (4.8) to hold even as η b is varied in Figures 5, 6, and 7. This is equivalent to assuming that if local magnetic structures on the Sun cause η b to increase or decrease by some factor relative to the Hackenberg/Banaszkiewicz model just described, then η(1.04R ⊙ ) increases or decreases by the same factor, but η(2.5R ⊙ ) does not change.

Illustration of the conduction-dominated and expansion-dominated regimes
If σ is fixed while η 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 σ is taken to have a power-law dependence on η b of the form σ = 5.7 × 10 −2 η 0.34 b .
(4.9) Equation (4.9) is used to determine σ in figures 5 through 7. The implication of (4.9) that σ is an increasing function of η b is plausible because increasing η b increases v 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 ( , which increases the fraction of P AWb that is dissipated at r < r 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 η 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 max,WS ranges from ≃ 3 to ≃ 30,Ṁ varies from 10 −14 M ⊙ yr −1 to 2 × 10 −14 M ⊙ yr −1 , values that are similar to the solarmass-loss rates inferred from Ulysses and PSP measurements (McComas et al. 2000;Schwadron & McComas 2008;Kasper et al. 2019). The fourth-order analytic approxi-mationṀ is also reasonably accurate at f max,WS 3, but deviates markedly from the numerical solution at f max,WS 2. A similar comment applies to the top-right panel of figure 5, which also shows ∞,4 , 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 fmax,WS data bins in that table. The vertical error bars correspond to one-half of the 100 km s −1 increment between the discretised U∞ values that define four of the five data bins. The quantity I1ρ b cs 3 in the left panel of the second row is the low-Mach-number approximation to q b given in (3.83).
that the model agrees fairly well with observational constraints on U ∞ (f max,WS ) when (4.9) holds. The left panel of the second row of figure 5 shows that the low-Machnumber approximation to q b given in (3.83) is quite accurate at small f max,WS , where U b /c s ≪ 1 (see the right panel in the second row of this figure). However, as f max,WS increases, U b increases, because the outflow is concentrated into a smaller cross-sectional area at the coronal base. This increase in U b leads to larger losses of internal energy within the transition region from p dV work and advection, which, in turn, causes q 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 b is significantly smaller than the AW energy flux at the coronal base , which increases with f max,WS because increasing f max,WS increases η b and hence B b . The lower-left panel of figure 5 shows that ρ b increases by a factor of ≃ 10 as f max,WS increases from 1 to 100, consistent with the ρ b ∝ η 1/2 b x 1/4 scaling in (3.57) and (3.62), given that x = c s 2 /v 2 esc is approximately constant over this range of f max . The lower-right panel shows that r c is typically several solar radii, whereas r A ranges from 12R ⊙ to 15R ⊙ 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 max,WS and δv ⊙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Ṁ is an increasing function of both f max,WS and δv ⊙eff , whereas U ∞ is a decreasing function of both f max,WS and δv ⊙eff . The second row shows that r A is a strongly decreasing function of δv ⊙eff and only weakly dependent on f max,WS for values of δv ⊙eff comparable to the Sun-like value in (4.6). The right panel of this row shows that r c is an increasing function of f max,WS and, for the most part, a decreasing function of δv ⊙eff . The lower-left panel of figure 6 shows that c s varies only weakly with f max,WS and δv ⊙eff . As shown in the lower-right panel, χ H varies by a factor of ≃ 2 over the parameter range shown. If δv ⊙eff is set equal to the value in (4.6) and f max is restricted to the interval (2, 8) so that U ∞ in the upper-right panel takes on fast-wind-like values of 600-800 km s −1 , then χ H ≃ 0.5 − 0.7, consistent with direct numerical simulations of reflection-driven AW turbulence (Perez et al. in press). Figure 7 displays contour plots of the same quantities as in figure 6, but this time as functions of f max,WS and B * , or, equivalently, B r (1 a.u.), using the parameters in (4.1), (4.6), (4.7), and (4.9). The top panels show thatṀ increases approximately linearly with B r (1 a.u.) at fixed f max,WS , and that U ∞ varies only weakly with B r (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 A depends more strongly on B r (1 a.u.) than on f max,WS , whereas the reverse is true for r c . As in figure 6, c s varies very weakly across the entire parameter range observed. The lower-right panel shows that, when (4.9) holds, χ H depends more strongly on f max,WS than on B r (1 a.u.).

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Ṁ , U ∞ , 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 dV 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 b and Alfvén critical point r A ; and (3) the finding in recent direct numerical simulations that most of the AW power at the coronal base P AW (r b ) is dissipated between r = r b and r = 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 soniccritical point, and the fact that part of the AW power reaches the super-Alfvénic region, where it enhances U ∞ 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: • the plasma density at the coronal base ρ(r b ) = ρ b is determined by equating the turbulent heating rate Q(r b ) and radiative cooling rate R(r b ); and • an analytic solution for the heat flux from the corona into the transition region q b is obtained by balancing, within the transition region, conductive heating against internal-energy losses from p dV work, advection, and radiative cooling.
The expression (3.57) for ρ b that results from setting Q(r b ) = R(r b ) contains a factor of x 1/4 , where x = c s 2 /v 2 esc is the dimensionless temperature. Thus, x must be found before the exact value of ρ b can be determined. However, becauseṀ 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 ρ 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 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 dV 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 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 δv 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Ṁ , q 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(Parker , 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 c , several R ⊙ 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 c from flowing outward at approximately the wave-enhanced effective sound speed, c s,eff ≡ [(p + p wave )/ρ] 1/2 = [c 2 s + 0.5(δv) 2 ] 1/2 , evaluated at r = r c , which is comparable to the square root of the right-hand side of (3.89). This is whyṀ ∼ A(r c )ρ c c s,eff (r c ).
On the other hand, localised motions near the transition region at speeds ≪ v esc are gravitationally bound, and the mass flux that they carry is not determinative ofṀ . For example, if at some time, such motions led to an overall mass outflow rate at r b exceeding the rate ∼ A(r c )ρ c c s,eff (r c ) at which mass flows past the critical point r 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Ṁ in two ways. First, q b reduces the net heating power within the quasi-isothermal sub-Alfvénic region, P net . As discussed following (2.16), the heating cost per unit mass for plasma to transit the quasiisothermal sub-Alfvénic region is c s 2 ln(ρ b /ρ(r A )) ≃ v 2 esc , andṀ ≃ P net /v 2 esc . By reducing P net , the conduction of heat from the corona into the transition region reducesṀ . The second way that the lower solar atmosphere affectsṀ is via the chromospheric/transitionregion transmission coefficient, f chr = P AW (r b )/P AW (R ⊙ ), whose value again influences the value of P net . To summarise this paragraph and the preceding paragraph, the chromosphere and transition region influenceṀ thermodynamically, but not dynamically.
The regulation of the mass flux at r b by the critical-point condition at r 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 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, ρ b is determined by the condition that Q(r b ) = R(r b ), without reference to conditions in the chromosphere or the value of q b . The sound speed at r = r b , which is ≃ c s , is approximately determined by balancing turbulent heating against internal-energy losses from p dV work within the corona and sub-Alfvénic solar wind. The outflow velocity at r = r b , U b , follows from the value ofṀ , which is controlled by the critical-point condition, as described above. The values of ρ b , c s , and U b jointly determine q b via (3.82), which embodies the requirement that conductive heating offset internal-energy losses (from radiation, p dV work, and advection) within the transition region. The values of ρ b and c s are sufficient to determine the approximate transitionregion pressure, p tr ≃ ρ b c s 2 . The pressure within the upper chromosphere, p chr (r), is an approximately exponentially decreasing function of altitude. The altitude of the transition region is determined by setting p chr (r) = p 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 b ), which is controlled by the balance between heating and p dV work in the corona, as discussed previously.
Although the condition R(r b ) = Q(r b ) determines the density ρ 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 Λ(T ), as shown in Figure 2. The density at the bottom of the transition region is instead largely determined by the values of ρ b and c 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 ∼ 10 4 K in the upper chromosphere.

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 p c s 2 /(2k B ) of the sub-Alfvénic region increases from ≃ 7×10 5 K when U ∞ ≃ 800 km s −1 to ≃ 8.4×10 5 K when U ∞ ≃ 400 km s −1 . In contrast, in measurements from the Ulysses spacecraft, the coronal freeze-in temperature increases from ≃ 9 × 10 5 K when U ∞ ≃ 800 km s −1 to ≃ 1.35 × 10 6 K when U ∞ ≃ 400 km 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énicregion 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 σ 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 σ 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 superradial 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. Appendix A. Solving for U (r) and r A Once y b , y c , and x are determined by solving (3.67), (3.94), and (3.97), the value of U (r) between r = r b and r = r A (which is as yet unknown) can be found by solving the Bernoulli integral (3.95) with Γ given by (3.96). An equation for r A can be obtained by evaluating the Bernoulli integral (3.95) at r = r A , setting U (r A ) = v A (r A ) and y = 1, and rewriting Γ using (3.96), which leads to in (A 1) and rewriting v Ab in (A 1) using (3.63), one obtains the following equation for r A : At r > r A , the outflow velocity U (r) cannot be determined from the Bernoulli integral, because the quasi-isothermal approximation does not apply. In addition, equation (3.40) yields a poor approximation for z − at r > r A , because 1 + y approaches a constant when r ≫ r A . A better approximation for z − in the super-Alfvénic region can be obtained from (3.25) and the simplifying approximation that |(d/dr) ln(v A )| = −(d/dr) ln y, which holds when B ∝ r −2 and ρ ∝ r −2 (the latter scalings being fairly accurate for 10R ⊙ r 60R ⊙ , a region in which the field lines are approximately radial and the flow speed is approximately constant). In this case, equation (3.25) becomes Integrating (3.21) from r = r A out to larger r then gives z 2 + = [z + (r A )] 2 4y 2 /(1 + y) 2 , where the value of z + (r A ) is obtained by setting y = 1 in (3.41). The approximate value of c 2 s (r) can then be found in terms of ρ by solving the internal-energy equation (3.32) with Q given by (3.22) and neglecting radiative cooling and thermal conduction. This leads to Total-energy conservation implies that the mechanical luminosity is independent of r, where the ratio of specific heats γ has been set equal to 5/3. The value of U (r) at r > r A can be obtained by setting and rewriting ρ(r) in terms of U (r) using (3.7); i.e., ρ(r)U (r)/B(r) Determining U (r) at r > r A via (A 7) through (A 9) requires evaluating r A and v A (r A ). An alternative method that avoids this requirement results from noting that ) where (3.66) has been used to eliminate q b . Replacing χ H in (A 10) with the expression on the right-hand side of (3.45) leads to the consistency check that L mech (r A ) = L mech (r b ).
(A 11) Combining (A 8) and (A 11), one can find U (r) at r > r A by setting L mech (r) = L mech (r b ).
(A 12) Equation (A 12) leads to a simple expression for the asymptotic wind velocity U ∞ , i.e., U (r) as r → ∞. As r → ∞, the kinetic-energy flux dominates the total energy flux, and

Appendix B. Approximate analytic solutions
As discussed in § 3.8, the core of the solar-wind model developed in § 3 is a set of three simultaneous equations, (3.67), (3.94), and (3.97), for the three unknowns y b , y c , and x. In this section, two different approximate analytic solutions to these equations are obtained in two different parameter regimes. Both solutions rely on the simplifying approximations in (3.100), which are repeated here: In particular, y b is taken to be sufficiently large that: (1) the U 2 b = v 2 Ab /y 2 b term in (3.96) can be dropped, which amounts to dropping the second-to-last term on the left-hand side of (3.97); and (2) O(y −1 b ) terms can be dropped in (3.45), so that The y −σ b term in (B 2) is retained, even though terms of order y −1 b are discarded, on the working assumption that σ ∼ 0.1-0.5, as is the case in the numerical examples in § 4. As mentioned in § 3.8, the last equality in (B 1) amounts to taking all of the super-radial expansion of the field lines to occur inside the wave-modified sonic critical point and the field lines to be completely radial at r r c . With these approximations, (3.67), (3.94), and (3.97) become, respectively, and B.1.

Conduction-dominated limit
When ǫ ⊙ is sufficiently small, an approximate solution to (B 4), (B 5), and (B 6) can be obtained through the method of dominant balance (Bender & Orszag 1978), in which two of the three terms in each equation are taken to be dominant, and the third term is taken to be much smaller in magnitude. Neglecting the smaller term in each equation yields the leading-order solution, with the smaller term producing higher-order corrections. In the present case, it is the second term on the left-hand side of each equation that can be neglected to leading order. In (B 4), this corresponds to balancing turbulent heating of the sub-Alfvénic region against the heat that is conducted from the corona into the transition region. In other words, conduction into the transition region rather than p dV work is the dominant sink of internal energy in the sub-Alfvénic region as a whole. Neglecting the second term in (B 5) amounts to assuming that the fluctuating velocity makes the dominant contribution to U c in (3.89), which is equivalent to taking δv(r c ) ≫ c s . The latter equality appears paradoxical, because the small-ǫ ⊙ limit corresponds to small values of the fluctuating velocity at the coronal base. However, as ǫ ⊙ → 0, the coronal temperature drops, the density scale height in the corona decreases, and the wave amplitude at the critical point grows as the AWs attempt to conserve wave action, which would lead to δv ∝ ρ −1/4 when U ≪ v A in the absence of reflection and dissipation (see (3.19) and (3.20)). Neglecting the second term in (B 6) amounts to taking the wave pressure force to have a larger cumulative effect than the thermal pressure force on the acceleration of the outflowing plasma between r b and r c , which can again be understood as a consequence of the wave amplitudes growing rapidly with increasing r when the density scale height in the corona is small.
As already noted in (3.101) and illustrated in Figure 4,Ṁ is exceedingly small in the conduction-dominated limit, and as a consequence U b ≪ c s . The heat flux at the coronal base is thus approximately given by (3.83), and (B 4) can be rewritten as Balancing the first and third terms in (B 7) and dropping the y −σ b term in (B 2) yields the leading-order solution for the dimensionless temperature in the sub-Alfvénic region: Anticipating the solution for y c , it is useful to predict at the outset (as will shortly be confirmed by (B 10) and (B 11)) that y c ≫ 1. .
(B 12) Equations (3.13), (3.57), (B 8), and (B 11) then yield the leading-order mass outflow rate in the conduction-dominated regime,Ṁ (cond) , given in (3.101). The asymptotic wind velocity U ∞ is obtained by settingṀ U 2 ∞ /2 equal to the mechanical luminosity at the coronal base L mech (r b ) as described in Appendix A, where L mech is defined in (A 7). When this procedure is carried out using (B 8) and (B 11), the two leading-order terms in L mech (r b ) cancel. To obtain the leading-order non-vanishing term in U 2 ∞ , one must account for the next largest term in (B 7), which results from the χ H correction; i.e., the second term on the right-hand side of (B 2). When this term is retained, one obtains the leading-order asymptotic wind speed in the conductiondominated regime, U (cond) ∞ , given in (3.102). The range of ǫ ⊙ values for which (3.101), (3.102), (B 8), (B 10), and (B 11) are approximately valid can be determined by requiring that the neglected terms in (B 5), (B 6), and (B 7) be small compared to the other terms when y b , y c , and x are given by (B 8), (B 10), and (B 11). The most stringent condition on ǫ ⊙ arises from carrying out this procedure for (B 6), which results in the requirement that where I 3 = I −4/7 1ρ 1/7 ⊙ (η b B * lb ) 2/7 , and I 4 equals the right-hand side of (B 10) without the ǫ ⊙ term; i.e., I 4 = ǫ 5/(7−7σ) ⊙ y b /y c , with y b /y c given by the right-hand side of (B 10). Upon neglecting the quantity 4 ln I 4 + 3 on the left-hand side of (B 13), which is smaller in magnitude than the remaining term when ǫ ⊙ is in the conduction-dominated regime but other parameters take on Sun-like values, one finds that the conduction-dominated regime corresponds to where W −1 is the lower branch of the Lambert W function.

B.2. Expansion-dominated limit
In the expansion-dominated limit, the first two terms on the left-hand sides of (B 4), (B 5), and (B 6) are treated as dominant. For this case, it is useful to define the variables The latter inequality is achieved when −1/ ln A 1 , |a|, b, and c are much smaller than 1. When (B 28) is satisfied, equations (B 23), (B 24), and (B 25) can be solved perturbatively through the recursion relations 1 u n − 1 = 0 if n = 0 S(u n−1 , a n−1 , b n−1 , c n−1 ) if n 1 (B 29) x n = (1 − a n−1 )u n −2 ln A 1 (B 30) where n = 0, 1, 2, . . . , and a −1 = c −1 = 0. The values of a n , b n , and c n are obtained by replacing (a, b, c, u, y b , y c , x) with (a n , b n , c n , u n , y b,n , y c,n , x n ) in (B 20), (B 21), and (B 22). For reference below, the values of ρ b , v Ab , and q b that result from this substitution (via (3.57), (3.63), and (3.82)) are denoted ρ b,n , v Ab,n , and q b,n . The values of y b,n and y c,n are obtained by replacing (y b , y c , x, u, p) with (y b,n , y c,n , x n , u n , p n ) in (B 15). The n th -order approximation forṀ in the expansion-dominated regime, denoteḋ M (exp) n , can be found by settingṀ =Ṁ For n 1, I define the n th -order approximation for U ∞ , denoted U (exp) ∞,n , to be the value of U ∞ in (A 13) when (x, y b , ψ, q b , ρ b , v Ab ) = (x n , y b,n , 1, q b,n , ρ b,n , v Ab,n ): where I have invoked (B 1) to approximate 3/2 + y b,n as y b,n . The leading-order approximation for U ∞ in the expansion-dominated regime, denoted U ∞,0 , is obtained from (B 33) with n = 0 after dropping terms that were neglected in the calculation of u 0 -in particular, the first, second, and last terms inside the brackets on the right-hand side of (B 33). This yields U (exp) ∞,0 = v esc .
(B 34) As in the conduction-dominated limit, the range of ǫ ⊙ values for which (B 29), (B 32), and (B 34) are approximately valid can be determined by imposing the constraint that the neglected terms in (B 4), (B 5), and (B 6) be small compared to the terms that are kept when y b , y c , and x are given by (B 15), (B 16), (B 23), (B 24), and u = 1. Carrying out this procedure for (B 4) and making the simplifying approximations that a is dominated by the last term on the right-hand side of (B 20), that q b is given by (3.83), and that ln A 1 ≃ ln ǫ 3/4 ⊙ , one obtains the requirement that ǫ ⊙ ≫ ǫ ⊙exp,min ≡ exp 7 2 W −1 − 4I where W −1 is, as above, the lower branch of the Lambert W function. Equation (B 35) corresponds to the requirement that the conductive losses from the sub-Alfvénic region into the transition region be negligible compared to p dV work in this region. Carrying out the above procedure for (B 5) and again making the simplifying approximation that ln A 1 ≃ ln ǫ (B 36) Equation (B 36) corresponds to the requirement that the sound speed make the dominant contribution to the outflow velocity U c at the critical point in (3.89). As illustrated by the shaded gray rectangle in figure (4), for Sun-like parameters (and in particular, for η b = 30), ǫ ⊙exp,max is approximately four orders of magnitude larger than ǫ ⊙exp,min . There is thus a finite range of ǫ ⊙ values that satisfy both (B 35) and (B 36). However, it should be noted that the approximations ln A 1 ≃ ln(ǫ 3/4 ⊙ ) and q b = I 1 ρ b c s 3 cause (B 35) to underestimate the lower bound on ǫ ⊙ in the expansion-dominated regime. Also, for Sun-like parameters, the assumption r c < r A that underlies the model of § 3 breaks down at values of ǫ ⊙ smaller than ǫ ⊙exp,max , as illustrated in Figure 4.