TESTING A PREDICTIVE THEORETICAL MODEL FOR THE MASS LOSS RATES OF COOL STARS

We begin with five fundamental parameters that are assumed to determine (nearly) all of the other relevant properties of a star: mass M_*, radius R_*, bolometric luminosity L_*, rotation period P_rot, and metallicity. We also assume that the star's iron abundance, expressed logarithmically with respect to hydrogen, is a good enough proxy for the abundances of other elements heavier than helium, i.e., [Fe/H] ≈log (Z/Z_☉). For spherical stars, the effective temperature T_eff and surface gravity g are calculated straightforwardly from

$$\begin{equation} \sigma T_{\rm eff}^{4} \, = \, \frac{L_{\ast }}{4\pi R_{\ast }^2} \,, \,\,\,\, g \, = \, \frac{G M_{\ast }}{R_{\ast }^2}\,, \end{equation} \tag{ 1 }$$

where σ is the Stefan–Boltzmann constant and G is the Newtonian gravitation constant.

We need to know the mass density in the stellar photosphere ρ_* in order to specify the properties of MHD waves at that height. The density was computed from the criterion such that the Rosseland mean optical depth should have a value of 2/3 in the photosphere. We used the AESOPUS opacity database,¹ which is a tabulation of the Rosseland mean opacity κ_R as a function of temperature, density, and metallicity (see Marigo & Aringer 2009). An approximate expression for the Rosseland optical depth τ_R in the photosphere,

$$\begin{equation} \tau _{\rm R} \, = \, \kappa _{\rm R} \rho _{\ast } H_{\ast } \, = \, 2/3, \end{equation} \tag{ 2 }$$

was solved for ρ_*, where H_* is the photospheric value of the density scale height (see below). We used straightforward linear interpolation to locate the relevant solutions for ρ_* as a function of T_eff, g, and [Fe/H].

Figure 1(a) shows how the photospheric density varies as a function of T_eff and log  g under the assumption of solar metallicity ([Fe/H] =0). For context we also show the location of the zero-age main sequence (ZAMS) from the models of Girardi et al. (2000), as well as a post-main-sequence evolutionary track for a 1 M_☉ star from the BaSTI² model database (Pietrinferni et al. 2004).

Zoom In Zoom Out Reset image size

We used the 2005 release of OPAL plasma equations of state³ (see also Rogers & Nayfonov 2002) to estimate the mean atomic weight μ in a partially ionized photosphere. For the range of parameters appropriate for cool stars, we found that μ is primarily sensitive to T_eff, and not to gravity or metallicity, so we produced a single parameter fit:

$$\begin{equation} \mu \, \approx \, \frac{7}{4} + \frac{1}{2} \tanh \left(\frac{3500 - T_{\rm eff}}{600} \right), \end{equation} \tag{ 3 }$$

where T_eff is expressed in K. Other quantities that will be needed later include the photospheric density scale height, which is given by

$$\begin{equation} H_{\ast } \, = \, \frac{k_{\rm B} T_{\rm eff}}{\mu m_{\rm H} g}, \end{equation} \tag{ 4 }$$

where k_B is Boltzmann's constant and m_H is the mass of a hydrogen atom. We also need to compute the equipartition magnetic field strength,

$$\begin{equation} B_{\rm eq} \, = \, \sqrt{8\pi P_{\ast }} \, = \, \sqrt{\frac{8\pi \rho _{\ast } k_{\rm B} T_{\rm eff}}{\mu m_{\rm H}}}, \end{equation} \tag{ 5 }$$

where P_* is the photospheric gas pressure. Because T_eff and μ do not vary over many orders of magnitude, it is roughly the case that B_eq∝ρ^1/2_*. However, in all calculations below we compute B_eq fully from Equation (5). In Section 4, we describe observations that show the photospheric magnetic field strength B_* is roughly linearly proportional to B_eq for many stars. The measurements determine the constant of proportionality, and we use

$$\begin{equation} B_{\ast } \, = \, 1.13 \, B_{\rm eq} \end{equation} \tag{ 6 }$$

in the remainder of this paper. Figure 1(b) shows B_* as a function of T_eff and log  g.

We consider MHD waves that are driven by turbulent convective motions in the stellar interior. The original models of wave generation from turbulence (e.g., Lighthill 1952; Proudman 1952; Stein 1967) dealt mainly with acoustic waves in an unmagnetized medium. More recently, however, it has been shown that when a stellar atmosphere is filled with magnetic flux tubes, the dominant carrier of wave energy should be transverse kink-mode oscillations (Musielak & Ulmschneider 2002a). When the magnetic flux tubes extend above the stellar surface and expand to fill the volume, the kink-mode waves become shear Alfvén waves (see Cranmer & van Ballegooijen 2005).

We utilize the model results of Musielak & Ulmschneider (2002a) to estimate the flux of energy in kink/Alfvén waves in stellar photospheres. For simplicity, we used only the simulations of Musielak & Ulmschneider (2002a) with their standard parameter choices: a mixing length parameter of α = 2 and a constant magnetic field strength that is 0.85 times the equipartition field strength. Our analytic fit to the results shown in their Figure 8 is

$$\begin{equation} F_{\rm A \ast } \, = \, F_{0} \left(\frac{T_{\rm eff}}{T_0} \right)^{\varepsilon } \exp \left[ - \left(\frac{T_{\rm eff}}{T_0} \right)^{25} \right], \end{equation} \tag{ 7 }$$

where the dependence on $\tilde{g} = \log \, g$ is given by

$$\begin{equation} \frac{F_0}{10^{9} \,\, \mbox{erg} \,\,\, \mbox{cm}^{-2} \,\,\, \mbox{s}^{-1}} \, = \, 5.724 \,\, \exp \left(- \frac{\tilde{g}}{11.48} \right), \end{equation} \tag{ 8 }$$

$$\begin{equation} \frac{T_0}{1000 \,\, \mbox{K}} \, = \, 5.624 + 0.6002 \, \tilde{g}, \end{equation} \tag{ 9 }$$

$$\begin{equation} \varepsilon \, = \, 6.774 + 0.5057 \, \tilde{g}. \end{equation} \tag{ 10 }$$

These fits are similar in form to those given by Fawzy & Cuntz (2011) for longitudinal MHD waves. Figure 2 shows a comparison between the above fitting formula and the plotted results of Musielak & Ulmschneider (2002a) for log  g = 3, 4, and 5. The behavior of F_A* for lower values of log  g was not given by Musielak & Ulmschneider (2002a), but similar results were found for a wider range of gravities by Ulmschneider et al. (1996) for acoustic waves.

Zoom In Zoom Out Reset image size

We used the kink-mode energy flux to determine the transverse velocity amplitude v_⊥ of Alfvén waves in the photosphere. The flux is defined as

$$\begin{equation} F_{\rm A \ast } \, = \, \rho _{\ast } v_{\perp \ast }^{2} V_{{\rm A} \ast } \end{equation} \tag{ 11 }$$

with V_A* = B_*/(4πρ_*)^1/2 being the photospheric Alfvén speed. The above expression is not exact for waves undergoing strong reflection (see, e.g., Heinemann & Olbert 1980), but it ends up giving a similar prediction for the height variation of v_⊥ in the corona that would come from a more accurate non-WKB model (Cranmer & van Ballegooijen 2005). Figure 1(c) shows how v_⊥* varies as a function of T_eff and log  g for solar metallicity stars.

For the well-observed case of the Sun, we know that most of the photospheric magnetic field is concentrated into small (100–200 km diameter) flux tubes concentrated in the intergranular downflow lanes (Solanki 1993; Berger & Title 2001). The field strength in these tubes is close to equipartition, with B_* ≈ 1400 G. However, these flux tubes have a filling factor f_* in the photosphere of about 0.1% to 1%, so the spatially averaged magnetic flux density B_*f_* is only of order 1–10 G (Schrijver & Harvey 1989).

Figure 3 illustrates the stellar magnetic field geometry that we assume to exist above the surface of a cool star. Flux tubes that are open to the stellar wind⁴ have a cross-sectional area A(r) that expands monotonically with increasing radial distance r from the star. The condition ∇ · B = 0 demands that the product of A and the magnetic field strength B remains constant. Thus, B(r) inside a flux tube decreases monotonically, from its photospheric value of B_*, with increasing distance. We normalize A such that at a given distance the total stellar surface area covered by open flux tubes is defined to be

$$\begin{equation} A \, = \, 4 \pi r^{2} f. \end{equation} \tag{ 12 }$$

The dimensionless filling factor f tends to increase with height to an asymptotic value of 1 as r → ∞ (see also Cuntz et al. 1999), but its increase is not necessarily monotonic. We do not explicitly consider the properties of closed magnetic "loops" on the stellar surface, but the radial variation of f(r) takes into account their presence.

Zoom In Zoom Out Reset image size

For each star, we intend to specify f_* on the basis of either direct measurements or empirical scaling relations. The model described in Section 3.1 also requires specifying the value of f at the sharp transition region (TR) between the cool chromosphere and hot corona. We generally know that f_* < f_TR < 1, so in the absence of better information we will apply the assumption that f_TR = f^θ_*, where θ is a dimensionless constant between 0 and 1. For the solar wind models of Cranmer et al. (2007) the exponent θ ranges between about 0.3 and 0.5.

