The Wolf-Rayet Stellar Response To The Iron Opacity Bump: Envelope Inflation, Winds, and Microturbulence

TL;DR

This study investigates how Wolf-Rayet stars respond to iron opacity bumps, revealing that most do not have inflated envelopes, and identifies conditions for wind-driven mass loss and sonic points related to the iron opacity bump.

Contribution

The paper provides semianalytical estimates and models showing the conditions under which Wolf-Rayet stars develop winds and sonic points due to iron opacity bumps, clarifying their envelope structures and mass-loss rates.

Findings

Most early-type WR stars lack inflated envelopes due to winds or low luminosity.

A minimum mass-loss rate is required for an opacity bump to harbor a sonic point.

WR winds driven by the iron opacity bump must have mass-loss rates above ~10^{-6} M_sun/yr.

Abstract

Early-type Wolf-Rayet(WR) stellar models harbor a super-Eddington layer in their outer envelopes due to a prominent iron opacity bump. In the past few decades, one-dimensional hydrostatic and time-steady hydrodynamic models have suggested a variety of WR responses to a super-Eddington force including envelope inflation and optically thick winds. In this paper, I study these responses using semianalytical estimates and WR models from both MESA and Ro & Matzner; four conclusions are present. First, early-type WR stars do not harbor inflated envelopes because they have either strong winds or insufficient luminosities. Second, the condition for an opacity bump to harbor a sonic point is expressible as a minimum mass-loss rate, $\dot{M}_{\rm sp}(L_*)$. In agreement with Grassitelli et al. and Ro, the majority of galactic early-type WR stars can harbor sonic points at the iron opacity bump. However, about half of those in the Large Magellanic Cloud cannot given typical wind parameters. Third, WR winds driven by the iron opacity bump must have mass-loss rates that exceed a global minimum of $10^{-5.8}-10^{-6}M_\odot\,{\rm yr}^{-1}$. Lastly, the observed early-type WR distribution follows a simple mass-loss relation derived here if the radiation-to-gas pressure ratio is approximately $p_r/p_g \sim 145$ in the wind; a value consistent with studies by Grafener et al. and Nakauchi & Saio.

Equations

$${\Gamma }_r\ge {\Gamma }_c\equiv \frac{8(4-3\beta )(1-\beta )}{8(4-3\beta )-3{\beta }^2}<1,$$

$$\frac{L_c}{L_{\ast }}\lesssim 0.04{\left(\frac{M_{\ast }}{25{\mathrm{M}}_{\odot }}\right)}^{-0.45}{\left(\frac{\rho }{10^{-10}\mathrm{g}\mathrm{c}{\mathrm{m}}^3}\right)}^{-0.5}.$$

$$\frac{1}{\rho }\frac{d{p}_r}{dr}=-\frac{{\Gamma }_r{v}_k^2}{r},$$

$$\frac{1}{\rho }\frac{d{p}_g}{dr}=-\frac{v_k^2(1-{\Gamma }_r)}{r},$$

$$\kappa (\rho ,T)\simeq {\kappa }_{\mathrm{Edd}}\equiv \frac{4\pi cG{M}_{\ast }}{L_{\ast }};$$

$$\frac{d\log (\rho )}{d\log (T)}=\frac{1-{\Gamma }_r}{\varphi {\Gamma }_r}-1,$$

$${\kappa }^{\prime }\equiv \frac{d\log (\kappa )}{d\log (r)}={\kappa }_{\rho }{\rho }^{\prime }+{\kappa }_T{T}^{\prime },$$

$${\kappa }_{\rho }\equiv \frac{\partial \log (\kappa )}{\partial \log (\rho )},{\kappa }_T\equiv \frac{\partial \log (\kappa )}{\partial \log (T)}.$$

$$r\simeq R_{\ast }{\left(1-R_{\ast }{\int }_{T(R_{\ast })}^{T(r)}\frac{4a_r{T}^3}{3G{M}_{\ast }{\rho }_{\mathrm{Edd}}(T)}dT\right)}^{-1}.$$

$$\alpha \equiv \frac{H_T}{r}=-\frac{T}{r}\frac{dr}{dT}=\frac{4p_r}{\rho {\Gamma }_r{v}_k^2},$$

$$\frac{d}{dr}(r^2\rho v)=0,$$

$$v\frac{dv}{dr}+\frac{1}{\rho }\frac{d{p}_g}{dr}=-\frac{v_k^2(1-{\Gamma }_r)}{r},$$

$$v^{\prime }=\frac{v_k^2\left(q_g+{\Gamma }_r(1+\varphi )-1\right)}{v^2-c_g^2},$$

$$v(R_{\mathrm{sp}})=c_g(R_{\mathrm{sp}})\simeq 33{\mathrm{km}\mathrm{s}}^{-1},$$

$${\Gamma }_r(R_{\mathrm{sp}})=\frac{1-q_g}{1+\varphi }<1.$$

$$\dot{E}=L_r+\dot{M}\mathcal{B},$$

$$w=\frac{5}{2}\frac{p_g}{\rho }+\frac{4p_r}{\rho },$$

$${\kappa }_{\mathrm{sp}}\equiv \frac{4\pi cG{M}_{\ast }(1-q_{g,\mathrm{sp}})}{L_{\ast }(1+{\varphi }_{\mathrm{sp}})}\simeq {\kappa }_{\mathrm{Edd}},$$

$$6{\Psi }^2$$

$$\xi =\frac{\dot{M}{v}_k^2}{2L_r},\mathcal{W}=\frac{1-q_g}{q_g},\Psi =\frac{\varphi }{1+\varphi }\mathcal{W},$$

$${\kappa }_{\rho }=\frac{d\log (\kappa )}{d\log (\rho )},{\kappa }_T=\frac{d\log (\kappa )}{d\log (T)}$$

$$q_g=\frac{2p_g}{\rho v_k^2},q_r=\frac{2p_r}{\rho v_k^2}\simeq \frac{\alpha }{2},\varphi =\frac{p_g}{4p_r},$$

$${\kappa }_{\rho }+\Psi {\kappa }_T<0.$$

$${\dot{M}}_{\mathrm{sp}}\equiv 4\pi R_{\mathrm{sp}}^2{\rho }_{\mathrm{sp}}{v}_{\mathrm{sp}}.$$

$${\dot{M}}_b\equiv 4\pi R_b^2{\rho }_b{v}_b.$$

$$v_b^{\prime }=-T^{\prime }\left(\frac{{\rho }^{\prime }}{T^{\prime }}\right)-2\simeq \frac{4}{\alpha {\Gamma }_r}-2,$$

$$v_b^2\simeq c_g^2+\frac{2c_g^2+v_k^2\left({\Gamma }_r(1+\varphi )-1\right)}{4(\alpha {\Gamma }_r{)}^{-1}-2}.$$

$${\Gamma }_r\ge \left(1+\frac{v^2}{v_k^2}{v}^{\prime }\right){\Gamma }_c.$$

$$\frac{{\tau }_0}{{\tau }_b}$$

$$t_{\mathrm{CAK}}=\frac{{\sigma }_e{v}_{\mathrm{th}}\rho }{dv/dr}\gtrsim 1,$$

Full text

1. Introduction