Alfvén waves propagate up from the stellar photosphere, partially reflect back down toward the Sun, develop into strong MHD turbulence, and dissipate gradually (Velli et al. 1991; Matthaeus et al. 1999; Cranmer & van Ballegooijen 2005). Temporarily ignoring the reflection and turbulent cascade, the overall energy balance of an Alfvén wave train is governed by the conservation of wave action. We define the flux of wave action $\tilde{S}$ as

$$\begin{equation} \tilde{S} \, \equiv \, \rho v_{\perp }^{2} V_{\rm A} (1 + M_{\rm A})^{2} A \, = \, \mbox{constant}, \end{equation} \tag{ 13 }$$

where M_A = u/V_A is the Alfvén Mach number and u is the radial outflow speed of the wind (see, e.g., Jacques 1977; Tu & Marsch 1995). Close to the stellar surface, where M_A ≪ 1, this condition is equivalent to energy flux conservation (F_AA = constant). In any case, the constant value of $\tilde{S}$ in Equation (13) is known for each star because the conditions at the photosphere are known (and it is also valid to assume M_A → 0 there as well). The behavior of the wave amplitude as a function of density varies from v_⊥∝ρ^−1/4 close to the star to v_⊥∝ρ^+1/4 at larger distances.

The waves gradually lose energy due to turbulent dissipation, but for locations reasonably close to the stellar surface—e.g., the region shown in Figure 3—it is not a bad approximation to use the undamped form of wave action conservation to compute the radial dependence of v_⊥ (see Cranmer & van Ballegooijen 2005). Wave damping gives rise to plasma heating, and we adopt a phenomenological heating rate that is consistent with the total energy flux that cascades from large to small eddies. This rate is constrained by the properties of the Alfvénic fluctuations at the largest scales, and it does not specify the exact kinetic means of dissipation once the energy reaches the smallest scales. Dimensionally, it is similar to the rate of cascading energy flux derived by von Kármán & Howarth (1938) for isotropic hydrodynamic turbulence. The volumetric heating rate is given by

$$\begin{equation} Q \, = \, \frac{\tilde{\alpha } \rho v_{\perp }^3}{\lambda _{\perp }} \end{equation} \tag{ 14 }$$

(Hollweg 1986; Hossain et al. 1995; Zhou & Matthaeus 1990; Matthaeus et al. 1999; Dmitruk et al. 2002). The dimensionless efficiency factor $\tilde{\alpha }$ depends on the local degree of wave reflection and is discussed further below. The perpendicular length scale λ_⊥ is an effective correlation length for the largest eddies in the turbulent cascade.

MHD turbulence occurs only when there exist counter-propagating Alfvén wave packets along a flux tube. The star naturally creates upward waves, and we assume that linear reflection gives rise to downward waves (Ferraro & Plumpton 1958). We specify the ratio of downward to upward wave amplitudes by the effective reflection coefficient ${\cal R}$, and the efficiency factor $\tilde{\alpha }$ is given by

$$\begin{equation} \tilde{\alpha } \, = \, \alpha _{0} \frac{{\cal R} (1 + {\cal R}) \sqrt{2}}{(1 + {\cal R}^{2})^{3/2}} \end{equation} \tag{ 15 }$$

(see, e.g., Cranmer et al. 2007). At the photospheric lower boundary, we assume total reflection with ${\cal R} = 1$ and thus $\tilde{\alpha } = \alpha _0$. Higher in the stellar atmosphere, we use the low-frequency limiting expression of Cranmer (2010),

$$\begin{equation} {\cal R} \, \approx \, \frac{V_{\rm A} - u_{\infty }}{V_{\rm A} + u_{\infty }}, \end{equation} \tag{ 16 }$$

where the wind's terminal speed is u_∞ and we also assume that M_A ≪ 1 in the atmosphere. This expression also assumes that the wind speed at the point where M_A = 1 (presumably far from the stellar surface) is roughly equal to u_∞. We also set α₀ = 0.5 based on the turbulent transport models of Breech et al. (2009).

If the turbulent heating given by Equation (14) is sufficient to produce a hot (T ≳ 10⁶ K) corona, then the plasma's high gas pressure gradient may provide enough outward acceleration to produce a transition from a subsonic (bound) state near the star to a supersonic (outflowing) state at larger distances (Parker 1958). In this section, we estimate the mass loss rate $\dot{M}$ of such a gas-pressure-driven stellar wind.

We begin by computing Q_* in the photosphere using ρ_* and v_⊥* in Equation (14). For the Sun, we have the observational constraint that λ_⊥* must be about the size of the granular motions that jostle the flux tubes (i.e., roughly 100–1000 km). For other stars we can assume that the horizontal scale of granulation remains proportional to the photospheric pressure scale height (Robinson et al. 2004). Thus, we use

$$\begin{equation} \lambda _{\perp \ast } \, = \, \lambda _{\perp \odot } \, \frac{H_{\ast }}{H_{\odot }}, \end{equation} \tag{ 17 }$$

where H_☉ = 139 km and the models of Cranmer & van Ballegooijen (2005) were used to set the solar normalization of the correlation length to λ_⊥☉ = 300 km.

In the photosphere, we assume the turbulent heating is swamped by radiative gains and losses that are determined by the conditions of local thermodynamic equilibrium (LTE), and the temperature is set by those processes alone. At larger heights in the flux tube, the turbulent heating Q begins to have an effect. We define the chromosphere as the region in which Q is balanced by radiative losses. As one increases in height, however, the density drops to the point where radiative losses alone can no longer balance the imposed heating rate; this occurs at the sharp TR between chromosphere and corona. (See Section 3.2 for cases where this transition does not occur at all.)

In the region between the photosphere and the TR, we assume that the wind flow speed is sufficiently sub-Alfvénic such that v_⊥∝ρ^−1/4. We also assume that λ_⊥ scales with the transverse size of the magnetic flux tube, so that λ_⊥∝A^1/2∝B^−1/2 (Hollweg 1986). Thus, Equation (14) can be rewritten as

$$\begin{equation} \frac{Q_{\rm TR}}{Q_{\ast }} \, = \, \frac{\tilde{\alpha }_{\rm TR}}{\tilde{\alpha }_{\ast }} \left(\frac{\rho _{\rm TR}}{\rho _{\ast }} \right)^{1/4} \left(\frac{B_{\rm TR}}{B_{\ast }} \right)^{1/2}, \end{equation} \tag{ 18 }$$

where $\tilde{\alpha }_{\ast } = 0.5$ and all other photospheric quantities are assumed to be known. We also know that

$$\begin{equation} \frac{B_{\rm TR}}{B_{\ast }} \, = \, \frac{f_{\ast }}{f_{\rm TR}} \, \approx \, f_{\ast }^{1-\theta }, \end{equation} \tag{ 19 }$$

where the last approximation holds if there is a universal relationship between f_* and f_TR as speculated in Section 2.2 above.

Just below the TR, the heating is just barely balanced by radiative cooling. In the optically thin limit, radiative cooling behaves as Q_cool = −n²Λ(T), where n is the number density in the fully ionized TR region. Let us then assume that Q_TR = max |Q_cool|, where

$$\begin{equation} \max | Q_{\rm cool} | \, = \, \frac{\rho _{\rm TR}^{2} \Lambda _{\rm max}}{m_{\rm H}^2}. \end{equation} \tag{ 20 }$$

The quantity Λ_max is the absolute maximum of the radiative loss curve Λ(T), and it occurs roughly at T_TR = 2 × 10⁵ K. The value of Λ_max depends on metallicity. To work out its dependence on Z, we computed a number of radiative loss curves for different metal abundances using version 4.2 of the CHIANTI atomic database (Young et al. 2003) with collisional ionization balance (Mazzotta et al. 1998). We started with a traditional (Grevesse & Sauval 1998) solar abundance mixture (Z/Z_☉ = 1) and then recomputed Λ(T) by varying the metal abundance ratio Z/Z_☉ between 0 and 10. We found that the maxima of the curves were fit well by the following parameterized function,

$$\begin{equation} \frac{\Lambda _{\rm max}}{10^{-23} \, \mbox{erg} \,\, \mbox{cm}^{3} \,\, \mbox{s}^{-1}} \, \approx \, 7.4 + 42 \left(\frac{Z}{Z_{\odot }} \right)^{1.13}. \end{equation} \tag{ 21 }$$

Other examples of the metallicity dependence of Λ(T) have been given by, e.g., Boehringer & Hensler (1989) and Gnat & Sternberg (2007). We have ignored any possible differences between a star's photospheric metal abundances and those in the low corona, although such differences have been measured in some cases (Testa 2010).

With the above assumptions, we solve for the TR density,

$$\begin{equation} \rho _{\rm TR} \, = \, \left[ \frac{\tilde{\alpha }_{\rm TR} Q_{\ast } m_{\rm H}^{2}}{\tilde{\alpha }_{\ast } \rho _{\ast }^{1/4} \Lambda _{\rm max}} \right]^{4/7} {f_{\ast }}^{2(1-\theta)/7} \end{equation} \tag{ 22 }$$

and we also derive the heating rate at the TR to be

$$\begin{equation} Q_{\rm TR} \, = \, \left(\frac{\tilde{\alpha }_{\rm TR} Q_{\ast }}{\tilde{\alpha }_{\ast }} \right)^{8/7} \left(\frac{m_{\rm H}^{2}}{\rho _{\ast }^{2} \Lambda _{\rm max}} \right)^{1/7} {f_{\ast }}^{4(1-\theta)/7}. \end{equation} \tag{ 23 }$$

A potential roadblock to solving Equations (22)–(23) is that we do not initially know the value of $\tilde{\alpha }_{\rm TR}$. This quantity depends on the reflection coefficient ${\cal R}$, which depends on the Alfvén speed V_A at the TR (see Equation (16)), which in turn depends on the unknown value of ρ_TR. In practice, we solve these equations iteratively. We start with an initial estimate of ${\cal R} = 0.5$, we compute $\tilde{\alpha }_{\rm TR}$, ρ_TR, and V_A at the TR, and then we recompute ${\cal R}$ for the next iteration. In all cases the process converges to a self-consistent set of values (with a relative accuracy of ∼10⁻⁷) in no more than 20 iterations.

The mass loss rate of the stellar wind is determined by the heating rate Q_TR as well as other sources and sinks of energy at the TR. The general idea that the solar wind's mass flux is set by the energy balance at the TR was first discussed by Hammer (1982). Hansteen & Leer (1995) worked out the basic scaling argument that is used below (see also Leer et al. 1982; Withbroe 1988; Schwadron & McComas 2003). In the low corona and wind, the time-steady equation of internal energy conservation is

$$\begin{equation} \frac{1}{A} \frac{\partial }{\partial r} \left\lbrace A \left[ F_{\rm H} - F_{\rm cond} + \rho u \left(\frac{u^2}{2} - \frac{{\it GM}_{\ast }}{r} \right) \right] \right\rbrace \, = \, 0, \end{equation} \tag{ 24 }$$

where F_H is the energy flux associated with the heating, F_cond is the energy flux transported by heat conduction along the field, and u is the outflow speed. The term in braces is constant as a function of radius, so it is straightforward to equate its value at the TR to its asymptotic value at r → ∞. The kinetic energy term proportional to u² is assumed to be negligibly small at the TR, but we assume it dominates the energy balance at large distances. Thus,

$$\begin{equation} A_{\rm TR} (F_{\rm H, TR} - F_{\rm cond}) - (\rho u A)_{\rm TR} \frac{{\it GM}_{\ast }}{R_{\ast }} \, = \, (\rho u A)_{\infty } \frac{u_{\infty }^2}{2}, \end{equation} \tag{ 25 }$$

where F_{H, TR} is the heat flux F_H at the TR, and we realize that the product ρuA is also constant via mass flux conservation. We also make the key assumption that u_∞ = V_esc = (2GM_*/R_*)^1/2, and thus we can write

$$\begin{equation} \dot{M} \, \equiv \, \rho u A \, = \, \frac{4\pi R_{\ast }^{2} f_{\rm TR}}{V_{\rm esc}^2} (F_{\rm H, TR} - F_{\rm cond}). \end{equation} \tag{ 26 }$$

To evaluate Equation (26) we need to estimate the value of F_{H, TR}. Formally, Q = |∇ · F_H|, so to determine the magnitude F_{H, TR} one would have to integrate Q(r) along the flux tube. Taking account of the expanding flux tube area A∝B⁻¹, and also assuming that r_TR ≈ R_*,

$$\begin{equation} F_{\rm H, TR} \, = \, \frac{1}{A_{\rm TR}} \int _{R_{\ast }}^{\infty } dr \, Q(r) \, A(r). \end{equation} \tag{ 27 }$$

Specifically, if Q∝r^−β and A∝r^γ, then

$$\begin{equation} F_{\rm H, TR} \, = \, \frac{Q_{\rm TR} R_{\ast }}{|\beta - \gamma - 1|} \, \equiv \, Q_{\rm TR} R_{\ast } h. \end{equation} \tag{ 28 }$$

Rather than specifying β and γ, we estimate the dimensionless scaling factor h by extracting both Q_TR and F_{H, TR} from the self-consistent solar wind models of Cranmer et al. (2007). For a range of fast and slow solar wind solutions, we found that Q_TR is typically between 1.5 × 10⁻⁵ and 4 × 10⁻⁵ erg cm⁻³ s⁻¹, and F_{H, TR} is typically between 8 × 10⁵ and 3 × 10⁶ erg cm⁻² s⁻¹. This results in h usually being between 0.5 and 1.5.

It is important to also verify that F_{H, TR} is less than the energy flux carried "passively" by the Alfvén waves as they propagate up from the photosphere. The latter quantity, which we call F_{A, TR}, represents the upper limit of available energy in the waves (at the TR) that can be extracted by the turbulent heating. Assuming that wave flux is conserved (i.e., that M_A ≪ 1 at the TR), then F_{A, TR} = f_*F_A*/f_TR. For the cool-star models discussed in Section 5, we found that the ratio F_{H, TR}/F_{A, TR} is usually around 0.1–0.5. Only in two cases did it exceed 1 (albeit with values no larger than 1.5), and in those cases we capped F_{H, TR} to be equal to F_{A, TR} to maintain energy conservation.

To evaluate Equation (26), we also need to estimate the magnitude of the downward conductive flux F_cond from the hot corona. For the solar TR and low corona, Withbroe (1988) found there to be an approximate balance between conduction and radiation losses. Withbroe (1988) determined that F_cond ≈ c_radP_TR, where P_TR is the gas pressure at the TR, and the constant of proportionality is

$$\begin{equation} c_{\rm rad} \, = \, \sqrt{\frac{\kappa _e}{2 k_{\rm B}^2} \int _{T_0}^{T_{\rm TR}} \Lambda (T) \, T^{1/2} \, dT}, \end{equation} \tag{ 29 }$$

where κ_e is the electron thermal conductivity, T₀ ≈ 10⁴ K is a representative chromospheric temperature, T_TR = 2 × 10⁵ K, and c_rad has units of speed. We evaluated the above integral to be able to scale out the metallicity-dependent factor given in Equation (21) above and found that

$$\begin{equation} c_{\rm rad} \, \approx \, 14 \, \sqrt{\frac{\Lambda _{\rm max}(Z)}{\Lambda _{\rm max}(Z_{\odot })}} \,\,\,\, \mbox{km} \,\, \mbox{s}^{-1}. \end{equation} \tag{ 30 }$$

We used this expression to estimate F_cond = c_radP_TR. For the specific case of the Sun, conduction is relatively unimportant in open flux tubes, since F_cond ≲ 0.05F_{H, TR}. For the other stars modeled in this paper, the ratio F_cond/F_{H, TR} spanned several orders of magnitude from 10⁻⁴ to 10⁻¹. In the eventuality that the estimated value of F_cond may exceed F_{H, TR}, one would need an improved description of the coronal temperature T(r) to compute a more accurate value of the conduction flux. In our numerical code, however, we do not allow F_cond to exceed a value of ξF_{H, TR}, where ξ is an arbitrary constant that we fixed to a value of 0.9. This condition was not met for any of the stars in the database of Section 5.

Once these energy fluxes are computed, we then compute $\dot{M}$ using Equations (23), (26), and (28), as well as the other definitions given above. Interestingly, this can be done without needing to know the temperature profile T(r). From a certain perspective, the corona's thermal response to the heating rate Q may be considered to be just an intermediate step toward the "final" outcome of a kinetic-energy-dominated outflow far from the star.⁵ However, it should be possible to estimate the maximum coronal temperature T_max by inverting scaling laws given by, e.g., Hammer et al. (1996) and Schwadron & McComas (2003).

The mass loss rate given by Equation (26) depends on our assumption that u_∞ = V_esc. Equation (25) shows that larger assumed values of u_∞ would give rise to lower mass loss rates, and smaller values of u_∞ would give larger mass loss rates. For the solar wind there is roughly a factor of three variation in u_∞, from about 0.4V_esc to 1.3V_esc. For other stars, it is rare to see observations where u_∞ exceeds V_esc, and Judge (1992) found generally that u_∞ < V_esc for luminous evolved stars. Even in the extreme case of u_∞ = 0, however, Equation (25) would give only two times the mass loss assumed by Equation (26). When compared to the larger typical observational uncertainties in $\dot{M}$, factors of two are not a major concern.

The Sun's mass loss rate of 2 × 10⁻¹⁴ to 3 × 10⁻¹⁴ M_☉ yr⁻¹ is modeled reasonably well with the model described here. The photospheric energy flux of Alfvén waves is F_A* ≈ 1.5 × 10⁸ erg cm⁻² s⁻¹, and the photospheric wave amplitude is v_⊥* ≈ 0.28 km s⁻¹ (Cranmer et al. 2007). Although magnetogram observations sometimes give filling factors f_* as large as 1% at solar maximum, values of 0.1% tend to better represent the coronal holes that are connected to the largest volume of open flux tubes (see Figure 3 of Cranmer & van Ballegooijen 2005). Assuming f_* = 0.001 and values for the other constants of α₀ = 0.5 and θ = 1/3, Equations (22)–(23) give ρ_TR ≈ 5 × 10⁻¹⁶ g cm⁻³ and Q_TR ≈ 4 × 10⁻⁵ erg cm⁻³ s⁻¹, which are in agreement with the models of Cranmer et al. (2007) and others. Thus, with h = 0.5, Equation (26) gives $\dot{M} \approx 3.5 \times 10^{-14}$ M_☉ yr⁻¹.