The inclusion of fine iron M-shell transitions in the Rosseland opacity calculations led to the discovery of a prominent opacity bump by the OPAL and the Opacity Project groups (Rogers & Iglesias, 1992; Seaton et al., 1994). Situated at $T\simeq 10^{5.2}$ K, the iron bump resolves many outstanding issues regarding the $\kappa$-mechanism in pulsating stars such as the “bump and beat” mass discrepancy in Cepheids (Moskalik et al., 1992), oscillations in $\beta$ Cephei and slowly pulsating B-type stars (Dziembowski & Pamiatnykh, 1993; Charpinet et al., 1997). For Wolf-Rayet (WR) stars, the iron bump has been suggested to trigger a variety of responses including envelope inflation (Ishii et al., 1999; Petrovic et al., 2006; Gräfener et al., 2012), strange-mode oscillations (Glatzel et al., 1993; Kiriakidis et al., 1996; Glatzel et al., 1999; Glatzel & Kaltschmidt, 2002), turbulence (Grassitelli et al., 2016), and quasi-steady winds (Nugis & Lamers, 2002; Gräfener & Hamann, 2005; Gräfener et al., 2017; Nakauchi & Saio, 2018; Grassitelli et al., 2018). These responses are challenging to test observationally since WR winds are optically thick which renders basic stellar parameters invisible (e.g., mass, radius, rotation).

WR spectral analyses employ CMFGEN (Hillier & Miller, 1998) and PoWR codes (Hamann, 1985; Gräfener et al., 2002) to study WR interiors. These codes compute radiative transfer calculations in the co-moving frame prescribed by a supersonic $\beta$-law velocity structure $v=v_{\infty}(1-R_{20}/r)^{\beta}$ (Castor et al., 1975). The radiative acceleration is treated self-consistently in PoWR (Sander et al., 2015). These studies find roughly half of stars have ‘hydrostatic’ radii, $R_{20}$, (defined where the Rosseland optical depth is $\tau\sim 20$) exceed stellar model predictions, $R_{*}$, by up to an order of magnitude. Consequently, the ‘hydrostatic surface’ temperatures, $T(R_{20})=(L_{*}/(4\pi\sigma R_{20}^{2}))^{1/4}$, and escape speeds, $v_{\rm esc}(R_{20})$, are significantly smaller than expected. This is especially problematic if we assume hydrogen-free stars are on the helium-burning main-sequence since their structures are simple (Langer, 1989; Schaerer & Maeder, 1992).

The stellar radius is not the only discrepancy seen in WR stars. WC stars are a WR sub-type with strong carbon spectral line features. Evolutionary models predict WC luminosities exceed $10^{5.4}$ $\mathrm{L}_{\odot}$ (Meynet & Maeder, 2005; Georgy et al., 2012) despite observations of galactic WC stars with $4.9\lesssim\mathrm{log_{10}}(L_{*}/\mathrm{L}_{\odot})\lesssim 5.4$ (Sander et al., 2012; Yoon, 2017). Stellar models predict the luminosities of early-type WN stars (which have nitrogen-dominant spectral lines) establishes the spectroscopic properties (Langer, 1989), which is not seen in observations (Hamann et al., 2006). While the luminosity discrepancy is significant, this investigation considers only the radius discrepancy problem.

Sub-surface convection is an inefficient transport mechanism in WR stellar models (e.g., see Eq. 2 and 29). As a result, the iron bump forces a hydrostatic envelope to develop a rarified layer of near-Eddington gas that is extended by several core-radii (Ishii et al., 1999) and, in some instances, a massive, dense, super-Eddington shell (Joss et al., 1973; Sanyal et al., 2015). A number of authors have proposed ‘envelope inflation’ as a partial resolution to the WR radius discrepancy problem (Kato & Iben, 1992; Schaerer, 1996; Heger & Langer, 1996; Petrovic et al., 2006; Gräfener et al., 2012). If realized, this phenomenon would have several consequences for stellar and transient observations. For instance, there would be earlier binary interactions (Kruckow et al., 2016), redder massive stars (Gräfener et al., 2012; Sanyal et al., 2015; McClelland & Eldridge, 2016) and ‘circumstellar’ imprints in supernova lightcurves and spectra (Moriya et al., 2015; Dessart et al., 2018).

An alternative response to a near- or super-Eddington force is an outflow. Nugis & Lamers (2002) argue the iron bump in WR envelopes can harbour a sonic point. The luminosity-mass relation for early-type WR stars is a good approximation for the radiative luminosity beneath the sonic point at the iron bump (Langer 1989; Schaerer & Maeder 1992; RM16) which permits the Eddington limit (or sonic point condition) near the surface to be written in terms of a critical opacity $\kappa_{*}=4\pi cGM_{*}L_{*}^{-1}$. Ro (2017) and Grassitelli et al. (2018) realize the critical opacity cannot be reached for arbitrarily small densities. Therefore, a sonic point exists above a minimum density and mass loss rate, $\dot{M}_{\rm sp}$. In agreement, these authors find galactic early-type WR stars can harbour a sonic point at the iron bump (i.e., $\dot{M}_{\rm gal}\gtrsim\dot{M}_{\rm sp}$). Both authors use the approximate sonic point condition from Nugis & Lamers (2002), which is that sonic points rest above the iron opacity peak (i.e., $\partial\kappa/\partial T>0$). In this paper, I revisit this calculation of $\dot{M}_{\rm sp}$ with the sonic point condition from Ro & Matzner (2016).

Petrovic et al. (2006) and Ro & Matzner (2016) (hereafter called RM16) question whether envelope inflation, a hydrostatic phenomenon, can manifest in a wind. RM16 found early-type WR stellar models with mass loss rates below $\dot{M}_{b}$ have structures that resemble inflated envelopes; although, the regions contain highly supersonic velocity fields. Since the authors consider only continuum-driven winds, additional line effects (e.g., Doppler-shifts and line-deshadowing) are neglected in the outflow acceleration; thus, $\dot{M}_{b}$ is likely an overestimate. Wind models with $\dot{M}>\dot{M}_{b}$ have sufficiently large densities and small pressure scale heights such that envelope inflation does not manifest.

In this paper, I compare the critical mass loss relations, $\dot{M}_{\rm sp}$ and $\dot{M}_{b}$, with observations of early-type WR stars in the Milky Way and Large Magellanic Cloud. I begin with an explanation for why radiative envelopes inflate (§2.1 and 2.1.2) followed by a semi-analytic description of their structure (§2.1.1). I review the equations for a continuum-driven wind in §2.2 along with the conditions for a sonic point. The critical mass loss relations for a sonic point and envelope inflation are derived in §2.2.2 and 2.2.3. Results and discussions follow in §3 and 4.

2. Stellar Responses to an Iron Opacity Bump

We begin by deriving the conditions for envelope inflation (§2.1) and winds (§2.2). These conditions may suggest there are only two stellar responses to a prominent iron bump; this is not the case. There are four possible responses partitioned by two critical mass loss relations in $(\dot{M},L_{*})$-space.

2.1. Hydrostatic Response

The Rosseland opacity tables show multiple bumps due to the partial ionization of Fe ($T_{\rm Fe}\equiv 10^{5.2}$ K), He ($10^{4.6-4.8}$ K), and H ($10^{4}$ K). An opacity bump increases the local Eddington ratio, $\Gamma_{r}=\kappa L_{r}/4\pi cGM$, until the onset of convection:

[TABLE]

where $\beta=p_{g}(p_{g}+p_{r})^{-1}$ is the ratio of gas, $p_{g}$, and radiation pressures, $p_{r}$ (and not related to the $\beta$-law). Convection reduces the radiative luminosity, $L_{r}$, until the opacity decreases or gas becomes optically thin.

Hydrogen-free WR stars are expected to be massive (7 - 35 $\mathrm{M}_{\odot}$), very luminous ($10^{4.8}$ - $10^{6}\,\mathrm{L}_{\odot}$), and compact ($\sim$ 1 $\mathrm{R}_{\odot}$) with radiation-dominated envelopes ($\beta\lesssim 0.3$) (Langer, 1989; Schaerer & Maeder, 1992). As a result, convective instability occurs at high Eddington ratios, $\Gamma_{c}\gtrsim 0.7$. The maximum convective luminosity $L_{c}\simeq 4\pi r^{2}v_{c}U$ from subsonic transport $v_{c}<c_{s}\simeq\sqrt{4p_{r}/(3\rho)}$ of thermal energy (density) $U\simeq 3p_{r}$ is a small fraction of the total WR luminosity,