Although the above calculation of $\dot{M}$ is relatively straightforward, it has not been boiled down to a simple scaling law such as that of Reimers (1975, 1977), Mullan (1978), or Schröder & Cuntz (2005). However, if we make the further assumptions that $\tilde{\alpha }$ and h are fixed constants, and that F_{H, TR} ≫ F_cond, we can isolate several interesting scalings.

• 1.  
  The ultimate driving of the wind comes from the basal flux of Alfvén wave energy F_A*. Schröder & Cuntz (2005) assumed that $\dot{M}$ scales linearly with F_A*, but in our case we can combine the above equations with the definition of Q_* to find $\dot{M} \propto F_{\rm A \ast }^{12/7}$, which is noticeably steeper than a pure linear dependence. This positive feedback is qualitatively similar to what occurs in radiatively driven winds of more massive stars, for which the mass loss rate is proportional to the radiative flux (or luminosity) to a power larger than one (i.e., $\dot{M} \propto L_{\ast }^{1.7}$; Castor et al. 1975; Owocki 2004).
• 2.  
  Extracting the dependence on magnetic filling factor, we found that $\dot{M} \propto f_{\ast }^{(4 + 3\theta)/7}$. Using the range of θ from solar models (0.3–0.5), this gives a relatively narrow range of exponents, $\dot{M} \propto f_{\ast }^{0.7}$ to f^0.8_*. Saar (1996a) estimated that f_*∝P^−1.8_rot for rotation periods P_rot > 3 days (see also Section 4). Thus, for stars in the unsaturated part of the age-activity-rotation relationship, it may be possible to estimate $\dot{M} \propto P_{\rm rot}^{-1.3}$.

These simple scaling relations are given only for illustrative purposes (see also Equation (45) below). The predictions of our "hot" coronal mass loss model should be considered to be the solutions of the full set of Equations (17)–(30).

In a high-density stellar atmosphere, it is possible that the turbulent heating described by Equation (14) could be balanced by radiative cooling even very far from the star. In that case, hot coronal temperatures may never occur (see, e.g., Suzuki 2007; Cranmer 2008). The density dependence of the heating rate Q determines whether radiative cooling remains important at large distances, and Figure 4 shows two examples of how Q may vary as a function of ρ. Near the stellar surface, where v_⊥∝ρ^−1/4 and B∝ρ^1/2, we can combine various assumptions to estimate that Q∝ρ^1/2. However, further from the star, where v_⊥∝ρ^+1/4 and B∝ρ, the density dependence becomes steeper, with Q∝ρ^9/4.

Zoom In Zoom Out Reset image size

Figure 4 compares the modeled heating rates with the "maximum cooling boundary" implied by Equation (20). The solar model crosses the boundary, and thus undergoes a transition to a hot corona. The solid curve was taken from a numerical model of fast solar wind from a polar coronal hole (Cranmer et al. 2007). On the other hand, a model for a representative late-type giant sits to the right of the cooling boundary, which implies that radiative losses can maintain the circumstellar temperature at chromospheric values of ∼10⁴ K even at large distances. The two models differ because the density at which M_A ≈ 1 (where the Q curves undergo a change in slope) for the giant is several orders of magnitude larger than the corresponding density for the solar case. This density is an output of a given mass loss model and cannot be specified a priori.

In this section, we develop a model for cool-star mass loss under the assumption of strong radiative cooling. In this case the Parker (1958) gas pressure driving mechanism cannot drive a significant outflow. However, when the flux of Alfvén waves is large, they can impart a strong bulk acceleration to the plasma due to wave pressure, which is a nondissipative net ponderomotive force exerted by virtue of wave propagation through an inhomogeneous medium (Bretherton & Garrett 1968; Jacques 1977). The subsequent calculation of $\dot{M}$ for a "cold wave-driven" stellar wind largely follows the development of Holzer et al. (1983) (see also Cranmer 2009).

Three key assumptions are (1) that the Alfvén wave amplitudes in the wind are larger than the local sound speeds, (2) that there is negligible wave damping between the stellar surface and the wave-modified critical point of the flow, and (3) that the critical point occurs far enough from the stellar surface that f ≈ 1 there (i.e., the flux tube expansion becomes radial). A fourth assumption from Holzer et al. (1983)—which was initially not applied here but later found to be valid—is that the stellar wind is sub-Alfvénic at the critical point (i.e., that M_A ≪ 1 at the critical point). Cranmer (2009) constructed a set of numerical models for the cold polar outflows of T Tauri stars that did not make the fourth assumption and found that for a wide range of parameters the assumption was justified. Thus, here we use the value of the critical radius given by Equation (35) of Holzer et al. (1983), in which all four of the above assumptions were applied, and

$$\begin{equation} \frac{r_{\rm crit}}{R_{\ast }} \, \approx \, \frac{7/4}{1 + (v_{\perp \ast } / V_{\rm esc})^2}. \end{equation} \tag{ 31 }$$

Once the critical point radius is known, it becomes possible to use the known properties of the Alfvén waves to determine the wind velocity and density at the critical point. Holzer et al. (1983) found analytic solutions for these quantities in the limiting case of M_A ≪ 1 at the critical point. Here, we describe a slightly more self-consistent way of computing M_A and the mass loss rate, but we also continue to use Equation (31) that was derived in the limit of M_A ≪ 1. There are three unknown quantities and three equations to constrain them. The three unknowns are the critical point values of the wind speed u, density ρ, and wave amplitude v_⊥. The first equation is the constraint that the right-hand side of the time-steady momentum equation must sum to zero at the critical point of the flow (e.g., Parker 1958). For the conditions described above, this gives

$$\begin{equation} \frac{2 u_{\rm crit}^{2}}{r_{\rm crit}} - \frac{{\it GM}_{\ast }}{r_{\rm crit}^2} \, = \, 0, \end{equation} \tag{ 32 }$$

and it is solved straightforwardly for u_crit. The second and third equations are, respectively, the definition of the critical point velocity in the "cold" limit of zero gas pressure,

$$\begin{equation} u_{\rm crit}^{2} \, = \, \frac{v_{\perp }^2}{4} \left(\frac{1 + 3 M_{\rm A}}{1 + M_{\rm A}} \right) \end{equation} \tag{ 33 }$$

and the conservation of wave action as given by Equation (13). The fact that V_A appears in these equations and depends on the (still unknown) density makes it difficult to find an explicit analytic solution for ρ_crit. We again use iteration from an initial guess to reach a self-consistent solution for u, ρ, and v_⊥ at the critical point. The stellar wind's mass loss rate is thus determined from $\dot{M} = 4\pi r_{\rm crit}^{2} u_{\rm crit} \rho _{\rm crit}$.

Because the mass loss rate is set at the critical point, we do not need to specify the terminal speed u_∞. For most implementations of the above model, the denominator in Equation (31) is close to unity and thus we have u²_crit ≈ V_esc²/7, or that u_crit is about 38% of the presumed value of u_∞. Of course, there have been stellar wind models with nonmonotonic radial variations of u(r), with u_crit > u_∞ (e.g., Falceta-Gonçalves et al. 2006). It is also possible for "too much" mass to be driven past the critical point, such that parcels of gas may be decelerated to stagnation at some height above r_crit and thus would want to fall back down toward the star. In reality, this parcel would collide with other parcels that are still accelerating, and a stochastic collection of shocked clumps is likely to result. Interactions between these parcels may result in an extra degree of collisional heating that could act as an extended source of gas pressure to help maintain a mean net outward flow. Situations similar to this have been suggested to occur in the outflows of pulsating cool stars (Bowen 1988; Struck et al. 2004), T Tauri stars (Cranmer 2008), and luminous blue variables (van Marle et al. 2009).

A proper treatment of a stellar wind powered by MHD turbulence—and accelerated by a combination of gas pressure and wave pressure effects—requires a self-consistent numerical solution to the conservation equations (e.g., Cranmer et al. 2007; Suzuki 2007; Cohen et al. 2009; Airapetian et al. 2010). However, in this paper, we explore simpler ways of estimating the combined effects of both processes.

Sections 3.1 and 3.2 gave us independent estimates for the mass loss rate assuming only gas pressure or wave pressure was active in the flux tube of interest. We refer to these two mass loss rates as $\dot{M}_{\rm hot}$ and $\dot{M}_{\rm cold}$, respectively. It seems clear that when one of these values is much larger than the other, then one process is dominant and the actual mass loss rate should be close to that larger value. For the manifestly "hot" example of the Sun, we found that $\dot{M}_{\rm hot} / \dot{M}_{\rm cold} \approx 20$, which correctly implies that gas pressure driving is dominant. For most examples of late-type giants with L_* > 10 L_☉, the ratio $\dot{M}_{\rm hot} / \dot{M}_{\rm cold}$ was found to decrease to values between about 0.1 and 3. This could mean that gas and wave pressure gradients are of the same order of magnitude for these stars.