[TABLE]

The above relation is evaluated using the empirical WR relations in Appendix A. Note that the opacity tables are unavailable for $\rho\lesssim 10^{-10.5}$ g cm3 near iron bump temperatures. Wave transport may contribute an additional $L_{w}\sim 0.5L_{c}$, although strong shock formation is inevitable at this point (Ro & Matzner, 2017). As a result, WR envelopes are effectively radiative in structure.

The WR envelope can be well approximated by the equations of continuity, radiative diffusion approximation

[TABLE]

and hydrostatic pressure balance,

[TABLE]

where $p_{r}=a_{r}T^{4}/3$ is the radiation pressure with $a_{r}$ being the radiative constant. Because convection in inefficient, the Eddington ratio, $\Gamma_{r}\simeq\kappa(\rho,T)L_{*}/(4\pi cGM_{*})$, becomes a function of density and temperature only. The Eddington limit (i.e., $\Gamma_{r}=1$) is a contour of constant opacity in $(\rho,T)$-space,

[TABLE]

I refer to this as the ‘Eddington contour’. The Eddington contour shifts to lower densities and gas pressures for higher luminosity-mass ratios, which increases with mass (see Schaerer & Maeder 1992, Langer 1989, or Appendix A). The dashed lines in Fig. 1 shows two solutions of the Eddington contour for $M_{*}=10$ and 25 $\mathrm{M}_{\odot}$ helium stars with solar metallicity. Luminosities are defined using Eq. (29).

Consider the ratio of Eq. (3) and (4) in the following form:

[TABLE]

where $\phi\equiv p_{g}/4p_{r}=\beta/[4(1-\beta)]$. Since the right hand side is a function of only density and temperature, the envelope structure around the iron bump is well-approximated by an ordinary differential equation in one variable, $\rho(T)$. Fig. 1 shows solutions to Eq. (6) for a range of densities starting at $T=10^{5}$ K. These solutions are qualitatively similar to Fig. 10 of G12 which shows envelope solutions for various outer boundary conditions. To understand the variety of solutions, let us consider Eq. (6) in more detail.

Suppose the envelope solution does not approach the Eddington contour for $T\gtrsim 10^{5.2}$ K and remains sub-Eddington. Near the iron peak, the Eddington contour extends to low densities and, consequently, small $\phi\ll 1$. A strongly radiation-dominated, sub-Eddington fluid demands the left side of Eq. (6) to be very large. As a result, $\rho\rightarrow 0$ and an atmosphere forms. These solutions are not shown in Fig. 1 since the divergence is extremely rapid and numerically challenging to resolve accurately.

If the envelope solution instead approaches the Eddington contour then a density inversion forms wherever $\Gamma_{r}=\Gamma_{\rm inv}\equiv(1+\phi)^{-1}<1$ (Joss et al., 1973). The opacity increases and the envelope becomes super-Eddington since $\Gamma_{\rm inv}\simeq 1-\phi\sim 1$. Both the density and Eddington ratio continue to rise until either $\phi$ is not small or the opacity gradient becomes negative:

[TABLE]

where the prime indicates a logarithmic radial derivative, ${}^{\prime}=d\mathrm{log}/d\mathrm{log}(r)$, and,

[TABLE]

This set of solutions is shown in Fig. 1 as solid lines.

A third solution may appear to exist where the Eddington ratio is neither small and a density inversion does not form: $0\ll\Gamma_{r}<\Gamma_{\rm inv}<1$. This is not possible beyond the iron peak since the Eddington contour swings to increasing density. There, the envelope must form a density inversion to continue tracing the Eddington contour, which forces $\Gamma_{r}\geq\Gamma_{\rm inv}$. In principle, these arguments apply to the hydrogen and helium opacity bumps under the same physical conditions: radiation-dominated, near-Eddington, and effectively radiative envelope. This is not seen in some models of massive main-sequence stars by Sanyal et al. (2015) since convection is not a negligible energy transport mechanism.

2.1.1 An Approximate Envelope Solution

Fig. 1 shows the structure of a hydrostatic, radiation-dominated, near-Eddington envelope with inefficient convective transport is well-approximated by the Eddington contour (dashed lines). Because the density, $\rho_{\rm Edd}(T)$, is a function of temperature along the Eddington contour (see Eq. 5), the temperature structure can be approximated for with the integration of Eq. (3) given the radius, $R_{*}$, and temperature $T(R_{*})$, at an inner boundary:

[TABLE]

Note that $\Gamma_{r}$ is one along the Eddington contour. Eq. (5) and (9) provide a solution for the density structure $\rho_{\rm Edd}(r)$.

Fig. 2 shows hydrostatic solutions for a 23 $\mathrm{M}_{\odot}$ pure helium star with solar metallicity from MESA and G12. A MESA inlist is available in Appendix B. The core radius is defined where $R_{*}=r(T=10^{5.6}\,{\rm K})$ in MESA. Also shown is a steady wind solution with $\dot{M}=10^{-5.2}$ $\mathrm{M}_{\odot}{\rm yr}^{-1}$from RM16. RM16 categorize this wind solution as ‘weak’ and ‘extended’ since the inertial terms do not affect the envelope structure significantly. It is clear from Fig. 2 that Eq. (9) is an excellent approximation for all of these solutions.

2.1.2 Envelope Inflation

RM16 found their wind solutions do not resemble inflated envelopes so long as the temperature scale height (see Eq. 3),

[TABLE]

remains small (i.e., $\alpha<1$). Here, we examine the quality of $\alpha$ as a measure of envelope inflation.

A sequence of $M_{*}=5-50\,\mathrm{M}_{\odot}$ hydrogen-free WR stellar models with solar metallicity ($Z_{\odot}=0.02$) and GN93hz abundances are solved for with the stellar evolution code, MESA (Paxton et al., 2018). Two chemical compositions are considered to represent WNE and WC interiors: (1) pure helium ($Y=1-Z$, $\mu=4/3$), and (2) a mixture of carbon, oxygen, and helium ($X_{\rm C}+X_{\rm O}+Y=1-Z$, $X_{\rm C}=0.4$, $X_{\rm O}=0.1$, $\mu\simeq 1.5$). I construct a second sequence of lower metallicity models ($Z_{\odot}=0.01$) to represent WNE and WC stars in the LMC. Empirical mass-radius-luminosity relations are constructed from these models at temperatures $T=10^{5.6}$ K well below the iron bump which are available in Appendix A. Solutions for Eq. (10) are shown in Fig. 3 for pure helium models with solar metallicity between $M_{*}=$15 - 23 $\mathrm{M}_{\odot}$ and fixed temperature domains, $T=10^{5}$ - $10^{5.6}$ K. These solutions are found using Eq. (5) and (9) with the empirical relations as inner boundary conditions.

Inferred radii from observed WR stars are found to exceed models without envelope inflation by factors of 3 or more (i.e., $R_{20}/R_{*}\gtrsim 3$). WR models with max($\alpha)>1$ have radii $R_{*}\gtrsim 1.4\,\mathrm{R}_{\odot}$ at $T=10^{5}\,$K. If stellar models with max($\alpha)>1$ are heuristically defined to be ‘inflated’ then the definition is conservatively small in regards to the radius discrepancy problem. Fig. 3 shows helium models with solar metallicity above $M_{b}\simeq 19\,\mathrm{M}_{\odot}$ or luminosity $L_{b}\simeq 10^{5.7}\,\mathrm{L}_{\odot}$ are considered to be inflated here. Using the relations from Schaerer & Maeder (1992) instead of those in Appendix A gives a minimum mass of $M_{b}\simeq 14\,\mathrm{M}_{\odot}$ as found by RM16. The critical mass for envelope inflation is sensitive to the empirical relations despite how similar they may appear (e.g., Fig. 6). A criterion for envelope inflation is introduced by Grassitelli et al. (2018); although, a precise definition is not necessary for the purposes of this investigation.

2.2. Dynamic Response

Consider a steady, spherically symmetric wind,

[TABLE]

with a constant rate of mass loss rate, $\dot{M}=4\pi r^{2}\rho v$, where $v$ is the radial wind velocity. We assume the gas and radiation are in thermal equilibrium and the wind optical depth is sufficiently large to apply the radiative diffusion approximation. The opacity in Eq. (3) can be approximated for with the Rosseland mean opacity $\kappa_{R}(\rho,T)$; the limitations of this are discussed in the next section (§2.3) and by RM16.

Combining Eq. (3) and (11) with the momentum equation,

[TABLE]

leads to a wind equation,

[TABLE]

with dimensionless variable $q_{g}\equiv 2c_{g}^{2}/v_{k}^{2}\ll 1$, gas sound speed $c_{g}^{2}=p_{g}/\rho=k_{b}T/\mu$, Boltzmann constant $k_{b}$, and mean molecular weight $\mu$.

2.2.1 Sonic Point Conditions

The wind equation is critical where the flow speed matches the gas sound speed,

[TABLE]

rather than the total sound speed, $c_{s}^{2}=\gamma p/\rho\simeq 4p_{r}/(3\rho)$, where $\gamma$ is the adiabatic index. The location of the sonic point, $R_{\rm sp}$, is found where the numerator of Eq. (13) is zero:

[TABLE]

For time-steady flow, energy conservation states

[TABLE]

where $\dot{E}$ is the rate of stellar energy loss, $\mathcal{B}=w+v^{2}/2-v_{k}^{2}$ is the Bernoulli factor, and

[TABLE]

is the specific enthalpy. The radiative luminosity is approximately constant near WR sonic points since gas pressure is negligible and $\dot{M}\mathcal{B}\ll L_{*}$ (RM16). Therefore, a constant radiative luminosity, $L_{r}\simeq L_{*}$, is a good approximation. This simplifies Eq. (15) to become

[TABLE]

which is a contour in ($\rho,T$)-space in close proximity to the Eddington contour (Eq. 5) since $q_{g,{\rm sp}},\phi\ll 1$.

2.2.2 Minimum Mass Loss Relation for a Transonic Wind: $\dot{M}_{\rm sp}$

A wind solution requires the velocity gradient to be positive at the sonic point ($v^{\prime}_{\rm sp}>0$). RM16 use l’H$\hat{\mathrm{o}}$pital’s rule to solve Eq. (13) exactly:

[TABLE]

where,

[TABLE]

[TABLE]

and,

[TABLE]

are evaluated at the sonic point. WR winds are far from photon tiring limit at the sonic point and, so, $\xi\simeq 0$. Since $\phi$, $q_{g}$, ${\cal W}^{-1}$, and $\Psi\ll 1$ are small in the domain of interest, Eq. (19) is well approximated by the expression

[TABLE]

The roots of Eq. (20) represent a contour in $\rho$ and $T$ for a given opacity table. This contour is shown as the black line in Fig. 4 along with a black dotted contour of the iron peak (i.e., $\kappa_{T}=0$). Also shown is the Eddington contour for a 23 $\mathrm{M}_{\odot}$ with $L_{*}=10^{5.795}$ $\mathrm{L}_{\odot}$ from RM16. As argued by Nugis & Lamers (2002), the sonic point must reside on the hot side of the iron peak to drive a wind. This is not true for sufficiently small sonic point densities due to the factor of $\Psi\propto\alpha^{-1}\propto\rho$. Instead, Eq. (20) reaches a minimum density of log${}_{10}(\rho\,/{\rm g\,cm^{-3}})\simeq-10.5$ and follows a path of nearly constant ratio of radiation-to-gas pressure or $\alpha\sim 3$. This feature is relevant for WR stellar models with inflated envelopes where $\alpha>1$.

Given the stellar mass, luminosity and radius relations and an opacity table (i.e., chemical composition), the intersection of contours defined by Eq. (18) and (19) provide the minimum sonic point density, $\rho_{\rm sp}$, and temperature (or sound speed, $v_{\rm sp}$). The intersection is shown as a black circle in Fig. 4. Lamers & Cassinelli (1999) argues the hydrostatic solution (i.e., $\dot{M}=0$) is a good approximation for the subsonic wind structure. Numerical solutions by Petrovic et al. and RM16 show this to be accurate in WR winds. Therefore, Eq. (9) can be used to provide an estimate the sonic point radius, $R_{\rm sp}$, at the minimum sonic point density or temperature. Combining $\rho_{\rm sp}$, $v_{\rm sp}$, and $R_{\rm sp}$ provides an estimate for the minimum mass loss rate for a WR wind to harbour a sonic point:

[TABLE]

For comparison, Ro (2017) and Grassitelli et al. (2018) assume $\kappa_{T}=0$ at the minimum sonic point for wind-like solutions (i.e., $v^{\prime}_{\rm sp}\geq 0$). Grassitelli et al. (2018) compute hydrodynamic stellar structure models to support their approximation of the sonic point radius as $R_{\rm sp}\simeq 1$ $\mathrm{R}_{\odot}$. Ro (2017) use the empirical mass-radius relation from Schaerer & Maeder (1992) for the sonic point radius (i.e., $R_{\rm sp}\simeq R_{*})$.

I find the approximations described above are fine for WNE stellar models with masses below $M_{*}\lesssim 19\,\mathrm{M}_{\odot}$ (or $L_{*}\lesssim 10^{5.7}\,\mathrm{L}_{\odot}$) and solar metallicity. Above these masses (or luminosities), however, the effects from envelope inflation and additional term, $\kappa_{\rho}$, in Eq. (20) are important. For instance, the previous minimum sonic point condition, $\kappa_{T}=0$, no longer corresponds to wind-like solutions. Instead, the solutions decelerate at the sonic point and are subsonic at larger radii. These effects also cause Eq. (21) to invert with stellar mass which suggests there is a global minimum mass loss relation for WNE stellar models harbouring sonic points at the iron bump; a discussion on this is in §4.2.

2.2.3 Maximum Mass Loss Rate for Extended Winds

Fig. 4 shows a MESA model of an (hydrostatic) inflated envelope along with three wind solutions from RM16: a compact ($\dot{M}>\dot{M}_{b}$), extended ($\dot{M}<\dot{M}_{b}$), and critical wind ($\dot{M}=\dot{M}_{b}$), where

[TABLE]

RM16 found their extended wind solutions resembled inflated envelopes so long as max($\alpha)\geq 1$. They solve the subsonic wind is structure between an inner boundary at $T=10^{5.6}\,$K, where the stellar properties are defined by Schaerer & Maeder (1992), and the outer boundary at the sonic point. Integration from the sonic point outwards gives the supersonic wind solution, which truncates when the flow is either subsonic again or reaches the He-bump.

To estimate the critical mass loss rate, $\dot{M}_{b}$, separating compact and extended wind solutions, RM16 use the density, $\rho_{b}$, and temperature, $T_{b}$, where $\alpha=1$ at the iron peak (red circle in Fig. 4). The critical wind solution is approximately tangential to the $\alpha=1$ contour (i.e., $\rho\simeq 4T^{\prime}$) which implies the velocity gradient follows

[TABLE]

using Eq. (10) and (11). With Eq. (13), the supersonic velocity at the iron peak is then

[TABLE]

Since the temperature scale height is small where max($\alpha)\leq 1$, the stellar radius, $R_{*}=r(T=10^{5.6{\rm\,K}})$ or Eq. (29) is taken for $R_{b}$. RM16 found the approximations for the critical mass loss relation $\dot{M}_{b}$ accurately distinguish their WR wind solutions.