One of the most straightforward things that can be done is to assume the combined effect of gas and wave pressure produces a mass loss rate equal to the sum of the two individual components, $\dot{M}_{\rm hot} + \dot{M}_{\rm cold}$. This preserves the idea that one dominant mechanism should determine $\dot{M}$ when the other would predict a negligibly small effect. It also makes sense based on Equation (24), which shows how the energy fluxes sum together linearly in the internal energy equation. If there were multiple sources of input energy flux, Equation (26) would show that the resulting mass loss rate should be proportional to their sum.

However, there is one complication that hinders us from simply adding together $\dot{M}_{\rm hot}$ and $\dot{M}_{\rm cold}$. The calculation of $\dot{M}_{\rm hot}$ from Section 3.1 contains the assumption that the TR turbulent heating always obeys the near-star density scaling Q∝ρ^1/2. It therefore predicts that all stars eventually undergo a transition to a hot corona. For some stars (like the late-type giant in Figure 4), however, we know that there should be no corona and it is erroneous to assume that $\dot{M}_{\rm hot}$ has any real meaning. Thus, for each model we compute the wind speed at the TR from mass flux conservation,

$$\begin{equation} u_{\rm TR} \, = \, \frac{\dot{M}_{\rm hot}}{4\pi R_{\ast }^{2} f_{\rm TR} \rho _{\rm TR}}, \end{equation} \tag{ 34 }$$

and we demand that for $\dot{M}_{\rm hot}$ to have a consistent interpretation, the TR Mach number M_{A, TR} = u_TR/V_{A, TR} should be much smaller than one. As expected, this condition was found to be violated for late-type giants having L_* ≳ 100 L_☉. Thus, in these cases we should replace $\dot{M}_{\rm hot}$ with either a drastically reduced value or zero—the latter in cases where the Q(ρ) curve always falls to the right of the maximum cooling boundary in Figure 4. After some experimentation, we found that reducing the initially computed value of Q_TR by a factor of exp (− 4M²_{A, TR}) does a reasonably good job of reproducing the result of using a more consistent Q(ρ) function. Thus, we propose that the summing of the "hot" and "cold" mass loss rates be done with the following approximate expression:

$$\begin{equation} \dot{M} \, \approx \, \dot{M}_{\rm cold} + \dot{M}_{\rm hot} \exp \big(- 4 M_{\rm A, TR}^{2} \big), \end{equation} \tag{ 35 }$$

where $\dot{M}_{\rm hot}$ and M_{A, TR} are computed using the assumptions of Section 3.1 and $\dot{M}_{\rm cold}$ is computed using the model given in Section 3.2.

Figure 8 is a broad overview of observed stellar mass loss. It plots the locations of individual stars in a Hertzsprung–Russell type diagram with their mass loss rates shown as symbol color (see also de Jager et al. 1988). A box illustrates the approximate regime of parameter space covered by the models developed in this paper; it extends from the main sequence up through the regime of the so-called hybrid chromosphere stars (Hartmann et al. 1980), and possibly also into the parameter space of cool luminous supergiants. In addition to the cool-star data discussed below, we also include in Figure 8 measured mass loss rates of hot, massive stars (Waters et al. 1987; Lamers et al. 1999; Mokiem et al. 2007; Searle et al. 2008), FGK supergiants (de Jager et al. 1988), AGB stars (Bergeat & Chevallier 2005; Guandalini 2010), red giants in globular clusters (Mészáros et al. 2009), and M dwarfs in pre-cataclysmic variable binaries (Debes 2006). Many of these stars are not included in the subsequent analysis because we have no firm rotation periods or magnetic activity indices for them.

Zoom In Zoom Out Reset image size

Table 2 lists the properties of 47 stars for which our knowledge appears to be complete enough to be able to compare theoretical and observed values of $\dot{M}$. For the Sun, the range of volume-integrated mass loss rates comes from Wang (1998). The sources for all listed values are given as numbered references that continue the sequence started in Table 1; the citations corresponding to numbers 1–32 are given in Table 1. In cases where T_eff, log  g, or [Fe/H] were estimated from the PASTEL database (Soubiran et al. 2010), we averaged together multiple measurements when more than one was given. In the few cases where the same star appears in both Table 1 and Table 2, for consistency's sake we will recompute f_* from the star's rotation period when calculating theoretical mass loss rates (see Section 5.2).

Table 2. Cool Star Mass Loss Data

Name	Teff (K)	log  g	M*/M☉	R*/R☉	L*/L☉	Prot (d)	[Fe/H]	$-\log (\dot{M} / [M_{\odot }\, \mbox{yr}^{-1}])$	Ref.
Sun	5770	4.44	1	1	1	25.3	0	13.5–13.7	⋅⋅⋅
α Cen A (G2 V)	5886	4.31	1.105	1.224	1.622	29	+0.197	13.40 (A+B)	33, 34, 23, 7
α Cen B (K0 V)	5473	4.54	0.934	0.863	0.603	36.2	+0.230	13.40 (A+B)	33, 34, 35, 7
70 Oph A (K0 V)	5300	4.52	0.89	0.86	0.53	19.7	+0.040	11.70 (A+B)	14, 34, 23
70 Oph B (K5 V)	4390	4.65	0.73	0.67	0.15	34	+0.040	11.70 (A+B)	14, 34, 36
Eri (K4.5 V)	5094	4.60	0.83	0.754	0.345	11.7	−0.097	12.22	18, 34, 23, 7
61 Cyg A (K5 V)	4425	4.63	0.69	0.665	0.153	35.4	−0.193	14.00	22, 34, 23, 7
Ind (K5 V)	4635	4.54	0.70	0.745	0.231	22	−0.088	14.00	37, 34, 23, 7
36 Oph A (K5 V)	5135	4.54	0.602	0.69	0.299	20.3	−0.206	12.52 (A+B)	19, 34, 7
36 Oph B (K5 V)	5103	4.58	0.486	0.59	0.213	22.9	−0.195	12.52 (A+B)	19, 34, 7
ξ Boo A (G8 V)	5551	4.57	0.86	0.801	0.550	6.2	−0.122	13.00 (A+B)	12, 34, 38, 7
ξ Boo B (K4 V)	4350	4.80	0.70	0.550	0.0977	11.5	−0.122	13.00 (A+B)	12, 34, 38, 7
61 Vir (G5 V)	5560	4.39	0.946	0.972	0.804	29	−0.002	14.22	39, 34, 36, 7
δ Eri (K0 IV)	5025	3.75	1.122	2.33	3.185	55.3	+0.069	13.10	40, 34, 23, 7
α Boo (K1.5 III)	4290	1.76	1.10	23	170	447	−0.53	9.60	41, 42, 7
α Tau (K5 III)	3898	1.33	1.5	44	394	648	−0.180	10.83	43, 44
γ Dra (K5 III)	3985	1.53	3.0	49	535	557	−0.150	11.06	43, 44
HR 6902 (G9 IIb)	4900	1.99	3.86	33	566	220	+0.430	10.68	45, 46, 47
β And (M0 III)	3742	1.55	1.5	34	204	188	−0.04	10.19	48, 49, 7
β UMi (K4 III)	4040	1.27	1.3	44	475	1030	−0.26	10.01	48, 50, 51, 7
μ UMa (M0 III)	3700	0.69	1.5	92	1430	488	0.00	8.92	48, 49, 7
α TrA (K2 II-III)	4150	1.50	23.7	143	5500	888	−0.06	9.77	52, 53, 7
λ Vel (K4 Ib-II)	3820	0.64	7.0	210	8510	1250	+0.23	8.52	54, 55, 7
BD +01 3070 (RGB)	5130	2.70	0.749	6.4	25.6	50.9	−1.85	8.76	56, 57
BD +05 3098 (RGB)	4930	2.00	0.746	14.3	109	109	−2.40	8.91	56, 57
BD +09 2574 (RGB)	4860	2.10	0.753	12.8	82.5	204	−1.95	8.94	56, 57
BD +09 2870 (RGB)	4600	1.40	0.864	30.7	381	235	−2.37	8.43	56, 57
BD +10 2495 (RGB)	4920	2.12	0.723	12.3	80.0	156	−1.83	9.06	56, 57
BD +12 2547 (AGB)	4610	1.50	0.780	26.0	275	265	−0.72	8.15	56, 57
BD +17 3248 (RHB)	5250	2.21	0.458	8.80	53.1	64.8	−2.02	8.95	56, 57
BD +18 2757 (AGB)	4840	1.43	0.446	21.3	225	127	−2.19	8.24	56, 57
BD +18 2976 (RGB)	4550	1.30	0.769	32.5	408	248	−2.40	7.65	56, 57
BD −03 5215 (RHB)	5420	2.60	0.884	7.80	47.4	42.5	−1.66	8.68	56, 57
HD 083212 (RGB)	4550	1.40	0.663	26.9	280	146	−1.49	8.14	56, 57
HD 101063 (SGB)	5070	3.40	1.19	3.60	7.73	28.6	−1.13	9.56	56, 57
HD 107752 (AGB)	4750	1.70	0.838	21.4	210	185	−2.88	8.24	56, 57
HD 110885 (RHB)	5330	2.50	0.757	8.10	47.8	39.3	−1.44	8.66	56, 57
HD 111721 (RGB)	5080	2.35	0.460	7.50	33.8	74.5	−1.26	9.20	56, 57
HD 115444 (RGB)	4750	1.62	0.584	19.6	176	169	−2.77	8.56	56, 57
HD 119516 (RHB)	5440	2.37	0.626	8.60	58.4	39.8	−2.50	8.38	56, 57
HD 121135 (AGB)	4925	1.90	0.789	16.5	144	76.3	−1.57	8.23	56, 57
HD 122956 (RGB)	4600	1.35	0.447	23.4	221	130	−1.78	8.03	56, 57
HD 135148 (RGB)	4275	0.80	0.239	32.2	312	164	−1.90	7.83	56, 57
HD 195636 (RHB)	5370	2.17	0.442	9.10	62.1	16.9	−2.83	8.05	56, 57
EV Lac (M3.5 V)	3168	4.80	0.315	0.369	0.0124	4.38	−0.200	13.70	26, 27, 34 23
V Hya (N6, AGB)	2160	−0.89	4.20	945	17540	576	+0.10	5.12	58, 59, 60
89 Her (F2 Ib)	6550	0.60	0.61	64.8	6970	143	−0.41	8.00	61, 62, 63, 7