The estimate for $\dot{M}_{b}$ is based on the supersonic conditions of the outflow where $c_{g}\lesssim v_{b}\lesssim 300$ km s${}^{-1}\sim(0.1-0.2)v_{\infty}$. The Rosseland mean approximation certainly fails at this point and Doppler effects become important. The neglect of radiative forces besides continuum acceleration implies $\dot{M}_{b}$ is an overestimate. In other words, envelope inflation, as defined in §2.1.2, is likely erased for mass loss rates smaller than $\dot{M}_{b}$ once all acceleration mechanisms are accounted for.

2.3. Validity of Approximations

2.3.1 Radiative Envelope Structure

The dynamic condition for convective instability is (RM16)

[TABLE]

RM16 found little acceleration is necessary to suppress convection in a radiation-dominated envelope since $\Gamma_{c}\simeq 1$ (Eq. 1). Because this study considers solutions where $v^{\prime}_{\rm sp}=0$, the convective instability criterion is satisfied in some region in the subsonic domain for the slowest winds. The minimum mass loss rate to suppress subsonic convection (solving Eq. 21 with 1, 13, and 25) differs from $\dot{M}_{\rm sp}$ by a few percent. Therefore, convection is not important for studies of WR sonic points.

Jiang et al. (2015) find massive hydrogen envelopes with a subsurface iron bump are inherently unstable to radiatively-driven hydrodynamic instabilities (RHIs; Blaes & Socrates 2003). Their StarTop model shows the density inversion from a hydrostatic one-dimensional model initially fragments, collapses, and washes out in multi-dimensional hydrodynamic simulations with radiative transfer. The authors found non-adiabatic mixing length models significantly overestimate the convective flux and find the effects of porosity and vertical oscillations to be important. While turbulent motions can transport kinetic energy, the flux remains small since the motions are comparable to the gas sound speed. The radiative acceleration remains larger than the gravitational acceleration and, as a result, a time- and spatially-averaged density inversion manifests at later times.

RM16 argue the criterion for RHIs is satisfied beneath the iron peak in hydrostatic WR envelopes around $T\sim 10^{5.6}$ K. In stellar oscillation theory, these motions are referred to as ‘strange’ modes (Gautschy & Glatzel, 1990; Blaes & Socrates, 2003) which have been found to excite pulsations in quasi-static helium envelopes models (Glatzel et al., 1993; Glatzel, 1994; Grassitelli et al., 2016). Turbulent energy transport is inefficient once the optical depth per pressure scale height, $\tau_{0}=\kappa\rho H$, is less than $\tau_{b}=c/c_{g}$ (Jiang et al., 2015). Since $H\simeq H_{T}/4$ in a radiation-dominated envelope, we have

[TABLE]

using Eq. (5), (10), (11), (29) and (30). Because $\tau_{0}/\tau_{b}$ is very small for typical WR parameters, the envelopes at the iron bump resemble the StarTop model by Jiang et al. (2015). Therefore, turbulent energy transport is not expected to be important in WR envelopes and the adiabatic estimate for convective transport, Eq. (2), is likely overestimated. This further supports the argument that convection is not energetically relevant around the iron bump in WR stars.

2.3.2 Rosseland Opacity and Microturbulence

Nugis & Lamers (2002) and RM16 argue the Doppler shift from flow acceleration is small in comparison to the iron line widths across the sonic point. They apply the Rosseland mean approximation in their wind calculations so long as the optical depth parameter (Castor et al., 1975) is large,

[TABLE]

where $\sigma_{e}=0.2\,$cm2 g*-1* is the electron scattering opacity of helium. A problem with the authors’ argument is the referenced thermal width, $v_{\rm th}=\sqrt{2}c_{g}$, is for hydrogen and not iron. Using the iron thermal width, $v_{\rm th,Fe}\simeq 7$ km s*-1*, reduces $t_{\rm CAK}\propto\mu^{-1/2}$ by about $7.5$. For weak or low density winds $\dot{M}\lesssim 10^{-5}\,\mathrm{M}_{\odot}\,$yr*-1*, the Rosseland approximation degrades close to the sonic point but not significantly: $\mathrm{min}(t_{\rm CAK})\gtrsim 0.3$, in comparison to O/B stellar atmospheres/winds where $t_{\rm CAK,O}\simeq 10^{-2}$.

Turbulent motions from sub-surface convection (Grassitelli et al., 2016) can broaden photospheric iron lines (Cantiello et al., 2009), increase $\mathrm{min}(t_{\rm CAK})$ and reduce the effects from Doppler motion. As discussed in §2.3.1, WR envelopes are susceptible to RHIs and may not be time-steady on scales smaller than the local scale height. If these motions are acoustic (barring non-linear acoustic effects), the photon mean free path, $\ell_{\rm ph}\sim(\kappa_{\rm Edd}\rho)^{-1}\propto\rho^{-1}$, increases more rapidly than the wavelength, $\lambda\propto c_{s}\propto\rho^{1/6}$ ($\lambda\propto\rho^{-1/2}$ for isothermal gas), from high to low temperatures. Therefore, inhomogeneities from RHIs are bound to become optically thin and become highly compressible. Indeed, Jiang et al. (2015) see strong shock formation in their StarTop models.

Jiang et al. (2015) show the inhomogeneous turbulent gas is radiatively porous in their StarTop models as photons preferentially flow through channels of smaller relative optical depth. The reduction in the effective radiative acceleration from porosity affects the stellar structure, although, weakly. These calculations employ the OPAL opacity tables which assumes the fluid is motionless on scales smaller than the photon mean free path. If WR envelopes harbour small-scale or ‘microturbulent’ motions then the desaturation of optically thick lines may increase the amount of line overlap, the total opacity, and effective radiative acceleration. Microturbulence may be an important factor in determining the structures of outer WR envelopes and winds. To employ the Rosseland mean approximation in the subsonic domain of weak winds, microturbulence must be the dominant contributor to the total line width of iron, $v_{\rm th,Fe}/[\mathrm{min}(t_{\rm CAK})]\sim 23$ km s${}^{-1}\simeq 0.7c_{g}$, such that $t_{\rm CAK}\gtrsim 1$. This is not a restrictive requirement considering the hydrostatic iron line widths are small and subsonic $v_{\rm th,Fe}\simeq 7$ km s${}^{-1}\simeq 0.2c_{g}$.

In §3, I show that an opacity enhancement is necessary for WR winds to have sonic points at the iron bump. Because opacity tables do not include the effects from microturbulent broadening, I use opacity tables with enhanced metallicity as a rough approximation for this effect. While a metallicity enhancement includes opacity contributions from other elements, iron is by far the dominant opacity source at this temperature regime. Implicit in these approximations is the assumption that the turbulent gas can be described using time-steady calculations. This is not justifiable since there are currently no multidimensional simulations of WR outflows.

2.3.3 Enhanced Iron Opacity Tables

Rosseland mean opacity tables are found to underestimate the contributions from iron under solar tachocline conditions ($T\simeq 2\times 10^{6}$ K, $n_{e}\simeq 10^{22-22.5}$ cm*-3*) in direct laboratory measurements by 55-75% (Bailey et al., 2015). Half of the helioseismic discrepancy is resolved if this correction is included. This discrepancy may be eliminated once other elements (e.g., Ni) are included. Corrected regions of stellar instability no longer show a discrepancy for the log $T_{\rm eff}$ vs. log $g$ locations of hybrid slowly pulsating B(SPB)-$\beta$ Cephei pulsators (Moravveji, 2016).

The conditions for wind and inflated envelope formation in WR stellar models are strongly sensitive to the iron opacity. While the iron bump considered here is at cooler temperatures than laboratory conditions, it has several more bound-bound and bound-free transitions available and, so, opacity enhancement corrections are possible.

3. Results

In summary, there are two critical mass loss relations for a WR stellar model to harbour a sonic point, $\dot{M}_{\rm sp}(L_{*})$, and inflated envelope, $\dot{M}_{b}(L_{*})$. These relations partition ($\dot{M},L_{*}$)-space into four regions that correspond to the following responses to the iron bump:

Compact envelope (no wind): $\dot{M}<\dot{M}_{\rm sp}$ and $L_{*}<L_{b}$ 2. 2.

Inflated envelope (no wind): $\dot{M}<\dot{M}_{\rm sp}$, $L_{*}>L_{b}$ 3. 3.

Compact wind: $\dot{M}>\dot{M}_{\rm sp}$ and $\dot{M}>\dot{M}_{b}(L_{*})$ 4. 4.

Extended wind: $\dot{M}>\dot{M}_{\rm sp}$ and $\dot{M}<\dot{M}_{b}(L_{*})$

Note that WR models with $L_{*}\lesssim L_{b}$ have max$(\alpha)<1$ and do not have inflated envelopes as defined in §2.1.2. Therefore, the mass loss rate necessary to erase an inflated envelope is definitively zero: $\dot{M}_{b}(L_{*}<L_{b})=0$.

Figures 5a,b show the critical mass loss relations for metallicity abundances in the Milky Way (MW) and Large Magellanic Cloud (LMC), respectively, and WNE and WC chemical compositions (§2.1.2). Also shown are observations of single hydrogen-free WNE and WC stars (Nugis & Lamers, 2000; Hamann et al., 2006; Sander et al., 2012; Hainich et al., 2014; Tramper et al., 2015, 2016) along with empirical mass loss relations from Tramper et al. (2016), $\dot{M}_{\rm TSK}$, and Yoon (2017), $\dot{M}_{\rm Yoon}$. The rates from Hamann et al. (2006) are adjusted so the clumping factor is $D=10$, instead of 4, to homogenize the data set (Yoon, 2017). Redundant observations between multiple surveys are connected by lines, which provide an impression of the systematic errors. Uncertainties are cited to be of order $\Delta[\dot{M},L_{\rm LMC},L_{\rm MW}]\sim[0.2,0.2,0.4]$ dex. Non-solid markers are WR stars with inferred radii, $R_{20}$, above a factor of three larger than predicted from hydrostatic models without envelope inflation at the same total luminosity. The chemical composition is found to have a relatively small influence. I use the critical mass loss elation from WNE stellar models in the remaining panels in Fig. 5.

4. Discussion

4.1. Are Wolf-Rayet Envelopes Inflated?

Models of MW and LMC WR stars are found to have either sufficiently small luminosities, $L_{*}<L_{b}$, or large mass loss rates, $\dot{M}>\dot{M}_{b}$, to prevent the formation of an extended wind (i.e., envelope inflation) at the iron bump. Therefore, the hydrostatic approximation is not valid for WR models with both envelope inflation and observed rates of mass loss. In other words, hydrostatic WR models with envelope inflation and empirical rates of mass loss are not consistent with one-dimensional hydrodynamic stellar equations. This is a conservative statement considering the neglect of additional radiative forces (e.g., Doppler shifts, line-deshadowing) and definition of ‘inflated envelopes’ used in this study (see §2.1.2). Decreasing the clumping factor from $D=10$ increases the observed mass loss rates and further degrades the hydrostatic approximation.

LMC WRs are observed to have slightly smaller mass loss rates and large inferred radii. These observations could be described by an extended wind with a larger iron opacity bump. This may be achieved with a substantial increase in either the iron abundance or opacity enhancement from microturbulent motions.

G12 propose macroscopic imhomogeneities as an opacity enhancement mechanism to increase the effects of envelope inflation. An opacity enhancement could also drive stronger outflows which would inhibit envelope inflation. Inhomogeneities in multidimensional simulations from Jiang et al. (2015) are found to render a medium porous to radiation which would deflate an inflated WR envelope. The effects of inhomgeneties remain speculative for now and until multidimensional simulations of WR envelopes arise. Without appeals to substantial opacity enhancement, one-dimensional WR models with observed rates of mass loss do not harbour inflated envelopes and, consequently, the radius discrepancy problem remains outstanding.

4.2. Global Minimum Mass Loss Rate

Fig. 5a,b show $\dot{M}_{\rm sp}(L_{*})$ for WNE and WC stellar models with solar and sub-solar metallicities. Less massive models have smaller $L_{*}/M_{*}$ ratios and require larger opacities, densities, and mass loss rates to harbour a sonic point. This is reflected in the rise of $\dot{M}_{\rm sp}(L_{*})$ for smaller $L_{*}\lesssim 10^{5.7}$ $\mathrm{L}_{\odot}$. The $\dot{M}_{\rm sp}(L_{*})$ relation is in agreement with the semi-analytical results by Ro (2017) and Grassitelli et al. (2018) for $L_{*}\lesssim L_{b}\sim 10^{5.7}$ $\mathrm{L}_{\odot}$ or $M_{*}\lesssim M_{b}\sim 19$ $\mathrm{M}_{\odot}$.

Stellar models with $L_{*}\gtrsim L_{b}$ develop inflated envelopes and have sonic points at radii larger than the core (i.e., $R_{\rm sp}/R_{*}\gtrsim 1.4$). Furthermore, the minimum sonic point density for a wind-like solution (i.e., $v^{\prime}_{\rm sp}\geq 0$) increases, rather than decreases, with stellar luminosity. As a result, $\dot{M}_{\rm sp}(L_{*})$ increases for $L_{*}\gtrsim L_{b}$ and a global minimum mass loss rate emerges at roughly min($\dot{M})\simeq 10^{-6}\,\mathrm{M}_{\odot}$ yr*-1*. This increases slightly to min($\dot{M})\simeq 10^{-5.7}\,\mathrm{M}_{\odot}$ yr*-1* for WR models with LMC metallicities which have $L_{b}\sim 10^{6}$ $\mathrm{L}_{\odot}$. The discovery of a global minimum mass loss rate for an iron bump to harbour a sonic point is only found when including the effects from envelope inflation and the conditions for wind-like sonic points.

4.3. Sonic Point Discrepancy

A small fraction of MW WRs in Fig. 5a appear to have insufficient mass loss rates or luminosity to be driven by the iron bump (i.e., $\dot{M}\lesssim\dot{M}_{\rm sp}(L_{*})$). This ‘discrepancy’ is small in comparison to the relatively large observational uncertainties. Therefore, stellar models of observed MW WR stars can harbour sonic points at the iron bump.

If there were a sonic point discrepancy at the iron bump for MW WR stars then a reduced clumping factor from $D_{\rm MW}=10$ to 4 (i.e., increase $\dot{M}_{\rm MW}$ by $\sqrt{10/4}\simeq 1.6$) is one possible resolution (Fig. 5c). Yoon (2017) suggests such a correction from their evolutionary calculations. The discrepancy can also be explained using an opacity table with 50% increase in metallicity (Fig. 5e). In §2.3.2 and §2.3.3, I describe two motivations for opacity enhancement regarding microturbulence and laboratory measurements. Employing clumping factors larger than $D_{\rm MW}=10$ or opacity tables with smaller metallicity/opacities for MW WR stars would widen the discrepancy.

Half of LMC WR stars have $\dot{M}\lesssim\dot{M}_{\rm sp}$ (Fig. 5b). A homogeneous wind (i.e., $D_{\rm LMC}=1$) does not resolve this discrepancy and would be inconsistent with clumped wind observations (Crowther, 2007; Smith, 2014). An empirical clumping factor of $D_{\rm LMC}=4$ requires an opacity enhancement of twice the LMC metallicity $Z\gtrsim 2Z_{\rm LMC}$ (Fig. 5f).