References. (1)–(32) See Table 1; (33) Porto de Mello et al. 2008; (34) Wood et al. 2005a; (35) DeWarf et al. 2010; (36) Baliunas et al. 1996; (37) Janson et al. 2009; (38) Wood & Linsky 2010; (39) Vogt et al. 2010; (40) Hekker & Aerts 2010; (41) Schröder & Cuntz 2007; (42) Carney et al. 2008; (43) Robinson et al. 1998; (44) Cayrel de Strobel et al. 2001; (45) Kirsch et al. 2001; (46) Griffin 1988; (47) Marshall 1996; (48) Judge & Stencel 1991; (49) Massarotti et al. 2008; (50) Tarrant et al. 2008; (51) de Medeiros & Mayor 1999; (52) Ayres et al. 2007; (53) Harper et al. 1995; (54) Carpenter et al. 1999; (55) Setiawan et al. 2004; (56) Dupree et al. 2009; (57) Cortés et al. 2009; (58) Bergeat & Chevallier 2005; (59) Lambert et al. 1986; (60) Barnbaum et al. 1995; (61) Sargent & Osmer 1969; (62) Danziger & Faber 1972; (63) Stasińska et al. 2006.

Download table as:  ASCIITypeset image

For binary systems with astrospheric measurements of $\dot{M}$ (see, e.g., Wood et al. 2002), the numbers given are assumed to be the sum of both stars' mass loss rates. We list that same value for both components and denote it with "(A+B)." We did not utilize the published astrospheric measurements of Proxima Cen and 40 Eri A, which gave only upper limits on $\dot{M}$, and λ And and DK UMa, which had uncertain detections of astrospheric H i Lyα absorption (Wood et al. 2005a, 2005b).

At the bottom of Table 2 we list three stars that have parameters at the outer bounds of what we intend to model. They are test cases for the limits of applicability of the physical processes summarized in Section 3. EV Lac is an active M dwarf and flare star that probably has a fully convective interior (e.g., Osten et al. 2010). Such stars may exhibit qualitatively different mechanisms of mass loss and rotation-activity correlation than do stars higher up the main sequence (Mullan 1996; Reiners & Basri 2007; Irwin et al. 2011; Martínez-Arnáiz et al. 2011; Vidotto et al. 2011). V Hya is an N-type carbon star with an extended and asymmetric AGB envelope and evidence for rapid rotation (Barnbaum et al. 1995; Knapp et al. 1999). 89 Her is a post-AGB yellow supergiant with multiple detections of circumstellar nebular material (Sargent & Osmer 1969; Bujarrabal et al. 2007). It is worthwhile to investigate to what extent the mass loss mechanisms proposed in this paper could be applicable to these kinds of stars.

Not all stars in Table 2 have precise measurements for their rotation period. For 61 Vir and 70 Oph B, we used published estimates of the rotation period that were obtained from the known correlation between rotation and chromospheric Ca ii activity (Baliunas et al. 1996). For essentially all stars more luminous than ∼5 L_☉ (with the exception of HR 6902; see Griffin 1988) we estimated P_rot via spectroscopic determinations of vsin i from the rotational broadening of photospheric absorption lines. The inclination angle i is the main unknown quantity. Chandrasekhar & Münch (1950) found that for an isotropically distributed set of inclination vectors, the mean value of sin i is π/4. Thus, we estimate a mean rotation period

$$\begin{equation} \langle P_{\rm rot} \rangle \, = \, \frac{2\pi R_{\ast }}{(4 / \pi) \, v \sin i}. \end{equation} \tag{ 40 }$$

Note that the median of sin i for an isotropic distribution is not equal to the mean; the former is given by $\sqrt{3}/2$. In order to encompass both values, as well as the majority of "most likely" values of P_rot, we can adopt generous uncertainty limits for which we will estimate f_* and the other derived quantities for mass loss. For the isotropic distribution of direction vectors, the quantity sin i falls between 0.5 and 1 approximately 87% of the time. This is a reasonably good definition for uncertainty bounds that would correspond to ±1.5σ if the distribution were Gaussian. Thus, for stars with only vsin i measurements, we use the following values as error bars on the derived rotation period:

$$\begin{equation} \frac{2}{\pi } \lesssim \frac{P_{\rm rot}}{\langle P_{\rm rot} \rangle } \lesssim \frac{4}{\pi }. \end{equation} \tag{ 41 }$$

We applied the combined model for mass loss that culminated in Equation (35) to the stars listed in Table 2. Below we show results of direct forward modeling, i.e., utilizing a known relationship for f_* as a function of Rossby number. First, however, we wanted to investigate whether or not a single monotonic relationship for f_*(Ro) could produce mass loss rates that were even remotely close to the measured values. Thus, we produced trial grids of models in which f_* was treated as a free parameter. For each star, we varied f_* from 10⁻⁵ to 1 and found the empirical value of the filling factor (f_emp) for which the modeled value of $\dot{M}$ matched the observed value given in Table 2. For the four binaries that have only systemic measurements of $\dot{M}$ we summed the model predictions for each component and made a single comparison with the observations.

Figure 9 shows the result of this process of "working backward" from the measured mass loss rates. The empirically constrained f_emp values are plotted against Rossby number, which is defined with (1) the simple Gunn et al. (1998) ZAMS value for τ_c (Equation (36)) and (2) a gravity-modified version of τ_c that gives giants larger convective overturn times.⁷ We varied the exponent in the gravity modification term and found that multiplying the ZAMS τ_c by (g_☉/g)^0.18 gives the narrowest distribution of f_* versus Ro. The optimal exponent of 0.18 is very close to the value of 0.23 that we found reproduced the results of Hall (1994) and Gondoin (2005, 2007). In Figure 9 we also show the same curves from Figure 7(b) that outline the measured range of filling factors. The lower envelope curve f_min, defined in Equation (38), appears to be a good match to the gravity-modified empirical values f_emp. This provides circumstantial evidence that f_min is indeed an appropriate proxy for the filling factor of open flux tubes as a function of Rossby number.

Zoom In Zoom Out Reset image size

We now put aside the empirical estimates for the filling factor and use only Equation (38) for f_* = f_min in the remainder of this paper. Table 3 shows some of the predicted properties of stellar coronae and winds for the 47 stars in our database. There were only seven stars for which Equation (38) gave a filling factor below the adopted "floor" value of f_basal = 10⁻⁴; we replaced f_min with f_basal in those cases. Table 3 also gives F_{H, TR}, the coronal heat flux deposited at the TR for each star. It may be useful to use this to predict the X-ray flux associated with open-field regions on these stars, but we should note that the closed-field regions (which we do not model) are likely to dominate the observed X-ray emission. We also list the various components of Equation (35) so that the contributions of gas pressure and wave pressure can be assessed (see below).

Table 3. Theoretical Wind Properties of Cool Stars

Name	Ro	log fmin	B* (G)	log FH, TR	$\dot{M}_{\rm hot} / \dot{M}_{\rm cold}$	log MA, TR	$-\log (\dot{M} / [M_{\odot } \,\mbox{yr}^{-1}])$
Sun	1.960	−2.996	1513.05	6.14	19.73	−2.42	13.44
α Cen A (G2 V)	2.074	−3.078	1308.79	6.17	11.43	−2.17	13.21
α Cen B (K0 V)	1.755	−2.835	1545.97	5.62	51.73	−2.93	14.09
70 Oph A (K0 V)	0.996	−2.025	1666.45	5.99	194.1	−3.25	13.43
70 Oph B (K5 V)	1.027	−2.067	2130.57	4.55	417.3	−4.50	15.17
Eri (K4.5 V)	0.519	−1.174	1832.80	5.97	977.4	−3.94	13.31
61 Cyg A (K5 V)	1.107	−2.173	2145.92	4.59	246.1	−4.43	15.15
Ind (K5 V)	0.755	−1.646	1941.88	5.15	551.7	−4.22	14.25
36 Oph A (K5 V)	0.923	−1.918	1927.79	5.67	228.1	−3.63	13.84
36 Oph B (K5 V)	1.021	−2.059	1973.54	5.51	173.1	−3.66	14.16
ξ Boo A (G8 V)	0.381	−0.853	1788.75	6.69	1335	−3.59	12.42
ξ Boo B (K4 V)	0.340	−0.755	2620.92	4.83	5964	−5.37	14.68
61 Vir (G5 V)	1.765	−2.843	1524.08	5.99	27.96	−2.63	13.56
δ Eri (K0 IV)	1.844	−2.907	1077.96	5.77	12.04	−2.18	12.72
α Boo (K1.5 III)	4.217	−4.000	411.54	5.23	0.57	−0.31	10.33
α Tau (K5 III)	4.183	−4.000	270.35	5.17	0.75	0.11	9.88
γ Dra (K5 III)	4.097	−4.000	301.38	5.24	0.85	−0.12	9.99
HR 6902 (G9 IIb)	3.157	−3.692	287.82	6.11	1.28	0.13	9.86
β And (M0 III)	1.229	−2.321	311.31	5.61	0.47	−0.85	8.86
β UMi (K4 III)	6.960	−4.000	262.45	5.32	0.85	0.28	9.72
μ UMa (M0 III)	2.182	−3.152	165.83	5.81	0.63	0.68	7.93
α TrA (K2 II-III)	7.010	−4.000	270.87	5.55	1.25	−0.13	9.32
λ Vel (K4 Ib-II)	5.808	−4.000	136.66	5.67	2.15	1.05	8.47
BD +01 3070 (RGB)	1.122	−2.192	792.36	6.49	2.36	−1.48	10.14
BD +05 3098 (RGB)	1.600	−2.700	553.26	6.37	0.50	−0.63	9.08
BD +09 2574 (RGB)	3.001	−3.618	632.02	5.74	0.49	−0.62	10.17
BD +09 2870 (RGB)	2.249	−3.196	476.21	6.00	0.37	−0.17	8.68
BD +10 2495 (RGB)	2.390	−3.284	602.67	6.01	0.51	−0.62	9.83
BD +12 2547 (AGB)	2.657	−3.439	378.07	5.98	0.60	0.10	9.20
BD +17 3248 (RHB)	1.259	−2.355	493.93	6.81	0.65	−0.53	9.05
BD +18 2757 (AGB)	1.400	−2.508	384.43	6.55	0.32	−0.03	8.05
BD +18 2976 (RGB)	2.218	−3.175	464.10	5.96	0.32	−0.12	8.54
BD −03 5215 (RHB)	1.095	−2.157	583.84	7.06	1.91	−0.90	9.36
HD 083212 (RGB)	1.361	−2.467	464.86	6.20	0.23	−0.44	8.07
HD 101063 (SGB)	0.813	−1.744	1312.86	6.26	40.43	−2.65	11.31
HD 107752 (AGB)	2.171	−3.145	524.60	6.08	0.42	−0.37	9.05
HD 110885 (RHB)	0.909	−1.897	575.05	7.06	2.02	−0.97	9.17
HD 111721 (RGB)	1.379	−2.486	609.98	6.42	0.61	−0.91	9.64
HD 115444 (RGB)	1.919	−2.965	491.45	6.14	0.34	−0.31	8.82
HD 119516 (RHB)	0.945	−1.950	483.08	7.24	1.29	−0.62	8.80
HD 121135 (AGB)	1.071	−2.126	502.44	6.68	0.60	−0.68	8.47
HD 122956 (RGB)	1.218	−2.309	444.03	6.27	0.19	−0.38	7.88
HD 135148 (RGB)	1.033	−2.075	394.38	5.99	0.07	−0.25	6.98
HD 195636 (RHB)	0.350	−0.779	435.75	7.64	3.38	−0.88	7.90
EV Lac (M3.5 V)	0.0706	−0.313	3005.25	2.72	2641	−8.15	17.78
V Hya (N6, AGB)	0.609	−1.367	12.39	6.53	1.32	1.99	4.34
89 Her (F2 Ib)	27.12	−4.000	43.59	1.92	0.05	−1.11	11.15

Download table as:  ASCIITypeset image