WR models with $\dot{M}=\dot{M}_{\rm sp}$ have peculiar ‘wind’ structures. First, their winds experience zero acceleration at the iron peak (i.e., $v^{\prime}_{\rm sp}=0$) and a maximum velocity that is the gas sound speed, which is very small, $c_{g}\ll v_{\infty}$. In one-dimension, this suggests the helium opacity bump may drive the wind instead. This could further explain the sonic point discrepancy for LMC WR stars since the density, temperatures, and $\dot{M}_{\rm sp}$ are smaller at the helium bump than the iron bump. There are issues with the prospect of steady winds launched from the helium bump. For example, a density inversion would remain between the iron and helium bumps. One-dimensional, time-steady calculations are likely inappropriate for this problem considering how density inversions are dynamically unstable in multidimensional simulations (Jiang et al., 2015). RM16 found wind models most wind models with $\dot{M}>\dot{M}_{\rm sp}({\rm Fe})$ do not reach escape velocities, but rather decelerate and form a density inversion. If such winds do launch from the helium bump then there exists three sonic points in the velocity field; a scenario where multi-dimensional effects are likely important.

4.4. WR Mass Loss Relation

RM16 found their extended and compact wind models are differentiable if max($\alpha)$ is above or below unity, respectively. A critical mass loss relation for extended winds, $\dot{M}_{b}(L_{*})$, where max($\alpha)=1$ was derived as a result. Suppose, instead, $\alpha$ did not exceed some smaller critical value than one due to other constraints set at, for example, the photosphere, infinity, or at a critical point like those in line-driven winds (Castor et al., 1975). A new mass loss relation, $\dot{M}_{\alpha}(L_{*})$, can be calculated in a similar fashion to $\dot{M}_{b}\equiv\dot{M}_{\alpha=1}$ for a different value of max($\alpha)$ at the iron peak. Fig. 5c-f show $\dot{M}_{\alpha}(L_{*})$ for various critical max($\alpha$) and metallicities. Empirical relations for $\dot{M}_{\alpha}$ are available in Appendix B.

The distribution of mass loss rates from WR stars and $\dot{M}_{\alpha}(L_{*})$ are qualitatively similar. Increasing the metallicity by $\Delta Z=0.01$ and reducing the clumping factor, $D$, such that $\dot{M}\gtrsim\dot{M}_{\rm sp}$ leads to a homogeneous range of $0.1\lesssim$ max($\alpha$) $\lesssim 0.5$ that spans both MW and LMC WR distributions (Fig. 5e,f). Smaller max$(\alpha)$ values shift the $\dot{M}_{\alpha}(L_{*})$ relation towards smaller luminosities (i.e., $L_{\alpha}<L_{b}$). WC stars typically have larger $\dot{M}$ than WNE stars and, so, are better fit with $\dot{M}_{\alpha}(L_{*})$ for smaller max$(\alpha)$.

The $\dot{M}_{\alpha}(L_{*})$ relation increases rapidly with $L_{*}$ for less luminous stellar models since their Eddington ratios are closer to unity and supersonic velocities, $v^{2}\simeq v_{k}^{2}(\Gamma_{r}-1)$, are smaller at the iron peak (see Eq. 24). The growth is more gradual for higher $L_{*}$ since $\Gamma_{r}-1\sim 1$. These dependencies lead to a steep and shallow component to $\dot{M}_{\alpha}(L_{*})$, respectively. The minimum mass loss relation to harbour a wind-like sonic point, $\dot{M}_{\rm sp}$, truncates the steep component for less luminous stellar models. As a result, less luminous stars with larger $\dot{M}$ have shallower overall gradients in $d\dot{M}/dL_{*}$ when compared to more luminous stellar models with smaller $\dot{M}$. This is a noted distinction between WC and WNE populations (Vink & de Koter, 2005; Hainich et al., 2014; Tramper et al., 2016; Yoon, 2017). The steep component reduces for larger $\dot{M}_{\rm sp}$ or smaller metallicities (e.g., Fig. 5d and f).

4.4.1 $\beta$-law Wind Models

The densities of supersonic winds are certainly small enough to expect other radiative forces not included here. For example, line-deshadowing is responsible for driving optically thin winds from O/B stars which follow a $\beta$-law velocity structure (Castor et al., 1975). Gräfener & Vink (2013) and Gräfener et al. (2017) use a $\beta$-law as a scaffold between the outer optically thin and inner optically thick wind to calculate the radiative forces. That is, the radiative forces are assumed to drive a $\beta$-law velocity structure. Their outer boundary condition is the ratio of wind speeds at infinity and $R_{20}$. Although, the physical meaning of $R_{20}$ is unclear if it is neither the sonic point at the iron bump, $R_{\rm sp}$, nor the core radius, $R_{*}$. Their wind models generate a mass loss relation that matches the galactic WNE distribution for $D=10$. The disagreement between their relation with the LMC WR distribution may come from the non-existence of a sonic point at the iron bump ($\dot{M}<\dot{M}_{\rm sp}$). As the authors suggest, this may further indicate the sonic point resides at the helium opacity bumps. Another resolution is to reduce the clumping factor and an enhance the iron opacity bump.

A prediction from Gräfener & Vink (2013), Gräfener et al. (2017), and Nakauchi & Saio (2018) is a nearly constant ratio of radiation and gas pressure $p_{r}/p_{g}\sim 100-160$ in the supersonic wind. Fig. 5c,d suggest galactic WNE winds are best fit with $\dot{M}_{\alpha}(L_{*})$ for max($\alpha)\sim 0.2$. This sets a ratio of radiation to gas pressure,

[TABLE]

at the iron peak using our empirical WR relations (Appendix A) and Eq. (10). This ratio is in agreement with the mentioned authors’ findings despite the different approaches to our wind model calculations. Envelope inflation, as defined in §2.1.2, is not consistent with their models considering $p_{r}/p_{g}\gtrsim 700$ for max$(\alpha)\gtrsim 1$. Our models suggests max$(\alpha)$ and $p_{r}/p_{g}$ in WC winds is about one-half of that in WNE winds.

5. Summary

This investigation explores several aspects regarding the structure of WR stellar models at the iron bump. First, a radiation-dominated, radiative, near-Eddington and hydrostatic envelope must develop either an atmosphere or density inversion at a prominent opacity bump (§2.1). The envelope structure is primarily determined by the Eddington contour (§2.1.1), which depends on the $L_{*}/M_{*}$ ratio and the chemical composition. Models with larger $L_{*}/M_{*}$ or iron abundance have Eddington contours at smaller densities which increases both the temperature scale height, $H_{T}\propto\rho^{-1}$, and envelope size. The radius discrepancy problem is the original motivation for this investigation and, so, a heuristic and conservative definition for envelope inflation is adopted: $\alpha\equiv H_{T}/R_{*}>1$, where $R_{*}$ is the radius at $T=10^{5.6}$ K.

Second, this investigation compares the observed distribution of WR mass loss rates to the critical rate, $\dot{M}_{b}$, to harbour extended winds (i.e., inflated envelopes with time-steady, transonic flow) assuming typical wind parameters and opacity tables. I show extended winds do not manifest since the observed WR distribution have either insufficient luminosities or strong winds. This suggests previous hydrostatic WR stellar models with both envelope inflation and mass loss are not consistent with the one-dimensional hydrodynamic stellar structure equations. In other words, the hydrostatic approximation is invalid near the iron bump in WR stellar models. This conclusion is conservative considering the effects from Doppler shifts and microturbulent line broadening, which are expected to increase the radiative acceleration, are not included.

Envelope inflation is a problematic resolution to the WR radius discrepancy problem for various other reasons. For example, the radius discrepancy problem would not be resolved but, rather, inverted for MW WR stars (i.e., why are some MW WR radii too small?). Another example is that WR wind models suggest the structures have constant radiation-to-gas pressure ratios of $p_{r}/p_{g}\simeq 100-160$ (Gräfener & Vink, 2013; Gräfener et al., 2017; Nakauchi & Saio, 2018). These values correspond to $\alpha\simeq 0.2$ and temperature scale heights (see §4.4.1) that are too small for envelope inflation to manifest using the definition in §2.1.2. WR stellar models with higher metallicities develop inflated envelopes at smaller luminosities. Fig. 5 suggests more than twice the host environment metallicity is required to resolve the radius discrepancy problem with this explanation.