Figure 10 compares the theoretical and measured mass loss rates with one another as a function of L_*. For the four binary systems listed in Table 2 with combined A+B mass loss rates, we separated the measured value into two pieces using the modeled $\dot{M}$ ratio for the two components. (This was done for this figure only because the measured rates are shown as a function of a single star's luminosity.) For stars with only vsin i rotation period estimates, we used Equation (40) to compute $\dot{M}$ for the central plotting symbol and the entries in Table 3, and we recomputed $\dot{M}$ for the lower and upper limits given by Equation (41) to obtain the error bars. Figure 10(b) shows a comparison with the semi-empirical scaling law proposed by Schröder & Cuntz (2005), with

$$\begin{equation} \dot{M}_{\rm SC} \, = \, \eta \frac{L_{\ast } R_{\ast }}{M_{\ast }} \left(\frac{T_{\rm eff}}{4000 \, \mbox{K}} \right)^{3.5} \left(1 + \frac{g_{\odot }}{4300 \, g} \right)\!, \end{equation} \tag{ 42 }$$

where L_*, R_*, and M_* are assumed to be in solar units. For this plot we computed the normalization constant η such that the average modeled mass loss rate would equal the average measured mass loss rate for all 47 stars. Averages were taken using the logarithm of $\dot{M}$ so that all stars would contribute to the average comparably. We found η = 8.5 × 10⁻¹⁴ M_☉ yr⁻¹, which is within the error bars of the Schröder & Cuntz (2005) value.

Zoom In Zoom Out Reset image size

Overall, our "standard model" (i.e., Equation (35) with α₀ = 0.5, h = 0.5, and θ = 1/3) appears to match the measured mass loss rates reasonably well. We emphasize that this model does not contain any arbitrary η normalization factors. For the three test-case stars at the bottom of Table 2, however, our model does not do as well. The model underpredicts the mass loss from the dM flare star EV Lac by at least four orders of magnitude, and it also fails for the F supergiant 89 Her by a slightly smaller amount. For EV Lac and other flare-active M dwarfs, it is possible that coronal mass ejections and other episodic sources of energy (Mullan 1996) could be responsible for the bulk of the observed mass loss. For the carbon star V Hya, the reasonably good agreement between the model and measurements is probably a coincidence, since our model does not include the dusty radiative transfer or strong radial pulsations that are likely to be important for AGB stars.

It is interesting to highlight the case of the moderately rotating K dwarfs Eri, 70 Oph, and 36 Oph, which Holzwarth & Jardine (2007) found to have anomalously high mass loss rates. They concluded that the observed magnetic fluxes for these stars were insufficient to produce their dense outflows. For these stars, our modeled mass loss rates tended to be about a factor of 10–20 below the measured values. However, these models were computed using f_min from Equation (38). The measured values of f_* given in Table 1 for these stars are larger than their corresponding f_min values by factors ranging from 3 to 20. If instead these values were used, our modeled mass loss rates would be in better agreement with the astrospheric observations of Wood et al. (2005a).

We also developed a statistical measure of how well a given model agrees with the measured database of mass loss rates. We defined a straightforward least-squares parameter

$$\begin{equation} \chi ^{2} \, = \, \frac{1}{N} \sum _{i=1}^{N} [ \log \dot{M}_{i} (\mbox{model}) - \log \dot{M}_{i} (\mbox{obs}) ]^{2}, \end{equation} \tag{ 43 }$$

where the total number of comparisons (N = 40) excludes the final three test cases in Table 2 and counts each of the four A+B binaries as one. When χ² < 1, then (on average) the modeled and measured mass loss rates are within an order of magnitude of one another. Table 4 summarizes the results, including comparisons with other published empirical prescriptions. The η normalization factors for each of these scaling laws were computed similarly as the factor in Equation (42) above. Note that our standard model appears to be a significant improvement over both the popular Reimers (1975, 1977) and Schröder & Cuntz (2005) scalings.

Table 4. Mass Loss "Goodness of Fit"

Model	χ2
This paper (standard model)	0.650
This paper (ZAMS τc)	1.575
This paper (h = 0.25)	0.794
This paper (h = 1)	0.564
This paper (h = 3)	0.504
This paper (θ = 0.2)	0.620
This paper (θ = 0.5)	0.707
This paper (all [Fe/H] =0)	0.647
This paper ($\dot{M} = \dot{M}_{\rm hot} + \dot{M}_{\rm cold}$)	0.703
This paper (α from Trampedach & Stein 2011)	0.770
Reimers (1975, 1977)	1.260
Mullan (1978), Equation (4a)	3.768
Nieuwenhuijzen & de Jager (1990)	2.356
Catelan (2000), Equation (A1)	1.924
Schröder & Cuntz (2005)	1.131

Download table as:  ASCIITypeset image

To further explore the proposed model, we varied some of the modeling parameters described in Section 3. Varying the TR filling factor exponent θ did not have much of an effect on χ². However, varying the flux height scaling factor h did change χ² significantly. We found that a larger value of h ≈ 3 gives much better agreement with the measured mass loss rates than does the standard value of h = 0.5 (see Table 4). We decided not to adopt this larger value, though, because it falls well outside the range of empirically determined h values for the Sun's corona.

We also tried removing some of the imposed complexity of the standard model to see if simpler assumptions would give adequate results. Removing the gravity-dependent modification factor of (g_☉/g)^0.18 from the definition of the convective overturn time resulted in significantly poorer agreement with the data (i.e., more than double the χ² of the standard model). We explored the importance of metallicity by replacing the published [Fe/H] by purely solar values ([Fe/H] =0). This actually improved the value of χ² from the standard model, but only by <1%. Removing the exponential factor in Equation (35) gave a slightly higher value of χ² (8% larger than the standard model).

We also noted that the theoretical photospheric Alfvén wave fluxes from Musielak & Ulmschneider (2002a) exhibited a strong dependence on the convective mixing length parameter (i.e., F_A*∝α^2.1). Thus, instead of simply assuming α = 2 as in the standard model, we created a linear regression fit to the tabulated simulation results of Trampedach & Stein (2011), who found empirical values of α between 1.6 and 2.2 depending on T_eff, log  g, and M_*. We used the following approximate fit

$$\begin{equation} \alpha _{\rm TS} \, \approx \, 1.91 - \frac{T_{\rm eff}}{6181 \, \mbox{K}} + \frac{\log \, g}{5.58} + \frac{M_{\ast }}{19.1 \, M_{\odot }} \end{equation} \tag{ 44 }$$

and did not allow α_TS to be less than 1.6 or greater than 2.2. Thus, we multiplied the value of F_A* from Equation (7) by a factor of (α_TS/2)^2.1. The predicted mass loss rates for the Table 2 stars had about an 18% higher value of χ² than the standard model, so we did not pursue this mixing length prescription any further.

For additional context about the hot and cold mass loss models described in Section 3.3, Figure 11 shows the ratio $\dot{M}_{\rm hot} / \dot{M}_{\rm cold}$ for the 47 modeled stars as well as the TR Mach number M_{A, TR} = u_TR/V_{A, TR}. It is clear that the dwarf stars are dominated by hot coronae, and the stellar wind outflow is still negligibly small at the coronal base. However, as the luminosity exceeds ∼50 L_☉ for the giant stars, the hot coronal contribution goes away and the acceleration becomes dominated by wave pressure.

Zoom In Zoom Out Reset image size

Figure 12 examines how the plasma number density at the TR, n_TR = ρ_TR/m_H, depends on stellar rotation. Holzwarth & Jardine (2007) assumed n_TR∝Ω^0.6∝P^−0.6_rot (see also Ivanova & Taam 2003). Although our models do not follow a single universal relation for both giants and dwarfs, the proposed scaling (or one slightly steeper) may be appropriate for certain sub-populations of stars. For the dwarf stars with well-determined rotation periods, it is interesting that the Sun's computed value of n_TR is larger than that of stars having higher magnetic activity. Equation (22) shows that the dependence on filling factor is weak (i.e., about f^0.19_* for θ = 1/3), so the variation comes mostly from the other stellar parameters. This seems to stand in contrast with other observational determinations of coronal electron densities, where n tends to increase with activity (Güdel 2004). However, X-ray determinations of number density are probably dominated by closed-field active regions, which are not necessarily correlated with the regions driving the stellar wind.

Zoom In Zoom Out Reset image size

In addition to the above comparisons with the individual stars of Table 2, we also created some purely theoretical sets of stellar models and computed $\dot{M}$ for them. We began with the ZAMS model parameters given by Girardi et al. (2000), and we assumed solar metallicity for a range of constant rotation rates. This gave rise to a two-dimensional grid of models (varying T_eff and P_rot) for main-sequence stars. Because the modeled stars are all high-gravity dwarfs, we used only Equation (36) for τ_c, in combination with Equation (38) for f_min as a function of Rossby number.

Figure 13 shows the resulting mass loss rates as a function of T_eff and P_rot. If we had not utilized a basal "floor" on f_*, we would have predicted a steep dropoff in mass loss for T_eff ≳ 7000 K, at which the Gunn et al. (1998) expression for τ_c decreases rapidly. However, because of the floor, there appear to be reasonably strong mass loss rates up to the point at which subsurface convection zones disappear at T_eff ≳ 9000 K. There is a slightly discontinuous dip in the predicted basal mass flux around T_eff ≈ 6100 K that arises because of the iteration for ${\cal R}$, ρ_TR, and Q_TR. If the calculation of these quantities is halted after only one iteration, the final value of $\dot{M}$ varies more smoothly as a function of T_eff. We plan to utilize a more self-consistent non-WKB model of Alfvén wave reflection in future versions of this work.

Zoom In Zoom Out Reset image size

The mass loss rates shown in Figure 13 are almost all due to the hot coronal processes discussed in Section 3.1. Thus, it is possible to simplify the components of Equation (26) in order to obtain an approximate scaling relation for $\dot{M}$ that is reasonable for these main sequence stellar models. Ignoring the weakest dependences on some stellar parameters (i.e., factors with exponents less than or equal to 1/7), we found

$$\begin{eqnarray} \frac{\dot{M}}{10^{-10} \, M_{\odot }\, \mbox{yr}^{-1}} \, &\sim &\, \left(\frac{R_{\ast }}{R_{\odot }} \right)^{16/7} \left(\frac{L_{\ast }}{L_{\odot }} \right)^{-2/7}\nonumber\\ && \times \left(\!\frac{F_{\rm A \ast }}{10^{9} \, \mbox{erg} \,\, \mbox{cm}^{-2} \,\, \mbox{s}^{-1}} \! \right)^{12/7} f_{\ast }^{(4 + 3\theta)/7} \,,\quad\quad \end{eqnarray} \tag{ 45 }$$

which reproduces the curves in Figure 13 to within about an order of magnitude. Despite the fact that this scaling formula is relatively easy to apply, we do not recommend its use in stellar evolution or population synthesis calculations. Once stars leave the main sequence, Equation (45) is no longer a good approximation.

We also computed a time-dependent mass loss rate for the evolutionary track of a star having M_* = 1 M_☉. There is evidence that the wind from the "young Sun" was significantly denser than it is today, and this more energetic outflow may have been important to early planetary evolution (e.g., Wood 2006; Güdel 2007; Sterenborg et al. 2011; Suzuki 2011). We used the BaSTI evolutionary track plotted in Figure 1 (Pietrinferni et al. 2004) for the time variation of R_* and L_*. We grafted on a model of rotational evolution for a solar-mass star from Figure 6(a) of Denissenkov et al. (2010). For late ages (t ≳ 100 Myr or log t ≳ 8), this model has approximately P_rot∝t^0.54. Such an age scaling is well within the range of empirically determined power laws (t^0.5 to t^0.6) obtained from young solar analogs (e.g., Barnes 2003; Güdel 2007).

Figure 14(a) shows how the luminosity and two dimensionless parameters related to the rotational dynamo (Ro and f_min) vary as a function of age for this model. Note that prior to about t ≈ 70 Myr the Rossby number is small enough that the filling factor appears to be saturated near its maximum assumed value of 0.5. At very late times, when the star begins to ascend the red giant branch, the Rossby number decreases again because of the increase in τ_c with decreasing T_eff and gravity. We utilized the (g_☉/g)^0.18 correction factor when computing τ_c, but it was relatively unimportant until the star left the main sequence.

Zoom In Zoom Out Reset image size

Figure 14(b) gives our prediction for the age variation of the a solar-type star's mass loss rate. For ages between about t ≈ 0.2 and 7 Gyr the decrease in mass loss appears to be fit approximately by a power law, with $\dot{M} \propto t^{-1.1}$. This is a significantly shallower age dependence than the t⁻² decline suggested by Wood et al. (2002) on the basis of astrosphere measurements. We note that if the rotation period was the only variable to change with time, Equation (45) would give something like $\dot{M} \propto P_{\rm rot}^{-2.4}$ (for f_* = f_min and θ = 1/3). Thus, a more rapid increase of P_rot with age—such as the t^0.85 dependence in the solar-mass rotational model of Landin et al. (2010)—would give rise to a steeper age–$\dot{M}$ relationship more similar to that of Wood et al. (2002).

For comparison, Figure 14(b) also shows that the Schröder & Cuntz (2005) scaling law predicts a much smaller range of mass loss variation for the young Sun than does the present model. We also show a model (Cranmer 2008, 2009) and measurements (Hartigan et al. 1995) for classical T Tauri stars (CTTSs) at the youngest ages. It is clear that for t ≲ 10 Myr some additional physical processes must be included (e.g., accretion-driven turbulence on the stellar surface) to successfully predict mass loss rates.