Third, Ro (2017) and Grassitelli et al. (2018) calculate the minimum mass loss rate for a WR stellar model to harbour a sonic point on the hot side of the iron bump (i.e., $\kappa_{T}<0$; Nugis & Lamers 2002). Respectively, the authors approximate the sonic point radius with either the mass-radius relation from Schaerer & Maeder (1992) or a constant value of 1 $\mathrm{R}_{\odot}$. Here, I show these approximations are valid for less massive stellar models without envelope inflation (i.e., $M_{*}\lesssim 19$ $\mathrm{M}_{\odot}$). In more massive stellar models, the minimum sonic point density for wind-like solutions is larger than suggested by conditions at the iron peak once the entire sonic point condition is included. The sonic point radius is also larger than previously approximated due to the effects of envelope inflation. This combination is found to increase, rather than decrease, the minimum mass loss relation with luminosity for hydrogen-free WR stars with $M_{*}\gtrsim 19$ $\mathrm{M}_{\odot}$ or $L_{*}\gtrsim 10^{5.7}$ $\mathrm{L}_{\odot}$ with solar metallicity. Likewise, $M_{*}\gtrsim 27$ $\mathrm{M}_{\odot}$ or $L_{*}\gtrsim 10^{5.9}$ $\mathrm{L}_{\odot}$ for models with LMC metallicities. As a result, hydrogen-free WR stars in the MW and LMC have global minimum mass loss rates of $10^{-6}$ and $10^{-5.7}$ $\mathrm{M}_{\odot}$ yr*-1*, respectively, if they harbour sonic points at the iron bump.

Fourth, half of early-type WR stars in the LMC have sufficiently small $\dot{M}$ such that their sonic points do not reside at the iron bump. The discrepancy vanishes given twice the LMC metallicity and a reduced clumping factor from $D=10$ to 4. Without a reduction in the clumping factor, super-solar metallicities are necessary. In this paper, I suggest microturbulence and laboratory corrections may enhance the opacity tables and resolve the sonic point discrepancy.

Indeed, early-type WR stars with mass loss rates above the minimum, $\dot{M}_{\rm sp}$, does not guarantee the winds escape via the iron bump. RM16 found nearly all of their wind solutions did not escape due to the steep potential; although, the effects from Doppler shifts were not included. Another scenario is that a subset of WR winds escape via the helium bump instead. That is, WR mass loss is not driven by radiation on UV lines but rather by the continuum field on the helium bump. The viability of this scenario requires an understanding of the deeper structure around the iron bump which may not be describable with classic mixing length theory, existing opacity tables, or one-dimensional time-steady calculations.

Lastly, the critical mass loss relation separating extended and compact winds is found to trace the observed WR distribution when smaller values of ${\rm max}(\alpha)=1$ are considered (§4.4). Nominal values of ${\rm max}(\alpha)\simeq 0.1-0.2$ correspond to $p_{r}/p_{g}\simeq 70-145$, which is the typical range of values argued by Gräfener & Vink (2013), Gräfener et al. (2017), and Nakauchi & Saio (2018). If the critical relation is accurate then it suggests the mass loss rate declines faster than a single power-law at late stages of evolution. This may have interesting implications for both the star and its circumstellar environment before core-collapse.

SR thanks Chris Matzner, Eliot Quataert, and Yan-Fei Jiang for useful conversations. SR appreciates the referee for their comments and suggestions. This research was funded by the Gordon and Betty Moore Foundation through Grant GBMF5076.

Appendix A

Empirical Fits to WR MESA Models

I use the test suite model make_massive_with_ uniform_composition from MESA-r9793 to construct models of pure helium and mixture of He/C/O. I assume these are representative of WNE and WC interiors, respectively. See §2.1.1 for the chemical mass fractions used. The evolution is stopped after the He core mass fraction is reduced by $2\%$. The luminosities and radii at $10^{5.6}$ K are shown in Figure 6 for stellar masses between 5 and 50 $\mathrm{M}_{\odot}$. The fitted empirical relations are

[TABLE]

where $\ell\equiv{\rm log_{10}}(L_{*}/\mathrm{L}_{\odot})$, $\mathcal{R}\equiv{\rm log_{10}}(R_{*}/\mathrm{R}_{\odot})$, and $m\equiv{\rm log_{10}}(M_{*}/\mathrm{M}_{\odot})$. See Fig. 6 for a comparison with Schaerer & Maeder (1992). Linear fits to the luminosities are shown for convenience:

[TABLE]

An inlist is available at the end of Appendix B that has been verified for MESA-r10108. The mixing length in our $23\,\mathrm{M}_{\odot}$ inlist is tuned to reproduce the G12 results.

Appendix B

Critical Mass Loss Relations

The minimum mass loss rate, $\dot{m}\equiv\mathrm{log_{10}}(\dot{M}/(\mathrm{M}_{\odot}\,\mathrm{yr}^{-1}))$, for an iron bump to harbour a sonic point in a hydrogen-free WR stellar model (Eq. 21) is approximately:

[TABLE]

with a global minimum at $\ell_{\rm sp}(Z)\equiv 6.2-20Z$, where $0.01\leq Z\leq 0.03$ is the metallicity.

An approximation for the WR mass loss relation $\dot{M}_{\alpha}$ (Eq. 22) follows a double logarithm of $L_{*}$,

[TABLE]

where

[TABLE]

for metallicities at $Z=[0.01,0.02,0.03]$ and $0.05\leq\alpha\leq 1$. Note that Eq. (32) is undefined for $\dot{m}_{\alpha}<\dot{m}_{\rm sp}$ since there is no sonic point. For $\alpha=1$, the maximum mass loss rate to harbour extended winds (i.e., inflated envelopes), $\dot{M}_{b}\equiv\dot{M}_{\alpha=1}$, asymptotically approaches zero at the minimum luminosity for envelope inflation in hydrostatic models $\ell_{b}$. The values for $\ell_{b}$ are fitted parameters from Eq. (32) and (22) which overestimate $L_{b}$ by $\sim 0.15$ dex.

⬇ &star_job

relax_initial_to_xaccrete = .true.

num_steps_to_relax_bomposition = 50

set_HELM_ion_neutral_blends = .true.

max_logRho_neutral_HELM = -20

relax_to_this_tau_factor = 1d0

dlogtau_factor = 0.1

relax_initial_tau_factor = .true.

set_initial_dt = .true.

years_for_initial_dt = 1d-8

/

&controls

max_number_backups = 150

max_number_retries = 600

initial_mass = 23d0

Zbase = 0.02d0

varcontrol_target = 1d-3

dX_nuc_drop_limit = 2d-2 ! 5d-3

delta_lgT_bntr_limit = 0.01

delta_lgT_bntr_hard_limit = 0.03

delta_lgRho_bntr_limit = 0.05

delta_lgRho_bntr_hard_limit = 0.15

cool_wind_RGB_scheme = ”

cool_wind_AGB_scheme = ”

accrete_same_as_surface = .false.

accrete_given_mass_fractions = .true.

! Pure He and stopping condition

num_accretion_species = 1

accretion_species_id(1) = ’he4’

accretion_species_xa(1) = 1d0

xa_bentral_lower_limit_species(1)=’he4’

xa_bentral_lower_limit(1) = 98d-2

! MLT++ off:

okay_to_reduce_gradT_excess = .false.

mixing_length_alpha = 0.325!tuned to G12

mesh_Pgas_div_P_exponent = 0.3

max_allowed_nz = 20000

/