11institutetext: Astronomical Institute Anton Pannekoek, University of
Amsterdam, 1098 XH Amsterdam, The Netherlands22institutetext: Dipartimento di
Fisica ‘Enrico Fermi,’ University of Pisa, I-56127 Pisa, Italy33institutetext:
TAPIR, Walter Burke Institute for Theoretical Physics, MC 350-17, California
Institute of Technology, Pasadena, CA 91125, USA44institutetext: Kavli
Institute for Theoretical Physics, University of California, Santa Barbara, CA
93106, USA

# A Systematic Survey of the Effects of Wind Mass Loss Algorithms on the
Evolution of Single Massive Stars

M. Renzo 11223344 C. D. Ott 33 S. N. Shore 22 S. E. de Mink 11

Mass loss processes are a key uncertainty in the evolution of massive stars.
They determine the amount of mass and angular momentum retained by the star,
thus influencing its evolution and presupernova structure. Because of the high
complexity of the physical processes driving mass loss, stellar evolution
calculations must employ parametric algorithms, and usually only include wind
mass loss. We carry out an extensive parameter study of wind mass loss and its
effects on massive star evolution using the open-source stellar evolution code
MESA. We provide a systematic comparison of wind mass loss algorithms for
solar-metallicity, nonrotating, single stars in the initial mass range of
$15\,M_{\odot}$ to $35M_{\odot}$. We consider combinations drawn from two hot
phase (i.e. roughly the main sequence) algorithms, three cool phase (i.e.
post-main-sequence) algorithms, and two Wolf-Rayet mass loss algorithms. We
discuss separately the effects of mass loss in each of these phases. In
addition, we consider linear wind efficiency scale factors of $1$, $0.33$, and
$0.1$ to account for suggested reductions in mass loss rates due to wind
inhomogeneities. We find that the initial to final mass mapping for each zero-
age main-sequence (ZAMS) mass has a $\sim 50\%$ uncertainty if all algorithm
combinations and wind efficiencies are considered. The ad-hoc efficiency scale
factor dominates this uncertainty. While the final total mass and internal
structure of our models vary tremendously with mass loss treatment, final
luminosity and effective temperature are much less sensitive for stars with
ZAMS mass $\lesssim 30\,M_{\odot}$. This indicates that uncertainty in wind
mass loss does not negatively affect estimates of the ZAMS mass of most
single-star supernova progenitors from pre-explosion observations. Our results
furthermore show that the internal structure of presupernova stars is
sensitive to variations in both main sequence and post main-sequence mass
loss. The compactness parameter $\xi\propto\mathcal{M}/R(\mathcal{M})$ has
been identified as a proxy for the “explodability” of a given presupernova
model. We find that $\xi$ varies by as much as $30\%$ for models of the same
ZAMS mass evolved with different wind efficiencies and mass loss algorithm
combinations. This suggests that the details of the mass loss treatment might
bias the outcome of detailed core-collapse supernova calculations and the
predictions for neutron star and black hole formation.

###### Key Words.:

stars: evolution – massive – mass loss – winds – supernovae: general

††offprints<EMAIL_ADDRESS>
##  1 Introduction

Mass loss is a key phenomenon for the co-evolution of massive stars ($M\gtrsim
8\,M_{\odot}$) and their environment, yet it is poorly understood. It plays an
important role throughout the stellar evolution and it may have a deciding
influence on the outcome of core collapse. Mass loss is responsible for a
large part of the chemical enrichment of the interstellar medium, and its
momentum input can trigger star formation, but it can also sweep away gas from
stellar clusters, preventing further star formation.

In the standard picture of single massive star evolution, mass loss influences
the duration of different evolutionary phases (e.g., Meynet et al., 2015),
especially the amount of time spent on the red supergiant (RSG) branch. It has
been suggested to be important for the solution of the so-called “red
supergiant problem” (Smartt et al., 2009; Smith et al., 2011), i.e. the
discrepancy between the observed maximum mass for type IIP supernovae (SNe)
and the theoretical predictions for the core collapse of RSG stars. This
discrepancy indicates an incomplete theoretical understanding of the evolution
and explosion of massive stars, especially in the mass range
$\sim$$16-30\,M_{\odot}$.

Mass loss also plays an essential role in the formation of Wolf-Rayet (WR)
stars. Two competing scenarios exist for the removal of (most of) the
hydrogen-rich envelope of stars: the single star picture (so-called “Conti
scenario”, Conti, 1975; Maeder & Conti, 1994; Lamers, 2013) and the binary
formation channel (e.g., Woosley et al., 1995; Wellstein & Langer, 1999; Smith
& Tombleson, 2015; Shara et al., 2017). Understanding mass loss phenomena is
necessary to discriminate between these two scenarios. Mass loss is invoked to
explain the variety of core-collapse SNe (e.g., Eldridge & Tout, 2004; Smith
et al., 2011; Georgy, 2012; Groh et al., 2013; Smith, 2014), because it can
modify the surface composition of the pre-SN star, possibly removing the
hydrogen-rich envelope (and perhaps some of the helium-rich shell) and leading
to type Ib/c SNe. Observations of SNe IIn and superluminous SNe may perhaps be
explained by the strong interaction between the SN shock and shells of
material ejected from the star in late stages of its evolution, (e.g., Smith
et al., 2011; Smith, 2014; Shiode & Quataert, 2014). It has also been proposed
that the coupling of mass loss and rotation might prevent pair instability SNe
for very massive, low metallicity stars (Ekström et al., 2008; Woosley, 2016).

Finally, mass loss plays a key role in shaping the (hydrostatic) internal
structure of massive stars at the pre-SN stage, which can influence the
expected SN outcome (see, e.g., Belczynski et al., 2010; O’Connor & Ott, 2011,
2013; Ugliano et al., 2012; Sukhbold et al., 2016): will the star successfully
explode and leave a neutron star (NS) remnant? Will the explosion be highly
asymmetric or weak, leading to fallback accretion and black hole (BH)
formation? Or will the explosion fail completely, leaving a BH with little
(Lovegrove & Woosley, 2013) or no electromagnetic counterpart?

Depending on the amount and geometry of the ejecta, which are governed by the
mass loss during the stellar lifetime and the SN energetics, also the SN kick
can vary (Zwicky, 1957; Blaauw, 1961; Boersma, 1961; Janka, 2013, 2016). The
kick can change the post-explosion orbital parameters if the star is in a
binary system and, therefore, our incomplete understanding of massive star
mass loss affects the predicted populations of sources for gravitational wave
astronomy (Abbott et al., 2016).

There exist three main channels of mass loss in the evolution of massive
stars:

  1. 1.

Steady-state winds. These are radiatively driven in hot stars (i.e., on the
main sequence; see Puls et al. 2008 for a review). In the supergiant phase,
the driving mechanism is uncertain. Winds could be radiatively driven via
lines or dust, or by other mechanisms, e.g., wave energy deposition (see,
e.g., Bennett 2010 for a brief review).

  2. 2.

Impulsive, pulsational and/or eruptive mass loss, e.g. disk shedding at
criticial rotation or giant eruptions such as those of Luminous Blue Variables
(LBVs) (e.g., Puls et al., 2008; Smith, 2014; Puls et al., 2015);

  3. 3.

Roche lobe overflow (RLOF) and possibly common envelope ejection in binary
systems, which can result in mass loss from the donor star and also mass loss
from the system as a whole in non-conservative cases.

Which of these processes dominates in terms of the total mass shed has been a
matter of debate in the literature (see Smith, 2014, and references therein).

Mass loss is an intrinsically dynamical phenomenon which involves bulk
acceleration of matter for escaping the star. Because of the dynamical nature
of mass loss, it is difficult to include in stellar evolution codes: most
simulations focus on single stars, and are carried out with hydrostatic codes
which cannot account for dynamical or impulsive events in a physical way. Even
hydrodynamical codes do not permit a self-consistent development of impulsive
events, since the physical processes triggering them are currently poorly
understood or even unknown. Impulsive outbursts of mass loss can, however, be
included using physically plausible prescriptions (e.g., Morozova et al.,
2015).

Most massive stars are found in binary systems (e.g., Mason et al., 2009; Sana
& Evans, 2011; Sana et al., 2012; Kiminki & Kobulnicky, 2012; Chini et al.,
2012; Kobulnicky et al., 2014; Almeida et al., 2016), where mass loss also
determines the angular momentum losses, thus the orbital evolution, and
ultimately the binary evolution path and its end point (merger, disruption of
the binary system, double compact object binary, etc.).

In this paper, our focus is on steady, radiatively-driven wind mass loss (Lucy
& Solomon, 1970) and our goal is to understand how different treatments of
this process affect massive star evolution and pre-SN structure. We consider
single massive stars and do not address the problem of mass loss in binaries
directly.

Most evolutionary calculations of massive stars only include wind mass loss,
which can be treated in the steady state approximation, although, strictly
speaking, a wind is dynamical as well.

The wind is in fact driven by a radiative acceleration which formally enters
into the momentum equation for the stellar plasma (Castor et al., 1975), the
stellar structure responds secularly, since wind mass loss rates are low and
change only slowly compared to impulsive mass loss events. At solar
metallicity, the wind mass loss rate is
$10^{-10}\,M_{\odot}\,\mathrm{yr^{-1}}\lesssim\dot{M}_{\mathrm{wind}}\lesssim
10^{-5}\,M_{\odot}\,\mathrm{yr^{-1}}$, while RLOF and LBV eruptions yield
$10^{-6}\,M_{\odot}\,\mathrm{yr^{-1}}\lesssim\dot{M}_{\mathrm{RLOF}}\lesssim\dot{M}_{LBV}\lesssim
10^{-2}\,M_{\odot}\,\mathrm{yr^{-1}}$, (e.g., de Jager et al., 1988; Vink et
al., 2001; van Loon et al., 2005; Smith & Owocki, 2006; Crowther, 2007; Puls
et al., 2008; Langer, 2012; Smith, 2014). These mass loss rates correspond to
timescales
$\tau_{\mathrm{wind}}=M/\dot{M}_{\mathrm{wind}}\gg\tau_{\mathrm{RLOF}},\tau_{\mathrm{LBV}}$.
Therefore, wind mass loss can be included in stellar evolution calculations
using parametric algorithms111These are often called “recipes” in jargon.
However, we prefer the term “algorithm” (Fibonacci, 1202) because it
underlines that these are mathematical representations of physical phenomena
relying on specific sets of assumptions..

Stellar winds are radiatively driven by the interaction of photons with
metallic ions (line-driven mass loss) or dust grains (dust-driven mass loss).
Therefore, they depend on the opacity and thus chemical composition,
ionization state, and density stratification, so indirectly also on the
equation of state of the outermost stellar layers. Metals effectively provide
all the opacity in stellar atmospheres because of their large number of lines.
Photons come out of the photosphere with well defined direction, interact with
metallic atoms/ions via bound-bound processes (absorption and line scattering)
and cede their momentum to the atoms/ions. When these de-excite, they emit
photons isotropically. The momentum of de-excitation photons averages out, and
the result is a net gain of momentum in the direction of the initial photon
(see, e.g., Puls et al., 2008).

Note that the radiation field in the stellar atmosphere, i.e. above the
photosphere, is _not_ isotropic. If there were not a net radial flux of
photons, also the momentum of the incoming photons would average out. In a
nutshell, the incident photons push metals outward, and metals drag hydrogen
and helium through collisional Coulomb coupling (see Puls et al., 2008; Vink,
2015, and references therein).

This simple theoretical picture of line-driven stellar winds is complicated by
two phenomena: the high nonlinearity of the driving mechanism, and the
possible presence of inhomogeneities (so called “clumpiness”) in the stellar
atmosphere. The high nonlinearity arises because the outflow of mass is driven
by the rate of interaction of photons with metals, but this in turn depends on
the local opacity, and therefore on the outflow properties, such as density,
and velocity (which can Doppler shift the lines), see Lamers & Cassinelli
(1999); Puls et al. (2008) and references therein. The presence of
inhomogeneities in the stellar atmosphere is both theoretically expected (see,
e.g., Owocki & Rybicki, 1984; Owocki et al., 1988; Feldmeier, 1995; Owocki &
Puls, 1999; Dessart & Owocki, 2005; Puls et al., 2008; Smith, 2014) and
observed comparing diagnostic spectral lines sensitive to the density $\rho$
(e.g., lines with P Cygni profiles) and $\rho^{2}$ (e.g., recombination lines,
such as H$\alpha$) in the same stellar wind (see e.g. Fullerton et al., 2006;
Bouret et al., 2005; Evans et al., 2004). This comparison shows that the
averaged $\langle\rho^{2}\rangle>\langle\rho\rangle^{2}$ (for reviews see Puls
et al. 2008 and Smith 2014). Therefore, the presence of over-dense clumps in
the wind causes an overestimation of the density inferred from observed
spectral features sensitive to $\rho^{2}$. This is not taken into account in
most wind mass loss algorithms in the literature. The overestimation of the
density directly results in an overestimated mass loss rate. Work by Crowther
et al. (2002); Hillier et al. (2003); Bouret et al. (2005); Fullerton et al.
(2006); Puls et al. (2008); Smith (2014); Abbott et al. (2016), suggests that
the algorithms used in stellar evolution calculations may yield mass loss
rates that are a factor of 2 to 10 too high. Puls et al. (2008) and Smith
(2014) suggest a factor of 3 as the most realistic overestimate. We refer the
interested reader to Puls et al. (2008), Smith (2014), and Renzo (2015) for
more details.

To date, there has been no systematic comparison of the various wind mass loss
algorithms, with varying corrections for clumpiness, and their combinations
and effects throughout the evolution and on the final structure of massive
stars. However, Eldridge & Tout (2004) compared different combinations of
semi-empirical mass loss rates to find the threshold in mass between type Ibc
and type II SNe. Yoon et al. (2010) discussed consequences of the revision of
mass loss rates because of the clumpiness during the WR phase.

In this study, we employ the open-source stellar evolution code MESA (Paxton
et al., 2011, 2013, 2015) to compare various combinations of mass loss
algorithms, using different efficiency factors to account for the
inhomogeneities in the wind (albeit in an _ad-hoc_ fashion). Our aim is to
understand the systematics of massive star evolution and pre-SN structure
caused by variations in the treatment of wind mass loss, focusing on the
differences in the evolution and pre-SN structure (effective temperature,
total mass, core masses, and interior structure).

The remainder of this paper is structured as follows. In Sec. 2, we discuss
some general aspects of the implementation of wind mass loss in stellar
evolution codes and give a very brief overview of the physical bases of the
wind mass loss algorithms we compare. A longer review of these, including the
limitations and the formulae implemented in MESA, can be found in Appendix A.
The more technical points not relevant to the physics of stellar winds are
discussed in Appendix B, with the explicit aim of making our result
reproducible. We compare our models when their final pre-SN mass is determined
in Sec. 3.1. We discuss the impact of winds during the hot, cool, and WR
phases (if reached) separately in Sec. 3.2, 3.3, and 3.4, respectively. We
compare a subset of our models at oxygen depletion in Sec. 3.5. The evolution
and pre-SN structure of the core is also sensitive to the mass loss history of
the stellar model, including the early mass loss during the main sequence, as
we show in Sec. 3.6. In Sec. 3.7, we discuss the evolution of a subset of our
models from oxygen depletion to the onset of core collapse. We discuss the
implications of uncertainties in wind mass loss and potential observational
constraints in Sec. 4, before concluding in Sec. 5.

##  2 Methods

###  2.1 Overview of the mass loss algorithms

Stellar evolution codes do not explicitly compute the acceleration of the gas
unbound in mass loss processes. The usual approach for including wind mass
loss is to use parametric algorithms prescribing a mass loss rate averaged
over each timestep. The time averaging is needed to compute each timestep with
a constant mass loss rate. Homogeneity of the wind is implicit in the standard
formalism, which is known to be a poor approximation and could cause a
significant overestimate of the mass loss rate. Stellar wind algorithms are
either parametric fits to observed mass loss rates, or theoretically derived
models with free parameters chosen empirically or heuristically. Each
algorithm gives a formula for the mass loss rate as a function of some
quantities characterizing the star
$\dot{M}\equiv\dot{M}(L,T_{\mathrm{eff}},Z,\dots)$. The precise set of
variables assumed to be independent varies between the algorithms. Most
algorithms do not include an explicit metallicity dependence, because they
either assume a specific chemical composition of the stellar atmosphere or are
based on observed samples with a specific metallicity $Z$. It is common
practice to impose a smooth scaling with $Z$ such as

$\dot{M}\propto Z^{a}\ \ ,$ (1)

with $a\simeq 0.5$ (e.g., Vink et al., 2000; Woosley et al., 2002). Eq. 1 is
in reasonable agreement with more sophisticated mass loss rate determinations
(see, e.g., Vink et al., 2001), but deviations should be expected both at very
low $Z$ (because of the lack of metal lines to drive the wind) and very high
$Z$ (because of line saturation preventing further driving of the wind). In
this study, we only consider solar metallicity. Note also that stellar
evolution codes usually neglect the errors on the coefficients of the
algorithms obtained as fits to observations.

Since mass loss rates have large uncertainties, it is common practice to
employ a linear efficiency factor $\eta$ to rescale rates to account for
various physical uncertainties. For example, $\eta\lesssim 1$ can be used to
account for wind clumpiness (e.g., Hamann & Koesterke, 1998; Woosley et al.,
2002; Woosley & Heger, 2007). Mass loss channels other than winds exist (e.g.,
LBV eruptions, binary interactions, etc.) and some authors (e.g., Dessart et
al., 2013; Meynet et al., 2015) use $\eta>1$ to explore the effects of
increased total mass loss. However, an averaged steady wind approximation
might not properly capture the readjustment of the stellar structure to non-
wind mass loss events, which may be sensitive to mass loss timing, and/or the
physical process triggering them, and may not be radiatively driven.

Most wind mass loss algorithms are tailored to a specific evolutionary stage.
To carry out a simulation of the entire evolution of the star, several mass
loss algorithms are commonly combined using computational definitions of the
evolutionary phases. This may introduce somewhat arbitrary switching points in
the evolution. Below, we list the physical basis and the abbreviations for the
two hot phase, three cool phase, and two WR phase wind mass loss algorithms
that we combine and compare here. We define each phase of the evolution in
Sec. 2.2.

Vink et al. (2000, 2001) (V): This wind scheme is a theoretical algorithm
obtained with numerical simulation of the line-driven process. It explicitly
includes the metallicity dependence and applies to OB stars during their hot
evolutionary phase. It also includes a detailed treatment of the so-called
“bistability jump”. This corresponds to a non-monotonic behavior of the mass
loss rate $\dot{M}$ as a function of effective temperature $T_{\mathrm{eff}}$
in certain temperature ranges (e.g., $T_{\mathrm{eff}}^{\mathrm{jump}}\simeq
25\,000\,\mathrm{K}$) because of the recombination of certain ions, which
provides more lines in the spectral domain relevant to drive the wind, cf.
Appendix A.1.

Kudritzki et al. (1989) (K): This analytical mass loss rate is obtained using
the Castor et al. (1975) model for line-driven acceleration. It assumes that
the wind is stationary, isothermal, spherically symmetric, and without
viscosity and heat conduction. The analytical solution is obtained assuming a
velocity structure $v\equiv v(r)$ as a function of the radius of the wind, and
solving self-consistently for density and radiative acceleration, cf. Appendix
A.2.

de Jager et al. (1988) (dJ): This empirical mass loss rate describes the
“averaged statistical behavior” of stars (excluding WR and Be stars) in the HR
diagram. It is commonly used for the cool (giant) phase of the evolution of
massive stars, cf. also Appendix A.3.

Nieuwenhuijzen & de Jager (1990) (NJ): This algorithm is also an empirical
mass loss rate drawn from the same sample of stars used by de Jager et al.
(1988). The two algorithms differ in the physical quantities the mass loss
rate is assumed to depend on: NJ used pre-computed stellar models to add a
dependence on the total mass ($\dot{M}\equiv\dot{M}(M)$), which is not a
directly observable quantity for single stars. It is also usually adopted for
the cool phase of stellar evolution. See Appendix A.4 for more details.

van Loon et al. (2005) (vL): This empirical mass loss rate is derived from a
sample of oxygen rich asymptotic giant branch (AGB) and red supergiants (RSG)
stars in the Large Magellanic Cloud. It assumes a dust-driven wind, i.e. mass
loss is driven by photons impinging on dust grains instead of metallic ions.
Note, however, that the presence of dust and its role as a wind driving agent
in supergiant stars is still debated in the literature (van Loon et al., 2005;
Ferrarotti & Gail, 2006), cf. Appendix A.5.

Nugis & Lamers (2000) (NL): This empirical mass loss rate is for WR stars.
Fitting the data for two populations of WR stars (one of known distance and
one for which they carry out a distance determination), they provide an
algorithm, which depends strongly on the surface chemical composition of the
star, cf. Appendix A.6.

Hamann et al. (1982, 1995) (H): This is a theoretical mass loss rate for WR
stars. It is derived assuming a spherically symmetric, homogeneous, and
stationary wind. They avoid solving for the dynamics with a complicated
radiative acceleration term by imposing a velocity structure $v\equiv v(r)$.
In this way, they are able to produce synthetic spectra, which they then fit
to observed WR stars to infer $\dot{M}$. Hamann & Koesterke (1998) suggest to
reduce the mass loss rate by a factor between 2 and 3 to account for wind
clumpiness, cf. Appendix A.7.

We report in Tab. 1 an approximate scaling of the mass loss rate with the
luminosity ($L$), which quantifies the amount of photons available to drive
the wind (neglecting the frequency dependence of the line transitions), and
the effective temperature ($T_{\mathrm{eff}}$), which can be considered as a
rough parametrization of the ionization state at the base of the wind and
therefore the opacity. We strongly recommend to consult Tab. 9 for the scaling
of $\dot{M}$ with physical stellar quantities for each of these algorithms for
anything beyond simple order of magnitude estimates.

Table 1: Approximate scaling of the mass loss rate with luminosity and effective temperature, $\log_{10}(\dot{M})\propto\alpha\log_{10}(L)+\beta\log_{10}(T_{\mathrm{eff}})$, predicted by the considered mass loss algorithms. Appendix A provides the full functional form of the algorithms implemented in MESA and Tab. 9 lists the scaling with all physical quantities.  | ID | $\alpha$ | $\beta$
---|---|---|---
Hot | V | $2.2$ | $1.0$
K | $1.2$ | $0.6$
Cool | dJ | $1.8$ | $-1.7$
NJ | $1.6$ | $-1.6$
vL | $1.0$ | $-6.3$
WR | NL | $1.3$ | $0.0$
H ($L>4.5L_{\odot}$) | $1.5$ | $0.0$
H ($L\leq 4.5L_{\odot}$) | $6.8$ | $0.0$

###  2.2 Combination of mass loss algorithms

Table 2: Combinations of wind mass loss algorithms employed in this study. The temperature threshold separating the hot phase and the cool phase is $T_{\mathrm{th}}=15\,000$ (10 000) K when using the Kudritzki _et al._ (Vink _et al._) algorithm. The WR phase is defined using the surface (outermost computational cell) hydrogen mass fraction $X_{s}$, without constraints on $T_{\mathrm{eff}}$. In the text, we do not mention the WR phase algorithm if the model discussed does not enter this phase. For a description of the algorithms, see Sec. 2.1 and the appendices listed in this table. We discuss the definition of each evolutionary phase in Sec. 2.2.  ID | Hot phase | Cool phase | WR phase
---|---|---|---
$T_{\mathrm{eff}}\gtrsim T_{\mathrm{th}}$ | $T_{\mathrm{eff}}\lesssim T_{\mathrm{th}}$ | $X_{s}<0.4$
V-dJ-NL | Vink _et al._ | A.1 | de Jager _et al._ | A.3 | Nugis & Lamers | A.6
V-dJ-H | Vink _et al._ | A.1 | de Jager _et al._ | A.3 | Hamann _et al._ | A.7
V-NJ-NL | Vink _et al._ | A.1 | Nieuwenhuijzen & de Jager | A.4 | Nugis & Lamers | A.6
V-NJ-H | Vink _et al._ | A.1 | Nieuwenhuijzen & de Jager | A.4 | Hamann _et al._ | A.7
V-vL-H | Vink _et al._ | A.1 | van Loon _et al._ | A.5 | Hamann _et al._ | A.7
V-vL-NL | Vink _et al._ | A.1 | van Loon _et al._ | A.5 | Nugis & Lamers | A.6
K-dJ-NL | Kudritzki _et al._ | A.2 | de Jager _et al._ | A.3 | Nugis & Lamers | A.6
K-dJ-H | Kudritzki _et al._ | A.2 | de Jager _et al._ | A.3 | Hamann _et al._ | A.7
K-NJ-NL | Kudritzki _et al._ | A.2 | Nieuwenhuijzen & de Jager | A.4 | Nugis & Lamers | A.6
K-NJ-H | Kudritzki _et al._ | A.2 | Nieuwenhuijzen & de Jager | A.4 | Hamann _et al._ | A.7
K-vL-NL | Kudritzki _et al._ | A.2 | van Loon _et al._ | A.5 | Nugis & Lamers | A.6
K-vL-H | Kudritzki _et al._ | A.2 | van Loon _et al._ | A.5 | Hamann _et al._ | A.7

The wind mass loss of massive stars is usually divided into three separate
phases, whose definition is somewhat arbitrary. When using the algorithm K for
the hot phase, we adopt the following thresholds based on the effective
temperature $T_{\mathrm{eff}}$ and surface (i.e., outermost computational
cell) hydrogen mass fraction $X_{s}$:

  * $\bullet$

Hot phase: $T_{\mathrm{eff}}\geq 15\,000\ \mathrm{K}$;

  * $\bullet$

Cool phase: $T_{\mathrm{eff}}<15\,000\ \mathrm{K}$;

  * $\bullet$

WR phase: $X_{s}<0.4\ \mathrm{regardless\ of}\ T_{\mathrm{eff}}$.

To follow the suggestions of Glebbeek et al. (2009), and to have a smoother
transition between the hot and cool phase wind algorithm, we use a slightly
different definition of the cool and hot evolutionary stages when using the V
mass loss algorithm for the hot phase:

  * $\bullet$

Hot phase: $T_{\mathrm{eff}}\geq 11\,000\ \mathrm{K}$;

  * $\bullet$

Cool phase: $T_{\mathrm{eff}}\leq 10\,000\ \mathrm{K}$;

and we use a linear interpolation between the hot phase wind and cool phase
wind in between. We choose this different threshold when using V to match how
the cool and hot phases are defined for the “Dutch” wind scheme in the MESA
code. Note that the interval from $T_{\mathrm{eff}}\simeq 10,\mathrm{kK}$ to
$\sim 15\,\mathrm{kK}$ is covered during a fraction of the Hertzsprung gap, in
a very short time.

The threshold dividing the cool and hot phases is qualitatively justified as
follows: the radiation pressure is determined by the product of opacity and
flux, which peaks between $10000$–$15000\ \mathrm{K}$ because of iron
recombination. Therefore, an effective temperature (i.e., in other words, a
radius, for each given luminosity) in this range is a physically meaningful
threshold to switch between wind mass loss algorithms.

The third phase is the WR phase and our criterion $X_{s}<0.4$ just requires a
hydrogen-poor stellar surface. This has very little in common with the
observational definition of what a WR star is, which is based on spectral
features that are not tracked by stellar evolution codes. Specifically, WR
stars are identified by their surface hydrogen depletion (Schmutz & Drissen,
1999) _and_ the presence of broad emission lines (van der Hucht, 2001;
Marchenko et al., 2010), indicating the presence of a wind with a steep
density and velocity gradient. Moreover, typical WR stars have high
$T_{\mathrm{eff}}$, and our definition might, in principle, produce
unrealistically cold (and red) WR stars. However, in absence of strong mixing
processes (e.g., due to rotation), the required surface hydrogen depletion can
only be reached by removing mass from the surface and revealing deep and hot
stellar layers. We do not find in our calculations cool but hydrogen depleted
models. WR stars are further subdivided into classes (WNH, WN, WC, WO, etc.)
based on the relative flux of specific lines. Here, we do not attempt to
distinguish between different WR sub-classes, because our simulations do not
produce the stellar spectra that would be necessary to distinguish these sub-
classes (see however Meynet & Maeder, 2003; Groh et al., 2014, and references
therein).

Our definitions of the evolutionary phases are commonly used in the literature
(see, e.g., Limongi & Chieffi, 2006; Eggenberger et al., 2007; Woosley &
Heger, 2007). We list in Tab. 2 the algorithm combinations explored here: each
combination is labeled by combining the abbreviations for each wind algorithm
introduced above.

To study the effects of the possible overestimate of the mass loss rate
caused, for example, by wind inhomogeneities (i.e., “clumpiness”), we use
three different values of the wind efficiency $\eta=1.0,\,0.33,\,0.1$. These
values span the range of observational estimates of the volume filling factor
of clumps (see Smith 2014 for a review). For simplicity, we use the same
$\eta$ during the entire evolution of a given model.

###  2.3 The Grid of Stellar Models

We employ release version 7624 of the open-source stellar evolution code MESA
(Paxton et al., 2011, 2013, 2015) and compute a grid of nonrotating solar
metallicity models. We choose $Z_{\odot}=0.019$ to match precisely the value
adopted in Vink et al. (2001). We consider initial masses of $15$, $20$, $25$,
$30$, and $35\,$ $M_{\odot}$. Higher initial mass models are more strongly
affected by numerical (and possibly physical) instabilities (see Appendix B
and references therein), which makes them less reliable for our purpose.
Appendix B describes the details of our MESA simulations that are not directly
related to mass loss. Here, we only mention that we use a 45-isotope nuclear
reaction network (mesa_45.net) until oxygen depletion, and a mixing length
parameter $\alpha_{\mathrm{mlt}}=2.0$ with exponential overshooting.

We run our grid of models in three steps and at each checkpoint we make a
selection of models to run at the next step222The models saved at each
checkpoint are available at https://zenodo.org/record/292924#.WK_eENWi60i .
The input files and customized routines are available at
https://stellarcollapse.org/renzo2017.. This selection is necessary to reduce
the total computational cost of our grid of models. In the first step, we
evolve our models from the zero age main sequence (ZAMS) to when the
temperature in the central computational cell rises above $T_{c}\geq
10^{9}\,\mathrm{K}$ (“end of the mass loss phase”). At this point, the star
has little time left to live (a few years) and lose mass through winds (for a
typical mass loss rate the star loses less than $10^{-4}\,M_{\odot}$). Also,
when $T_{c}\geq 10^{9}\,\mathrm{K}$ is reached, MESA starts to artificially
damp the mass loss rate for stability reasons. Mass loss is completely shut
off when $T_{c}\geq 2\times 10^{9}\,\mathrm{K}$.

Note, however, that the observation of SN impostors, type IIn SNe, and intra-
night flash spectroscopy of normal type II SNe suggests that non-wind mass
loss phenomena (neglected here) may occur in the very late phase of the
stellar life (see, e.g., Quataert & Shiode, 2012; Smith & Arnett, 2014;
Quataert et al., 2016; Khazov et al., 2016).

In the second step, we restart a subset of MESA models and run to oxygen
depletion (see Sec. 3.5). This is defined as the first time when the mass
fraction of ${}^{16}\mathrm{O}$ in the central cell drops below 0.04 and the
mass fraction of ${}^{28}\mathrm{Si}$ is higher than 0.01 (indicating that
some oxygen burning has already occurred, following Sukhbold & Woosley, 2014).
Note that these thresholds are an artificial choice, since the evolution of
the star is continuous throughout its lifetime.

Finally, we restart a select subset of stars at oxygen depletion and run to
the onset of core collapse, defined by

$\mathrm{max}\\{|v|\\}\geq 10^{3}\ \mathrm{km\ s^{-1}}\ \ ,$ (2)

where $v$ is the radial infall velocity (see, e.g., Heger et al., 2005;
Sukhbold & Woosley, 2014). For this last phase of the evolution, we switch
from the 45-isotopes nuclear reaction network to a customized 203-isotopes
nuclear reaction network (see also Appendix B), to capture the details of the
weak interaction physics (electron captures and $\beta$-decays) occurring
during silicon burning and before collapse. This physics determines the final
number of electrons per nucleon $Y_{e}$, thus the effective Chandrasekhar mass
of the stellar core, and ultimately the core structure at the onset of core
collapse. Experiments with a single-zone model for silicon burning show that
at least $\sim 200$ isotopes are required to obtain a converged final value of
$Y_{e}$ in the core (i.e., one that is independent of the size of the nuclear
reaction network; see also Farmer et al., 2016). Since MESA solves the fully
coupled set of equations for the chemical composition and structure of the
star, the increase in the number of isotopes forces us to reduce the number of
spatial mesh points in each model because of memory limitations (see Appendix
B.1.3).

##  3 Results

Figure 1: Uncertainty in the mapping between $M_{\mathrm{ZAMS}}$ and the
relative final mass $M/M_{\mathrm{ZAMS}}$ due to wind mass loss. Each colored
bar corresponds to a specific wind algorithm combination defined in Tab. 2.
The pluses, crosses, and circles correspond to $\eta=1.0,\,0.33,\,0.1$,
respectively. We employ the vertical bars to emphasize the spread. The
uncertainty in wind mass loss limits the predictive power of stellar evolution
studies for the final mass of stars of given initial mass. Only models with
$M_{\mathrm{ZAMS}}=35\,M_{\odot}$ and wind efficiency $\eta=1$ reach the WR
stage (cf. Tab. 4 and Sec. 3.4). They are shown in the rightmost panel and we
list the WR mass loss algorithm only for them. Note that the maximum relative
mass for the four algorithm combinations using the H WR mass loss algorithm in
the rightmost panel are the results obtained using the corresponding hot and
cool mass loss combination with $\eta=0.1$ (not reaching the WR phase).

###  3.1 Overview at the end of the mass loss phase, $T_{c}\geq
10^{9}\,\mathrm{K}$

Table 3: Model summary at the end of the mass loss phase (when $T_{c}\geq
10^{9}\,\mathrm{K}$, roughly corresponding to neon core ignition). We provide
total mass $M$, helium core mass $M_{\mathrm{He}}$, and carbon-oxygen (CO)
core mass $M_{\mathrm{CO}}$. We omit models differing only by the WR wind
scheme if they do not reach the WR stage (i.e. $X_{s}>0.4$ at all times). We
define the edge of the CO core as first location going inward where
$X(^{4}\mathrm{He})<0.01$, without requiring a minimum mass fraction of carbon
or oxygen. Analogously, we define the outer edge of the helium core as the
first location going inward where $X(^{1}\mathrm{H})<0.01$.  End of mass loss
phase: $T_{c}\geq 10^{9}\,\mathrm{K}$
---
$M_{\mathrm{ZAMS}}=15\,M_{\odot}$ | $M_{\mathrm{ZAMS}}=20\,M_{\odot}$
$\eta$ | ID | $M\ [M_{\odot}]$ | $M_{\mathrm{He}}\ [M_{\odot}]$ | $M_{\mathrm{CO}}\ [M_{\odot}]$ | $\eta$ | ID | $M\ [M_{\odot}]$ | $M_{\mathrm{He}}\ [M_{\odot}]$ | $M_{\mathrm{CO}}\ [M_{\odot}]$
0.1 | V-dJ | 14.66 |  04.99 |  03.20 | 0.1 | V-dJ | 19.23 |  07.04 | 4.91
V-NJ | 14.64 |  04.99 |  03.20 | V-NJ | 19.23 |  07.04 | 4.91
V-vL | 13.60 |  04.98 |  03.20 | V-vL | 18.10 |  07.03 | 4.91
K-dJ | 14.67 |  05.00 |  03.22 | K-dJ | 19.38 |  07.01 | 4.90
K-NJ | 14.66 |  05.00 |  03.22 | K-NJ | 19.37 |  07.01 | 4.90
K-vL | 13.61 |  05.00 |  03.21 | K-vL | 18.67 |  07.01 | 4.90
0.33 | V-dJ | 13.94 |  04.93 |  03.16 | 0.33 | V-dJ | 17.47 |  07.01 | 4.88
V-NJ | 13.90 |  04.93 |  03.16 | V-NJ | 17.48 |  07.01 | 4.88
V-vL | 10.39 |  04.93 |  03.15 | V-vL | 13.47 |  06.99 | 4.87
K-dJ | 13.92 |  04.98 |  03.19 | K-dJ | 17.62 |  07.00 | 4.87
K-NJ | 13.88 |  04.98 |  03.19 | K-NJ | 17.62 |  07.00 | 4.88
K-vL | 10.11 |  04.97 |  03.19 | K-vL | 13.90 |  06.98 | 4.87
1.0 | V-dJ | 12.86 |  04.65 |  02.95 | 1.0 | V-dJ | 11.81 |  07.06 | 4.92
V-NJ | 12.74 |  04.65 |  02.95 | V-NJ | 12.04 |  07.06 | 4.92
V-vL |  05.25 |  04.64 |  02.94 | V-vL |  08.80 |  07.02 | 4.89
K-dJ | 11.97 |  04.92 |  03.15 | K-dJ | 12.62 |  06.95 | 4.84
K-NJ | 11.87 |  04.92 |  03.14 | K-NJ | 12.77 |  06.95 | 4.84
K-vL |  05.70 |  04.90 |  03.13 | K-vL |  08.81 |  06.90 | 4.80
$M_{\mathrm{ZAMS}}=25\,M_{\odot}$ | $M_{\mathrm{ZAMS}}=30\,M_{\odot}$
$\eta$ | ID | $M\ [M_{\odot}]$ | $M_{\mathrm{He}}\ [M_{\odot}]$ | $M_{\mathrm{CO}}\ [M_{\odot}]$ | $\eta$ | ID | $M\ [M_{\odot}]$ | $M_{\mathrm{He}}\ [M_{\odot}]$ | $M_{\mathrm{CO}}\ [M_{\odot}]$
0.1 | V-dJ | 23.85 |  09.14 |  06.58 | 0.1 | V-dJ | 28.11 | 10.97 | 7.97
V-NJ | 23.86 |  09.14 |  06.59 | V-NJ | 28.17 | 10.97 | 7.94
V-vL | 23.56 |  09.14 |  06.59 | V-vL | 28.70 | 10.96 | 7.95
K-dJ | 23.74 |  09.23 |  06.66 | K-dJ | 28.13 | 10.91 | 7.90
K-NJ | 23.76 |  09.24 |  06.67 | K-NJ | 28.19 | 10.91 | 7.86
K-vL | 23.24 |  09.23 |  06.67 | K-vL | 28.75 | 10.90 | 7.91
0.33 | V-dJ | 21.98 |  08.87 |  06.37 | 0.33 | V-dJ | 24.16 | 11.11 | 8.13
V-NJ | 22.02 |  08.87 |  06.38 | V-NJ | 24.34 | 11.11 | 8.13
V-vL | 22.13 |  08.87 |  06.36 | V-vL | 26.02 | 11.10 | 8.12
K-dJ | 21.27 |  09.11 |  06.57 | K-dJ | 25.03 | 10.89 | 7.88
K-NJ | 21.34 |  09.11 |  06.57 | K-NJ | 25.18 | 10.89 | 7.87
K-vL | 20.26 |  09.11 |  06.55 | K-vL | 27.12 | 10.87 | 7.80
1.0 | V-dJ | 13.29 |  09.05 |  06.50 | 1.0 | V-dJ | 15.36 | 11.28 | 8.23
V-NJ | 13.67 |  09.05 |  06.48 | V-NJ | 15.78 | 11.29 | 8.26
V-vL | 12.97 |  09.04 |  06.51 | V-vL | 18.29 | 11.31 | 8.27
K-dJ | 15.57 |  08.89 |  06.38 | K-dJ | 18.51 | 10.92 | 7.96
K-NJ | 15.84 |  08.89 |  06.37 | K-NJ | 18.80 | 10.92 | 7.94
K-vL | 14.69 |  08.90 |  06.38 | K-vL | 22.53 | 10.89 | 7.94
$M_{\mathrm{ZAMS}}=35\,M_{\odot}$ | a
$\eta$ | ID | $M\ [M_{\odot}]$ | $M_{\mathrm{He}}\ [M_{\odot}]$ | $M_{\mathrm{CO}}\ [M_{\odot}]$
0.1 | V-dJ-NL | 32.61 | 10.25 |  07.19
V-NJ-NL | 32.69 | 10.24 |  07.19
V-vL-NL | 33.99 | 10.03 |  07.03
K-dJ-NL | 32.43 | 11.87 |  08.54
K-NJ-NL | 32.53 | 11.87 |  08.56
K-vL-NL | 33.97 | 11.83 |  08.51
0.33 | V-dJ-NL | 27.36 | 12.72 |  09.32
V-NJ-NL | 27.62 | 12.72 |  09.32
V-vL-NL | 31.40 | 12.67 |  09.23
K-dJ-NL | 27.50 | 12.68 |  09.21
K-NJ-NL | 27.78 | 12.68 |  09.28
K-vL-NL | 31.76 | 12.63 |  09.21
1.0 | V-dJ-NL | 20.03 | 13.30 |  09.81
V-NJ-NL | 20.20 | 13.30 |  09.94
V-vL-NL | 24.85 | 13.30 |  09.97
K-dJ-NL | 20.32 | 13.82 | 10.33
K-NJ-NL | 20.45 | 13.82 | 10.36
K-vL-NL | 25.13 | 13.82 | 10.41
V-dJ-H | 19.73 | 13.30 |  09.93
V-NJ-H | 19.93 | 13.30 |  09.99
K-dJ-H | 19.48 | 13.81 | 10.41 |  |  |  |  |
K-NJ-H | 19.69 | 13.82 | 10.41 |  |  |  |  |
Table 4: Final color and surface properties of the computed models. We follow
the definitions of Georgy (2012) to classify models as RSGs, YSGs, or BSGs
(see also text). WR stars have $X_{s}\leq 0.4$, regardless of their surface
temperature. The first and second columns indicate the wind efficiency and the
mass loss algorithm combination, respectively. The WR stars are computed
twice, once with the NL mass loss algorithm (see Sec. A.6) and once with the H
algorithm (see Sec. A.7).  End of mass loss phase: $T_{c}\geq
10^{9}\,\mathrm{K}$
---
$M_{\mathrm{ZAMS}}=15\,M_{\odot}$ | $M_{\mathrm{ZAMS}}=20\,M_{\odot}$
$\eta$ | ID | $R\ [R_{\odot}]$ | $\log_{10}(L/L_{\odot})$ | $\log_{\mathrm{10}}(T_{\mathrm{eff}}/\mathrm{[K]})$ | color | $\eta$ | ID | $R\ [R_{\odot}]$ | $\log_{10}(L/L_{\odot})$ | $\log_{\mathrm{10}}(T_{\mathrm{eff}}/\mathrm{[K]})$ | color
0.1 | V-dJ | 911 | 5.06 | 3.55 | RSG | 0.1 | V-dJ | 992 | 5.27 | 3.58 | RSG
V-NJ | 911 | 5.06 | 3.55 | RSG | V-NJ | 992 | 5.27 | 3.58 | RSG
V-vL | 927 | 5.06 | 3.54 | RSG | V-vL | 996 | 5.27 | 3.58 | RSG
K-dJ | 914 | 5.06 | 3.55 | RSG | K-dJ | 987 | 5.27 | 3.58 | RSG
K-NJ | 914 | 5.06 | 3.55 | RSG | K-NJ | 989 | 5.27 | 3.58 | RSG
K-vL | 930 | 5.06 | 3.54 | RSG | K-vL | 991 | 5.27 | 3.58 | RSG
0.33 | V-dJ | 916 | 5.05 | 3.54 | RSG | 0.33 | V-dJ | 998 | 5.27 | 3.58 | RSG
V-NJ | 916 | 5.05 | 3.54 | RSG | V-NJ | 999 | 5.27 | 3.58 | RSG
V-vL | 968 | 5.05 | 3.53 | RSG | V-vL | 985 | 5.27 | 3.58 | RSG
K-dJ | 920 | 5.06 | 3.54 | RSG | K-dJ | 998 | 5.27 | 3.58 | RSG
K-NJ | 921 | 5.06 | 3.54 | RSG | K-NJ | 998 | 5.27 | 3.58 | RSG
K-vL | 959 | 5.06 | 3.54 | RSG | K-vL | 991 | 5.27 | 3.58 | RSG
1.0 | V-dJ | 895 | 5.02 | 3.54 | RSG | 1.0 | V-dJ | 950 | 5.28 | 3.59 | RSG
V-NJ | 896 | 5.02 | 3.54 | RSG | V-NJ | 956 | 5.28 | 3.59 | RSG
V-vL | 630 | 5.02 | 3.62 | YSG | V-vL | 672 | 5.28 | 3.67 | YSG
K-dJ | 947 | 5.05 | 3.54 | RSG | K-dJ | 972 | 5.27 | 3.59 | RSG
K-NJ | 949 | 5.05 | 3.54 | RSG | K-NJ | 976 | 5.27 | 3.58 | RSG
K-vL | 644 | 5.05 | 3.62 | YSG | K-vL | 712 | 5.27 | 3.65 | YSG
$M_{\mathrm{ZAMS}}=25\,M_{\odot}$ | $M_{\mathrm{ZAMS}}=30\,M_{\odot}$
$\eta$ | ID | $R\ [R_{\odot}]$ | $\log_{10}(L/L_{\odot})$ | $\log_{\mathrm{10}}(T_{\mathrm{eff}}/\mathrm{[K]})$ | color | $\eta$ | ID | $R\ [R_{\odot}]$ | $\log_{10}(L/L_{\odot})$ | $\log_{\mathrm{10}}(T_{\mathrm{eff}}/\mathrm{[K]})$ | color
0.1 | V-dJ | 899 | 5.41 | 3.64 | YSG | 0.1 | V-dJ | 698 | 5.50 | 3.72 | YSG
V-NJ | 900 | 5.42 | 3.64 | YSG | V-NJ | 697 | 5.51 | 3.72 | YSG
V-vL | 902 | 5.42 | 3.64 | YSG | V-vL | 696 | 5.51 | 3.72 | YSG
K-dJ | 891 | 5.43 | 3.64 | YSG | K-dJ | 705 | 5.53 | 3.72 | YSG
K-NJ | 891 | 5.44 | 3.65 | YSG | K-NJ | 704 | 5.53 | 3.72 | YSG
K-vL | 893 | 5.43 | 3.64 | YSG | K-vL | 708 | 5.52 | 3.72 | YSG
0.33 | V-dJ | 926 | 5.42 | 3.63 | YSG | 0.33 | V-dJ | 670 | 5.54 | 3.73 | YSG
V-NJ | 924 | 5.41 | 3.63 | YSG | V-NJ | 670 | 5.54 | 3.73 | YSG
V-vL | 932 | 5.42 | 3.63 | YSG | V-vL | 675 | 5.55 | 3.73 | YSG
K-dJ | 896 | 5.42 | 3.64 | YSG | K-dJ | 719 | 5.53 | 3.72 | YSG
K-NJ | 897 | 5.42 | 3.64 | YSG | K-NJ | 717 | 5.52 | 3.71 | YSG
K-vL | 899 | 5.43 | 3.64 | YSG | K-vL | 711 | 5.53 | 3.72 | YSG
1.0 | V-dJ | 797 | 5.43 | 3.67 | YSG | 1.0 | V-dJ | 614 | 5.56 | 3.76 | YSG
V-NJ | 811 | 5.43 | 3.67 | YSG | V-NJ | 628 | 5.56 | 3.75 | YSG
V-vL | 782 | 5.43 | 3.67 | YSG | V-vL | 610 | 5.55 | 3.76 | YSG
K-dJ | 874 | 5.42 | 3.65 | YSG | K-dJ | 755 | 5.54 | 3.71 | YSG
K-NJ | 878 | 5.42 | 3.65 | YSG | K-NJ | 758 | 5.55 | 3.71 | YSG
K-vL | 857 | 5.42 | 3.65 | YSG | K-vL | 723 | 5.53 | 3.72 | YSG
$M_{\mathrm{ZAMS}}=35\,M_{\odot}$ | a
$\eta$ | ID | $R\ [R_{\odot}]$ | $\log_{10}(L/L_{\odot})$ | $\log_{\mathrm{10}}(T_{\mathrm{eff}}/\mathrm{[K]})$ | color
0.1 | V-dJ | 860 | 5.53 | 3.68 | YSG
V-NJ | 861 | 5.53 | 3.68 | YSG
V-vL | 920 | 5.52 | 3.66 | YSG
K-dJ | 726 | 5.59 | 3.73 | YSG
K-NJ | 730 | 5.59 | 3.73 | YSG
K-vL | 755 | 5.58 | 3.72 | YSG
0.33 | V-dJ | 529 | 5.62 | 3.80 | BSG
V-NJ | 531 | 5.62 | 3.80 | BSG
V-vL | 569 | 5.61 | 3.79 | YSG
K-dJ | 539 | 5.62 | 3.80 | BSG
K-NJ | 541 | 5.62 | 3.80 | BSG
K-vL | 595 | 5.62 | 3.78 | YSG
1.0 | V-dJ-NL | 258 | 5.53 | 3.94 | WR
V-NJ-NL | 255 | 5.66 | 3.97 | WR
V-vL | 398 | 5.66 | 3.88 | BSG
K-dJ-NL | 167 | 5.75 | 4.09 | WR
K-NJ-NL | 168 | 5.73 | 4.08 | WR
K-vL | 313 | 5.67 | 3.93 | BSG
V-dJ-H | 253 | 5.60 | 3.96 | WR
V-NJ-H | 260 | 5.65 | 3.97 | WR
K-dJ-H | 157 | 5.72 | 4.09 | WR |  |  |  |  |  |
K-NJ-H | 158 | 5.72 | 4.09 | WR |  |  |  |  |  |
Table 5: Maximum spread in total mass, core masses, and radius at the end of the mass loss phase, for models differing in mass loss algorithm combination and efficiency $\eta$. We also list the maximum and minimum total mass for each $M_{\mathrm{ZAMS}}$. $M_{\mathrm{ZAMS}}$ | $\max\Delta M$ | $\max M$ | $\min M$ | $\max\Delta M_{\mathrm{He}}$ | $\max\Delta M_{\mathrm{CO}}$ | $\max\Delta R$
---|---|---|---|---|---|---
$[M_{\odot}]$ | $[M_{\odot}]$ | $[M_{\odot}]$ | $[M_{\odot}]$ | $[M_{\odot}]$ | $[M_{\odot}]$ | $[R_{\odot}]$
15 |  09.42 | 14.66 |  05.25 | 0.36 | 0.28 | 338
20 | 10.58 | 19.38 |  08.80 | 0.16 | 0.12 | 327
25 | 10.89 | 23.86 | 12.97 | 0.37 | 0.31 | 150
30 | 13.39 | 28.75 | 15.36 | 0.44 | 0.47 | 148
35 | 14.51 | 33.99 | 19.48 | 3.79 | 3.38 | 763

The lifetime remaining for the star after $T_{c}\geq 10^{9}\,\mathrm{K}$ is
short ($\sim 15$ years for $15\,M_{\odot}$, $\sim 4$ years for
$30\,M_{\odot}$), and the photosphere of the star remains frozen at the same
effective temperature $T_{\mathrm{eff}}$ and luminosity $L$ as long as no
impulsive mass loss events or other instabilities take place. The subsequent
evolution will lead to changes of the internal structure only. Therefore, we
discuss the final mass and pre-SN appearance of our models already at this
stage.

Figure 1 shows the overall impact of the uncertainty in stellar wind mass loss
rates on the final mass. We plot the final mass relative to the initial mass
for each considered $M_{\mathrm{ZAMS}}$ and the large vertical spread
illustrates the uncertainty. As an example, we can consider a star of
$M_{\mathrm{ZAMS}}=20\,M_{\odot}$. It can, in principle, reach the onset of
collapse with $M=19.38\,M_{\odot}$ if evolved with the V-NJ mass loss
combination with reduced efficiency $\eta=0.1$. But it might also evolve to
$M=8.81\,M_{\odot}$ with the K-vL combination and efficiency $\eta=1.0$. If
these are truly limiting cases and anything in between is unconstrained, then
the uncertainty in the final mass is greater than $50\%$. For more
quantitative results, see Tab. 5 where we report the maximum spreads of the
ZAMS to pre-SN mass mapping for our entire grid of models. Tab. 3 and Tab. 4
list the main physical quantities for all the stars in our grid at the end of
the mass loss phase.

It is important to note that the vertical spread in the ZAMS to pre-SN mass
mapping shown in Fig. 1 is dominated by the highly uncertain wind efficiency
$\eta$. At fixed $\eta$, variations due to different wind algorithm
combinations are minor. _This makes $\eta$ the most important free parameter
for wind mass loss in stellar evolution calculations._

In Fig. 1, the spread in $M_{\mathrm{ZAMS}}$ to pre-SN mass decreases for
higher $M_{\mathrm{ZAMS}}$. This is, however, only because we show the
relative final mass, i.e. at higher $M_{\mathrm{ZAMS}}$ we divide the final
mass by a larger number. In absolute numbers, the uncertainty in the final to
initial mass relation increases for more massive stars. We summarize this in
Tab. 3. As an example, the maximum spread between the final masses of
$35\,M_{\odot}$ models is $\max{\Delta M}\equiv
33.99\,M_{\odot}-19.48\,M_{\odot}=14.51\,M_{\odot}$. In the $15\,M_{\odot}$
case, it is $\max{\Delta M}\equiv
14.66\,M_{\odot}-5.25\,M_{\odot}=9.42\,M_{\odot}$. As expected, all wind mass
loss algorithm combinations yield a higher mass loss rate for more massive
(and thus more luminous) stars. For stars with $M_{\mathrm{ZAMS}}\lesssim
20\,M_{\odot}$, the models using the vL algorithm (dust driven mass loss in
the cool evolutionary phase) produce much higher mass loss than all other
algorithms (see also Sec. 3.3): for example, the $15\,M_{\odot}$ model using
the V-vL combination with full efficiency $\eta=1.0$ results in a pre-SN mass
of only $5.25\,M_{\odot}$, while using the combination V-dJ or V-NJ, we obtain
a final mass of $\sim 12.7\,M_{\odot}$. We discuss this effect in more detail
in Sec. 3.3.

From Fig. 1, we also note that for $M_{\mathrm{ZAMS}}\gtrsim 20\,M_{\odot}$ we
obtain a range $1.0\lesssim M/M_{\mathrm{ZAMS}}\lesssim 0.5$. For any given
$\eta$, the range is smaller. In any case, the size of this interval does not
decrease going to higher masses, indicating that the different mass loss
prescriptions and efficiencies do not appear to be converging with
$M_{\mathrm{ZAMS}}$ in the mass range we consider here.

Figure 2: Total stellar mass as a function of time for the $15M_{\odot}$
models. Each color corresponds to a given combination of wind algorithms (see
Tab. 2). None of these $15\,M_{\odot}$ models reach the WR phase. The solid,
dashed, and dot-dashed curves correspond to efficiency
$\eta=0.1,\,0.33,\,1.0$, respectively, regardless of color. The red vertical
line indicates roughly the terminal age main sequence (TAMS, i.e.
$X_{c}<0.01$). The enhanced mass loss of models using the Vink et al. (2000,
2001) (V) algorithm close to TAMS in the left panel is due to the bistability
jump, see Sec. 3.2. Note that the left panel has a different vertical scale to
magnify the differences in the tracks during the main sequence evolution.

In Fig. 2, we show, as a representative example, the time evolution of the
total mass for our $15\,M_{\odot}$ models. The overall qualitative behavior is
the same for all the other considered $M_{\mathrm{ZAMS}}$. Fig. 2 shows that
the amount of mass lost during the main sequence (i.e. prior to the vertical
dot-dashed line) is relatively small, only a few percent of the total mass for
$M_{\mathrm{ZAMS}}=15\,M_{\odot}$, even when using $\eta=1.0$. Most of the
mass is lost after hydrogen core burning. Both the mass loss rate and the
spread between the predictions of different algorithms increases dramatically
after the end of the main sequence. Hence, uncertainties in the post-main-
sequence cool phase mass loss rates have by far the greatest effect on the
final mass of the star.

Using our nonrotating single-star models at the end of the mass loss phase and
ignoring any potential subsequent non-wind mass loss events, we can attempt to
classify their pre-SN color. We follow Georgy (2012) by assuming
$\log_{10}(T_{\mathrm{eff}}/\mathrm{[K]})\lesssim 3.6$ for RSG,
$3.6\lesssim\log_{10}(T_{\mathrm{eff}}/\mathrm{[K]})\lesssim 3.8$ for YSG,
$\log_{10}(T_{\mathrm{eff}}/\mathrm{[K]})\gtrsim 3.8$ for BSG, and we only
require the surface abundance of hydrogen to be $X_{s}<0.4$ without
temperature thresholds for WR stars. We list the outcome of this
classification in Tab. 4. Almost all of our $15\,M_{\odot}$ and
$20\,M_{\odot}$ models end their evolution as RSGs. The exceptions are those
using the vL cool mass loss algorithm and $\eta=1.0$, which end their lives as
YSGs. All of our $25\,M_{\odot}$ and $30\,M_{\odot}$ models end as YSGs. The
$35\,M_{\odot}$ models computed with $\eta=0.1$ also end as YSGs. With
$\eta=0.33$ they instead become BSGs, because of the higher mass loss rate,
except when using the vL cool mass loss rate (cf. Sec. 3.3), which produces
YSGs. With $\eta=1.0$ we find WR pre-SN models, unless vL is used during the
cool phase, in which case our $35\,M_{\odot}$ models would explode as BSGs.
However, these results are highly dependent on the somewhat arbitrary
temperature thresholds assumed to divide the categories. For example, assuming
$\log_{10}(T_{\mathrm{eff}}/\mathrm{[K]})\leq 3.68$ as the threshold dividing
RSG and YSG, all models with $M_{\mathrm{ZAMS}}\leq 25\,M_{\odot}$ would be
RSGs.

Figure 3: Mass loss rate during the hot evolutionary phase (including the main
sequence, see Sec. 2.2) as a function of time for all computed
$M_{\mathrm{ZAMS}}$ and wind efficiency $\eta=1$. The solid and dashed curves
are computed using the V (Vink et al., 2000, 2001) and K (Kudritzki et al.,
1989) mass loss algorithm, respectively. The rapid rise of the solid curves is
due to the inclusion of the bistability jump (see Sec. 3.2). The qualitative
behavior of the curves shown does not change substantially when varying
$\eta$. The curves end when $T_{\mathrm{eff}}=15\,000\,\mathrm{K}$.

###  3.2 The “hot phase” mass loss

Table 6:  Stellar properties at the end of hot phase, when $T_{\mathrm{eff}}$ decreases below $15\,000\,\mathrm{K}$ for the first time. For the descriptions of the V and K algorithms see Sec. A.1 and Sec. A.2, respectively.  Hot wind algorithm | $M_{\mathrm{ZAMS}}\ [M_{\odot}]$ | $\eta$ | $R\ [R_{\odot}]$ | $L\ [10^{4}L_{\odot}]$ | $M\ [M_{\odot}]$ | $M_{\mathrm{He}}\ [M_{\odot}]$ | age [Myr]
---|---|---|---|---|---|---|---
V | 15 | 1.00 |  034.85 |  05.44 | 14.54 |  03.69 | 12.7917
K | 15 | 1.00 |  035.44 |  05.64 | 14.73 |  03.71 | 12.8239
V | 15 | 0.33 |  035.91 |  05.76 | 14.85 |  03.76 | 12.7503
K | 15 | 0.33 |  035.97 |  05.82 | 14.91 |  03.77 | 12.7612
V | 15 | 0.10 |  036.27 |  05.85 | 14.95 |  03.78 | 12.7363
K | 15 | 0.10 |  035.96 |  05.83 | 14.97 |  03.78 | 12.7359
V | 20 | 1.00 |  051.97 | 12.20 | 19.06 |  05.83 |  08.8937
K | 20 | 1.00 |  052.75 | 12.32 | 19.42 |  05.87 |  08.8947
V | 20 | 0.33 |  053.10 | 12.67 | 19.69 |  05.95 |  08.8373
K | 20 | 0.33 |  053.35 | 12.76 | 19.81 |  05.95 |  08.8374
V | 20 | 0.10 |  053.46 | 12.80 | 19.91 |  05.99 |  08.8176
K | 20 | 0.10 |  054.39 | 13.10 | 19.94 |  05.98 |  08.8183
V | 25 | 1.00 |  067.58 | 20.63 | 23.24 |  08.03 |  07.0367
K | 25 | 1.00 |  068.70 | 21.41 | 24.00 |  07.91 |  07.0094
V | 25 | 0.33 |  070.59 | 22.55 | 24.44 |  07.97 |  06.9653
K | 25 | 0.33 |  070.05 | 22.04 | 24.66 |  08.07 |  06.9546
V | 25 | 0.10 |  070.88 | 22.39 | 24.84 |  08.10 |  06.9390
K | 25 | 0.10 |  070.50 | 22.19 | 24.90 |  08.13 |  06.9355
V | 30 | 1.00 |  080.82 | 29.52 | 27.00 | 10.18 |  05.9700
K | 30 | 1.00 |  084.32 | 31.82 | 28.46 | 10.00 |  05.9215
V | 30 | 0.33 |  084.64 | 32.11 | 29.10 | 10.16 |  05.8853
K | 30 | 0.33 |  085.14 | 32.83 | 29.48 | 10.11 |  05.8689
V | 30 | 0.10 |  085.58 | 32.85 | 29.77 |  09.96 |  05.8558
K | 30 | 0.10 |  085.34 | 32.93 | 29.84 |  09.89 |  05.8503
V | 35 | 1.00 |  095.27 | 40.28 | 30.32 | 12.29 |  05.2807
K | 35 | 1.00 |  097.93 | 41.83 | 32.83 | 12.51 |  05.2164
V | 35 | 0.33 |  098.15 | 43.52 | 33.62 | 12.19 |  05.1879
K | 35 | 0.33 |  099.84 | 44.24 | 34.27 | 12.09 |  05.1661
V | 35 | 0.10 | 100.81 | 45.59 | 34.63 | 11.58 |  05.1559
K | 35 | 0.10 | 100.35 | 45.03 | 34.78 | 11.78 |  05.1484

“Hot phase” evolution, i.e. $T_{\mathrm{eff}}\gtrsim 15\,000\,\mathrm{K}$ if
using the Kudritzki et al. (1989) (K) mass loss algorithm, or
$T_{\mathrm{eff}}\gtrsim 11\,000\,\mathrm{K}$ if using Vink et al. (2000,
2001) (V), roughly covers the main sequence evolution, the subsequent overall
contraction caused by hydrogen depletion (known as the main sequence “hook”),
and the initial part of the Hertzsprung gap.

In Tab. 6, we summarize key quantitative results for our models at the end of
the hot phase. While this is the longest phase of the evolution, covering
$\gtrsim 90\%$ of the stellar lifetime, the amount of mass lost during this
phase is relatively small regardless of the algorithm used. For example,
initially $35\,M_{\odot}$ stars computed with $\eta=1.0$ lose $\sim 15\%$ of
their mass, while initially $15\,M_{\odot}$ stars only lose a few percent

We can also infer from Tab. 6 that the amount of mass lost during the hot
phase is always higher with V than with K, and the difference increases with
increasing $\eta$ and $M_{\mathrm{ZAMS}}$: it is only $\sim 0.02\,M_{\odot}$
for $\eta=0.1$ and $M_{\mathrm{ZAMS}}=15\,M_{\odot}$, and grows to $\sim
0.2\,M_{\odot}$ for $\eta=1.0$ and the same initial mass. For the
$35\,M_{\odot}$ models and $\eta=1.0$, the difference between the total mass
shed using V or K reaches $\sim 0.5\,M_{\odot}$ with $\eta=1.0$.

In Fig. 3, we plot the mass loss rates $\dot{M}$ given by V and K as a
function of time in the hot phase for models with $\eta=1$. The reason for the
higher total mass lost with V is that this algorithm includes a detailed
treatment of the bistability jump, which is an increase in the cross section
for photon interactions caused by the recombination of ions driving the mass
loss. This enhancement of the cross section happens when the effective
temperature drops below $T_{\mathrm{eff}}^{\mathrm{jump}}\simeq
25\,000\,\mathrm{K}$, Vink et al. (2000). This is what causes the sudden
tremendous increase in the mass loss rate in V models seen in Fig. 3. The
subsequent drop in models with masses higher than $20\,M_{\odot}$ happens
because these cross the bistability jump region twice during the contraction
following core hydrogen depletion. Overall, the average mass loss rate
$\langle|\dot{M}_{\mathrm{hot}}|\rangle$ of V models is driven up and
surpasses that of K models on the main sequence, resulting in higher mass loss
in V models. Note that the use of two different thresholds to separate the hot
and cool phase of evolution for V and K does not significantly influence the
duration of the hot phase: the gap between the two thresholds is covered in a
fraction of the Herzsprung gap duration.

Different values of $\eta$ produce small ($\lesssim 2\%$) age differences:
models remaining more massive (i.e. computed with lower $\eta$) evolve
slightly faster. However, these differences are too small to potentially be
used as observational tests for the wind efficiency.

Figure 3 also shows that at any point in time, the mass loss rate of more
massive stars is higher, because they produce a higher photon flux to drive
the wind. A factor of $\sim 2.3\left(\simeq 35/15\right)$ in initial mass
translates to a difference of almost two orders of magnitude in the mass loss
rate. The difference between the V and K rates is a non-monotonic function of
the mass, because of the different functional dependencies of the two
algorithms: for $20\,M_{\odot}$ models they are roughly equal before the
surface cools enough for the bistability jump to occur; for $15\,M_{\odot}$
models, the K algorithm gives an initially higher mass loss rate, and, for
$30\,M_{\odot}$ models, K mass loss is instead lower.

Although only a small amount of mass is lost, the hot phase mass loss can
significantly influence the core evolution. This is because no shell sources
decouple the surface from the convective core during most of this phase. The
use of different algorithms during the hot phase can thus create small (seed)
differences in the core structure, which may then be amplified by the
subsequent evolution of the star and contraction of the core. These
differences are small at the end of the hot phase, and we will discuss them at
later stages in the evolution in Sec. 3.6 and 3.7.

In Tab. 6, we also list the helium core masses $M_{\mathrm{He}}$ of our models
at the end of the hot phase. The effects of hot phase wind mass loss on
$M_{\mathrm{He}}$ are less straightforward to interpret than its influence on
the total mass. First, the value of $M_{\mathrm{He}}$ depends on the
definition of the helium core. We define the outer edge of the helium core as
the first location going inward where $X(^{1}\mathrm{H})<0.01$. Second, where
this interface is at the end of the hot phase is very sensitive to mixing:
depending on the mass and metallicity of the star (and the convective
stability criterion adopted), deep convective shells can develop at the
beginning of the Hertzsprung gap, and they shape the chemical composition
profile and determine $M_{\mathrm{He}}$. Tab. 6 shows that the maximum spread
in $M_{\mathrm{He}}$ increases with $M_{\mathrm{ZAMS}}$, starting from
$\max(\Delta M_{\mathrm{He}})\simeq 0.1\,M_{\odot}$ for $15\,M_{\odot}$
models, up to $\max(\Delta M_{\mathrm{He}})\simeq 0.9\,M_{\odot}$ for
$35\,M_{\odot}$ models. The spread in these values is almost entirely due to
variations in $\eta$. The difference in $M_{\mathrm{He}}$ between models of
same mass and $\eta$ (thus differing only in the use of V or K) is also
increasing with increasing $M_{\mathrm{ZAMS}}$ but remains below $\sim
0.2\,M_{\odot}$. This trend directly reflects the larger uncertainties in the
modeling of winds from more massive stars.

###  3.3 The “cool phase” mass loss

Regardless of efficiency $\eta$ and mass loss algorithm, most of the mass loss
through stellar winds happens during the cool phase of the evolution. Fig. 2
demonstrates this clearly for the $15\,M_{\odot}$ models, which lose
$2\%-60\%$ of their total initial mass during this phase, depending on the
algorithm combination and wind efficiency. The increase in the mass loss rate
from the hot phase can be understood in terms of the effective gravity of the
star (although we stress that the algorithms compared here do not depend
explicitly on it): for any given luminosity of a massive star, if the stellar
surface is cool, necessarily its radius must be large, and thus it will be
easier for matter to leave the gravitational potential well of the star.
Moreover, at lower temperature the opacity tends to be higher because of
recombination of ions and possibly dust formation, thus enhancing the wind
driving.

Table 3 summarizes key quantitative results of our models at the end of the
mass loss phase, including the total mass and the core masses. One striking
result is that the dust-driven van Loon et al. (2005) (vL) mass loss algorithm
results in significantly different total masses and core masses than the de
Jager et al. (1988) (dJ) and the Nieuwenhuijzen & de Jager (1990) (NJ)
algorithms. This can also be seen in Fig. 1, where wind combinations using vL
for the cool phase produce different vertical spreads. These differences are
strongest for $\eta=1.0$. The most extreme example are $15\,M_{\odot}$ models
computed with $\eta=1.0$: regardless of the hot phase mass loss algorithm,
they end their evolution with masses of $\sim 5-6\,M_{\odot}$ with vL in the
cool phase, while they remain as massive as $\sim 11-12\,M_{\odot}$ with dJ or
NJ. We find the opposite trend at the upper end of our mass range: a
$30\,M_{\odot}$ star computed with $\eta=1.0$ and V during the hot phase
reaches the end of the mass loss phase with a total mass of $\sim 15M_{\odot}$
if using either dJ or NJ in the cool phase, while it ends its life with $\sim
18\,M_{\odot}$ with vL.

On the one hand, the similarities between the dJ and the NJ rates are expected
(see also Mauron & Josselin, 2011; Eldridge & Tout, 2004): both are semi-
empirical rates derived from the same sample of observed stars. They differ
only in the choice of the stellar variables used to parametrize $\dot{M}$. On
the other hand, the vL algorithm is also semi-empirical, but based on the
analysis of a different sample of stars assuming a dust-driven model of the
wind, i.e. wind mass loss is not driven by photons impinging on metallic ions,
but rather on dust particles. If a dust-driven (instead of line-driven) model
of the wind is assumed, the resulting mass loss rate is generally higher, and
much more $T_{\mathrm{eff}}$-dependent (see Tabs. 1 and 9).

The very steep dependence of the vL rate on $T_{\mathrm{eff}}$ causes the
different evolution of models above and below $\sim 25\,M_{\odot}$. During the
early stage of the cool phase, the vL rate is always much higher than the
others, and the stellar wind described by this algorithm reveals the deeper
and hotter layers of the star (see also Tab. 4). As $T_{\mathrm{eff}}$
increases, the vL mass loss rate decreases rapidly ($\propto T^{-6.3}$, see
Tabs. 1 and 9), which is attributed to the temperature sensitivity of the
microscopic dust formation processes (Wachter et al., 2002) that we of course
do not track explicitly in our calculations. Also, for any mass loss process
at a given luminosity, higher $T_{\mathrm{eff}}$ correspond to smaller radii,
and thus higher effective gravity at the stellar surface. For
$M_{\mathrm{ZAMS}}\gtrsim 25\,M_{\odot}$, the steep
$T_{\mathrm{eff}}$-dependence leads to a self regulation of the vL rate. For
lower initial masses, the vL rate is also initially higher than the dJ or NJ
rate, but not high enough to reach the self-regulating regime: the vL rate
remains higher than the dJ and NJ rates for the whole evolution and produces
pre-SN structures of a much lower final mass than when dJ or NJ are used. This
is summarized in Tab. 3.

The comparison of $M_{\mathrm{He}}$ listed in Tab. 3 for models with the same
$\eta$ and hot phase mass loss (either V or K) reveals that the effect of the
cool phase mass loss on the helium core mass is very small and almost
negligible. We find the only appreciable differences when using $\eta=1.0$ and
the vL cool mass loss rate, and they are only of order $\sim 0.01\,M_{\odot}$.

###  3.4 Models reaching the WR stage

Out of the 94 models computed to $T_{c}\geq 10^{9}\,\mathrm{K}$, only 8 reach
the conditions to switch to the WR wind scheme. These are all $35\,M_{\odot}$
models computed with $\eta=1.0$, and none use the vL algorithm in the cool
phase that precedes the WR phase (cf. Tab. 4). The lack of WR models using vL
is explained by the self-damping of this mass loss scheme for more massive and
thus more luminous stars (see Sec. 3.3).

The typical duration of the WR phase is $\sim 0.02-0.05\,\mathrm{Myr}$. The
differences in duration of the WR phase with the Nugis & Lamers (2000) (NL)
and Hamann et al. (1982, 1995); Hamann & Koesterke (1998) (H) algorithms are
negligible. However, models computed with the dJ algorithm in the cool phase
have WR phases that are systematically longer by a few ten thousand years than
the corresponding models computed with the NJ algorithm. Moreover, the
duration of the WR phase of models computed with the K algorithm in the hot
phase is about a factor of $\sim 2.3$ longer than in models using the V
algorithm. This is because $M_{\mathrm{ZAMS}}=35\,M_{\odot}$ models computed
with $\eta=1.0$ and the V algorithm reach the end of the hot phase with a
helium core that is $0.22\,M_{\odot}$ less massive than models computed with
the K algorithm (cf. Tab. 6). Therefore, the subsequent evolution is slowed
down, and the WR phase is reached slightly later.

The NL WR phase algorithm produces higher final masses than the H algorithm:
the difference is $\sim 0.3$ $(\sim 0.8)\,M_{\odot}$ for models using the V
(K) hot phase algorithm and does not depend strongly on cool phase mass loss
(cf. Tab. 3). Note that these differences are very small fractions of the
initial mass $M_{\mathrm{ZAMS}}=35\,M_{\odot}$ of these models (cf. Fig. 1).

The WR wind does not have a strong influence on $M_{\mathrm{He}}$ at the end
of the mass loss phase. For stars with $M_{\mathrm{ZAMS}}\lesssim
40M_{\odot}$, the He core mass is determined well before the beginning of the
WR phase, and the wind mass loss is not strong enough to dig into the He core
directly. More massive stars with stronger winds, may become hydrogen depleted
already during the core hydrogen burning phase. In that case, the WR mass loss
rate might have an impact on $M_{\mathrm{He}}$ through the quasi-static
response of the convective core to mass loss (Meynet et al., 1994; de Koter et
al., 1997; Crowther et al., 2010; Bestenlehner et al., 2011). Also, in even
more massive stars, the entire hydrogen-rich envelope can be lost to winds,
making $M_{\mathrm{He}}$ the total mass of the star, and winds can then
further reduce it (e.g. Woosley, 2016).

While we find no systematic effect of WR mass loss on $M_{\mathrm{He}}$, the
situation is different for $M_{\mathrm{CO}}$ (cf. Tab. 3). NL models yield
$M_{\mathrm{CO}}$ that are systematically lower by $\sim 0.05-0.1\,M_{\odot}$
than H models and the largest differences are between models that also use
different hot phase mass loss algorithms. We find that models differing only
in the WR algorithm have lower $M_{\mathrm{CO}}$ for higher final masses. For
example, the V-dJ-NL model has $M_{\mathrm{CO}}=9.81\,M_{\odot}$ and final
mass $M=20.03\,M_{\odot}$. In contrast, the V-dJ-H model has a higher
$M_{\mathrm{CO}}=9.93\,M_{\odot}$, but a lower final mass of
$M=19.73\,M_{\odot}$. However, we note that the trend that lower final masses
correspond to higher $M_{\mathrm{CO}}$ holds for most of our $35\,M_{\odot}$
models, independent of if they become WR stars or not.

The differences in both total mass and $M_{\mathrm{CO}}$ found varying the WR
mass loss algorithm are more sensitive to the previously employed hot phase
mass loss algorithm than to the cool phase algorithm. Moreover,
$M_{\mathrm{He}}$ is also almost insensitive to the cool phase (cf. Sec. 3.3)
and WR phase mass loss. Therefore, the differences found here are most likely
related to the differences in $M_{\mathrm{He}}$ at the end of the hot phase
(and in the total mass for a given $M_{\mathrm{He}}$), and consequently the
position of the hydrogen burning shell, which indirectly influences the helium
burning, the mixing processes shaping the composition profile, and ultimately
the resulting $M_{\mathrm{CO}}$ and the amount of mass lost during the WR
phase.

###  3.5 Models at oxygen depletion

Table 7: Stellar properties at core oxygen depletion, i.e. $X_{c}(\mathrm{{}^{16}O})<0.04$ and $X_{c}(\mathrm{{}^{28}Si})>0.01$. The last column shows the maximum difference in the compactness parameter $\xi_{2.5}^{\mathrm{O\ depl}}$ for each choice of $\eta$ and $M_{\mathrm{ZAMS}}$. These runs are re-started from the corresponding MESA models saved at the end of the mass loss phase, when $T_{c}\geq 10^{9}\,\mathrm{K}$.  Oxygen depletion: $X_{c}(\mathrm{{}^{16}O})<0.04$ and $X_{c}(\mathrm{{}^{28}Si})>0.01$ |
---|---
$M_{\mathrm{ZAMS}}\ [M_{\odot}]$ | $\eta$ | ID | $R\ [R_{\odot}]$ | $M_{\mathrm{tot}}\ [M_{\odot}]$ | $M_{\mathrm{He}}\ [M_{\odot}]$ | $M_{\mathrm{CO}}\ [M_{\odot}]$ | $\xi_{2.5}^{\mathrm{O\ depl}}$ | $\max\Delta\xi_{2.5}^{\mathrm{O\ depl}}$
15 | 0.13 | V-dJ |  0908 | 14.66 |  04.99 | 3.17 | 0.155 | 0.001
V-vL |  0924 | 13.60 |  04.98 | 3.17 | 0.156
K-NJ |  0911 | 14.66 |  05.00 | 3.17 | 0.156
0.33 | V-NJ |  0914 | 13.90 |  04.93 | 3.13 | 0.153 | 0.001
V-vL |  0967 | 10.39 |  04.93 | 3.13 | 0.152
1.03 | V-NJ |  0895 | 12.74 |  04.65 | 2.95 | 0.141 | 0.014
V-vL |  0629 |  05.25 |  04.64 | 2.94 | 0.139
K-NJ |  0946 | 11.87 |  04.92 | 3.12 | 0.153
K-vL |  0643 |  05.70 |  04.90 | 3.11 | 0.152
20 | 0.13 | V-dJ |  0994 | 19.23 |  07.04 | 4.83 | 0.182 | 0.007
V-vL |  0998 | 18.10 |  07.03 | 4.83 | 0.178
K-NJ |  0991 | 19.37 |  07.01 | 4.82 | 0.179
K-vL |  0991 | 18.67 |  07.01 | 4.82 | 0.175
0.33 | V-NJ | 1001 | 17.48 |  07.01 | 4.81 | 0.161 | 0.052
V-vL |  0987 | 13.47 |  06.99 | 4.80 | 0.175
K-dJ |  0999 | 17.62 |  07.00 | 4.80 | 0.213
K-vL |  0993 | 13.90 |  06.98 | 4.80 | 0.175
1.03 | V-dJ |  0951 | 11.81 |  07.06 | 4.85 | 0.182 | 0.010
V-vL |  0673 |  08.80 |  07.02 | 4.82 | 0.176
K-NJ |  0978 | 12.77 |  06.95 | 4.77 | 0.176
K-vL |  0712 |  08.81 |  06.90 | 4.73 | 0.172
25 | 0.13 | V-dJ |  0898 | 23.85 |  09.14 | 6.46 | 0.180 | 0.023
V-vL |  0899 | 23.56 |  09.14 | 6.44 | 0.179
K-NJ |  0888 | 23.76 |  09.24 | 6.52 | 0.159
K-vL |  0889 | 23.24 |  09.23 | 6.50 | 0.157
0.33 | V-NJ |  0922 | 22.02 |  08.87 | 6.25 | 0.161 | 0.049
V-vL |  0929 | 22.13 |  08.87 | 6.20 | 0.164
K-dJ |  0894 | 21.27 |  09.11 | 6.42 | 0.210
1.03 | V-NJ |  0812 | 13.67 |  09.05 | 6.32 | 0.164 | 0.050
V-vL |  0786 | 12.97 |  09.04 | 6.29 | 0.211
K-dJ |  0875 | 15.57 |  08.89 | 6.21 | 0.161
K-vL |  0860 | 14.69 |  08.90 | 6.23 | 0.200
30 | 0.13 | V-dJ |  0701 | 28.11 | 10.97 | 7.89 | 0.242 | 0.001
V-NJ |  0697 | 28.17 | 10.97 | 7.90 | 0.243
K-NJ |  0706 | 28.19 | 10.91 | 7.83 | 0.242
K-vL |  0719 | 28.75 | 10.90 | 7.74 | 0.242
0.33 | V-dJ |  0668 | 24.16 | 11.11 | 7.99 | 0.243 | 0.023
V-vL |  0723 | 26.02 | 11.10 | 8.00 | 0.243
K-NJ |  0716 | 25.18 | 10.89 | 7.85 | 0.228
K-vL |  0709 | 27.12 | 10.87 | 7.79 | 0.220
1.03 | V-dJ |  0610 | 15.36 | 11.28 | 8.20 | 0.179 | 0.065
K-dJ |  0756 | 18.51 | 10.92 | 7.85 | 0.244
K-vL |  0723 | 22.53 | 10.89 | 7.74 | 0.243

To reduce the computational cost of our model grid, we select a subset of 44
stars to continue until oxygen depletion, which we define as the time when
$X_{c}(^{16}\mathrm{O})\leq 0.04$ and $X_{c}(^{28}\mathrm{Si})\geq 0.01$.
These models span the range of properties found at the end of the mass loss
phase, and are listed in Tab. 7 together with their properties at oxygen
depletion. This selection allows us to avoid running multiple models that have
very similar evolutionary paths from the end of the mass loss phase onward.
The duration of the evolution between the end of the mass loss phase
($T_{c}\geq 10^{9}\,\mathrm{K}$) and oxygen depletion is of order years to
decades, depending on the total mass and core masses at the end of the mass
loss phase.

In the very short time to oxygen depletion, neither total mass nor helium core
mass change appreciably. This can be inferred by comparing the entries of Tab.
3 and Tab. 7, which also reveals that the stellar radii vary only within $\pm
3\,R_{\odot}$ in most models.

The CO core masses at oxygen depletion summarized in Tab. 7 are systematically
a few percent lower than those listed in Tab. 3 at the end of the mass loss
phase. This seems counter-intuitive, since one would expect the CO core to
grow in mass because of the ashes of helium shell burning. However, the
boundary location for the CO core is determined by convective mixing within
and above the He burning shell, which brings helium-rich material inward and
moves the CO core boundary to a smaller mass coordinate. _This implies that
the core mass is very sensitive to the mixing parameters._ The maximum spread
in $M_{\mathrm{CO}}$ obtained varying the wind algorithm is of order $\sim
0.1\,M_{\odot}$ and it increases with $M_{\mathrm{ZAMS}}$ up to about
$\max(\Delta M_{\mathrm{CO}})\simeq 0.5\,M_{\odot}$ for $30\,M_{\odot}$ models
(see also Tab. 5). Note that the $15\,M_{\odot}$ models are outliers in that
they show a larger spread of CO core masses of up to $0.28\,M_{\odot}$ between
models using the V and K hot mass loss algorithm. This is a consequence of a
combination of (_i_) the V algorithm leading to a lower total mass and a lower
He core mass at the end of the hot phase and (_ii_) the relatively low mass
loss during the cool phase (compared to more massive stars), which results in
deeper convective episodes. Together, these lead to small CO cores: the two
$15M_{\odot}$ models using the V algorithm with $\eta=1.0$ have
$M_{\mathrm{CO}}\simeq 2.95M_{\odot}$ at the end of the mass loss phase, which
is about $0.2\,M_{\odot}$ smaller than the average for $15\,M_{\odot}$ models.

###  3.6 The Compactness Parameter $\xi_{\mathcal{M}}$

Although the internal structure of our models is not yet final at oxygen
depletion, some quantities (e.g., core masses) are already close to their
final values, and it becomes possible to draw first connections between
internal structure and the potential final outcome of core collapse (Sukhbold
& Woosley, 2014). For this, we include the compactness parameter
$\xi_{2.5}^{\mathrm{O\,depl}}$ in Tab. 7. O’Connor & Ott (2011) define the
compactness parameter as

$\xi_{\mathcal{M}}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\frac{\mathcal{M}/M_{\odot}}{R(\mathcal{M})/1000\
\mathrm{km}}\ \ .$ (3)

This parameter provides a single measure for the complex inner core structure
(mass coordinate smaller than $\mathcal{M}$) of a star, allowing for a
simplified discussion of the differences in the internal structure produced by
the various mass loss algorithms we compare.

We set $\mathcal{M}=2.5\,M_{\odot}$ because this is the typical mass above
which the proto-NS that will form during core collapse will become a BH
(O’Connor & Ott, 2011). This mass cut remains well outside of the typical iron
core mass and includes the layers of the star that the shock will encounter
after core bounce. These layers determine the accretion ram pressure that the
shock has to overcome for a successful explosion. The qualitative behavior of
the supernova dynamics and outcome with $\xi_{\mathcal{M}}$ is known to be
robust against different choices of $\mathcal{M}$ (O’Connor & Ott, 2011, 2013;
Ugliano et al., 2012; Sukhbold & Woosley, 2014).

One-dimensional parametric core-collapse SN explosion simulations (O’Connor &
Ott, 2011, 2013; Ugliano et al., 2012; Ertl et al., 2016; Sukhbold et al.,
2016) show that the value of $\xi_{2.5}$ at the onset of core collapse
indicates the most probable remnant. High values of $\xi_{2.5}$ indicate a
more compact pre-SN structure that is harder to explode and that will more
likely result in a BH remnant. Conversely, low values of $\xi_{2.5}$ indicate
a steeper density gradient and an easier to explode structure, suggesting that
the remnant will more likely be a NS (O’Connor & Ott, 2011; Ugliano et al.,
2012; Clausen et al., 2015). Sukhbold & Woosley (2014) suggest that the value
of the compactness parameter evaluated at oxygen depletion,
$\xi_{2.5}^{\mathrm{O\ depl}}$, can already be used to infer the most likely
outcome of the core collapse event. The evolution from oxygen depletion to
core collapse tends to increase the compactness and amplify the differences
between different stellar models, but the key features that determine the
interpretation of $\xi_{2.5}$ appear to be set already at oxygen depletion.
Other parameters to relate the pre-SN structure to the most-likely remnant can
be defined in the context of neutrino-driven explosions (see Ertl et al.,
2016), but they rely on physical quantities that are not at all set at earlier
stages of the evolution (e.g., the entropy profile throughout the silicon
layer and iron core), and therefore they are not useful diagnostics before the
onset of core collapse.

####  3.6.1 Evolution of $\xi_{2.5}$ until oxygen depletion

Figure 4: Top panel: Time evolution of the compactness parameter $\xi_{2.5}$
for all $M_{\mathrm{ZAMS}}=25\,M_{\odot}$ models computed until oxygen
depletion. Bottom panel: central temperature and density for the same models
from when the central density increases above $10^{6}\,\mathrm{g\ cm^{-3}}$.
All curves end at oxygen depletion ($X_{c}(^{16}\mathrm{O})\leq 0.04$ and
$X_{c}(^{28}\mathrm{Si})\geq 0.01$). The color indicates the wind algorithm
combination, labeled according to Tab. 2. Dot-dashed, dashed, and solid curves
are calculated using $\eta=0.1,\,0.33,\,1.0$, respectively. Note that we only
show models listed in Tab. 7. The vertical dot-dashed lines in the top panel
approximately indicate core hydrogen exhaustion (TAMS), core helium exhaustion
($Y_{c}\simeq 0$), and end of the mass loss phase ($T_{c}\geq
10^{9}\,\mathrm{K}$).

The compactness parameter is a function of time,
$\xi_{\mathcal{M}}\equiv\xi_{\mathcal{M}}(t)$, because of the changes in the
radius of a given mass coordinate $\mathcal{M}$. These can be caused by
contraction of the core, onset of partial electron degeneracy within the mass
coordinate $\mathcal{M}$, and by episodes of convective mixing and shell
burning (Sukhbold & Woosley, 2014). The top panel of Fig. 4 shows examples of
the evolution of $\xi_{2.5}$ until oxygen depletion in our $25\,M_{\odot}$
models. Note that we use a reversed logarithmic scale on the x axis to
emphasize the late evolutionary stages.

Figure 4 shows that the compactness parameter is constant during the main
sequence evolution. During this phase, it is also almost independent of
$M_{\mathrm{ZAMS}}$ and mass loss algorithm because all stars considered here
have convective main-sequence cores that are always much larger than the mass
coordinate at which we evaluate the compactness. After core hydrogen
exhaustion, $\xi_{2.5}$ increases because of the overall contraction, reaching
$\xi_{2.5}\simeq 0.02$ in our $25\,M_{\odot}$ models. Then it slowly continues
to increase during the hydrogen-shell burning and helium core burning phases.
The increase speeds up significantly during core carbon burning, reaching
values of $\xi_{2.5}\simeq 0.1$ in our $25\,M_{\odot}$ models. Neon core
burning ignition and the onset of carbon shell burning mark a critical point
in the evolution of $\xi_{2.5}$ at which the various curves in Fig. 4 begin to
diverge. Sukhbold & Woosley (2014) find the same and point out that the
subsequent evolution of $\xi_{2.5}$ is highly sensitive to the details of
carbon shell burning (i.e. the number, locations, and durations of shell
burning episodes). It is important to note from Fig. 4 that the effects of
mass loss on core structure (represented by $\xi_{2.5}$) are delayed: At the
time mass loss ends in our models (at $T_{c}>10^{9}\,\mathrm{K}$), differences
in $\xi_{2.5}$ are minute. These seed differences grow and become substantial
only in the last decade before core collapse.

The bottom panel of Fig. 4 depicts central density–temperature tracks for our
$25\,M_{\odot}$ models that are evolved to oxygen depletion. The tracks start
roughly at neon core ignition and show that the mass-loss history (i.e. the
choice of wind mass loss algorithm combination) also influences the innermost
core thermodynamics and structure. This is because the nuclear burning
processes in the core are regulated by the amount of mass that needs to be
sustained by the core itself, i.e. the mass below the innermost burning shell
above the core. This in turn depends on the location and luminosity of the
shell burning regions and therefore on the total mass of the star.

In Fig. 5, we show the values of $\xi_{2.5}^{\mathrm{O\ depl}}$ for all models
that we run to oxygen depletion. The spread in each panel is due to the
different algorithmic treatments of wind mass loss: for example,
$25\,M_{\odot}$ models show values ranging between 0.210 and 0.157. Generally
speaking, the spread in $\xi_{2.5}^{\mathrm{O\ depl}}$ increases with
increasing $\eta$ and $M_{\mathrm{ZAMS}}$, i.e. the stronger the stellar wind,
the more it influences the core structure. We emphasize that a few percent
variation of $\xi_{2.5}^{\mathrm{O\ depl}}$ can result in important
differences in the core structure at the onset of core collapse: the
subsequent contraction of the core and the details of carbon, oxygen, and
silicon shell burning, amplify the differences between models that are still
relatively similar at oxygen depletion (see Sec. 3.7 and Sukhbold & Woosley,
2014).

Figure 5: Compactness parameter at oxygen depletion
$X_{c}(^{16}\mathrm{O})\leq 0.04$ and $X_{c}(^{28}\mathrm{Si})\leq 0.01$. Each
color corresponds to a wind algorithm combination (top axis), and each panel
shows a different initial mass, indicated on the bottom. Crosses, pluses, and
circles correspond to $\eta=0.1,\,0.33,\,1.0$, respectively. Note that we show
only the models listed in Tab. 7. The vertical spread indicates different core
structures, which will evolve differently to the onset of core collapse, and
possibly result in different SN outcomes.

####  3.6.2 Effects of the wind efficiency on $\xi_{2.5}^{\mathrm{O\ depl}}$

The effects of varying $\eta$ on $\xi_{2.5}^{\mathrm{O\ depl}}$ can be
inferred from the comparison of models in Tab. 7 of the same
$M_{\mathrm{ZAMS}}$ using the same wind algorithm combinations but different
efficiencies. The variations of $\xi_{2.5}^{\mathrm{O\ depl}}$ with $\eta$ are
typically non-monotonic: for example, the $25\,M_{\odot}$ model computed using
the V-vL combination reaches oxygen depletion with $\xi_{2.5}^{\mathrm{O\
depl}}=0.179$, 0.164, and 0.211 for $\eta=0.1$, 0.33, and 1.0, respectively.
Within the framework of each wind algorithm, higher values of $\eta$
correspond to higher mass loss rates and thus to a progressive shift of the
evolution toward that of lower initial mass. However, the compactness
parameter is known to be a highly non-monotonic function of
$M_{\mathrm{ZAMS}}$ (Sukhbold & Woosley, 2014). Therefore, a higher mass loss
rate can sometimes result in a decrease, an increase, or even almost no
variation of $\xi_{2.5}^{\mathrm{O\ depl}}$. For example, in $15\,M_{\odot}$
models computed with V-vL, $\xi_{2.5}^{\mathrm{O\ depl}}$ decreases when going
from $\eta=0.1$ to $\eta=1.0$. In $25\,M_{\odot}$ models computed with V-vL,
we instead find that $\xi_{2.5}^{\mathrm{O\ depl}}$ increases with the same
change in $\eta$. And in the case of $20\,M_{\odot}$ models computed with
V-dJ, we find only tiny variations of $\xi_{2.5}^{\mathrm{O\ depl}}$ when
$\eta$ is varied from $0.1$ to $1.0$. See Tab. 7 for more examples and
details.

####  3.6.3 Effects of varying the mass loss algorithm on
$\xi_{2.5}^{\mathrm{O\ depl}}$

The last column of Tab. 7 shows that the effects of different wind mass loss
algorithms on $\xi_{2.5}$ are in most cases small until oxygen depletion for
$\eta=0.1$ and $\eta=0.33$. For $\eta=0.1$, the V algorithm generally results
in higher values of $\xi_{2.5}^{\mathrm{O\ depl}}$ than the K algorithm. This
holds for all studied ZAMS masses with the exception of the $30\,M_{\odot}$
models. These do not exhibit this trend because all their burning shells are
outside the mass coordinate $\mathcal{M}=2.5\,M_{\odot}$, as can be seen from
the values of $M_{\mathrm{CO}}$ listed in Tab. 7.

The spread in $\xi_{2.5}^{\mathrm{O\ depl}}$ between different mass loss
algorithms increases for models with $\eta=0.33$. For example, Fig. 5 shows
that the $25\,M_{\odot}$ model with the K-dJ combination reaches
$\xi_{2.5}^{\mathrm{O\ depl}}\simeq 0.21$ (similar to its $20\,M_{\odot}$
counterpart, cf. Tab. 7), which is $\sim 30\%$ higher than with other
algorithm combinations.

Models computed with $\eta=1.0$ are most suitable to discuss the effect of
different wind mass loss algorithm combinations. Both the hot phase and the
cool phase mass loss algorithms influence $\xi_{2.5}^{\mathrm{O\ depl}}$, but
their detailed effect varies with $M_{\mathrm{ZAMS}}$ and is strongest in the
$25\,M_{\odot}$ and $30\,M_{\odot}$ models. For example, from Tab. 7, we find
$\xi_{2.5}^{\mathrm{O\ depl}}\sim 0.16$, for the $\eta=1.0$, $25\,M_{\odot}$
models with the V-NJ and K-dJ combinations, while models with V-vL and K-vL
result in $\xi_{2.5}^{\mathrm{O\ depl}}\sim 0.2$. At lower
$M_{\mathrm{ZAMS}}$, even for $\eta=1.0$, differences in
$\xi_{2.5}^{\mathrm{O\ depl}}$ due to the choice of mass loss algorithm
combination are overall (with few exceptions) rather small and typically at
the few percent level. While these differences will be amplified by the
subsequent evolution toward core collapse, they are small compared to the
tremendous differences in total mass resulting from the different algorithm
combinations (cf. Tab. 7).

An interesting question to address is the relative importance of hot phase
(i.e. main sequence) and cool phase (i.e. post main sequence) mass loss for
$\xi_{2.5}^{\mathrm{O\ depl}}$. Naively, one would think that by the time the
core and envelope are essentially decoupled, loss of envelope mass should have
limited impact on the subsequent evolution of the core. Our results suggest
that this is not generally the case.

The limited set of models run to oxygen depletion and listed in Tab. 7 give a
complex, but necessarily incomplete picture of the relative importance of each
mass loss phase for the core structure. Which phase is most relevant depends
on $M_{\mathrm{ZAMS}}$ and $\eta$. For brevity, we focus here on the $\eta=1$
case and compare models evolved with the same hot phase mass loss algorithm (V
or K) and different cool phase algorithms.

For $15\,M_{\odot}$ models with $\eta=1$, the tremendous mass loss with the vL
algorithm in the cool phase has little effect on He and CO core masses and on
$\xi_{2.5}^{\mathrm{O\ depl}}$. For example, the final masses of K-vL and K-NJ
are $5.70\,M_{\odot}$ and $11.87\,M_{\odot}$, respectively. Yet their
$\xi_{2.5}^{\mathrm{O\ depl}}$ are very close to each other, $0.152$ and
$0.153$, respectively. Qualitatively, the same is true for the V-NJ and V-vL
combinations. On the other hand there is a larger spread between combinations
using V and K, suggesting that the small differences seeded by hot phase mass
loss dominate in the $15\,M_{\odot}$ $\eta=1$ case. Similarly, we find for
$30\,M_{\odot}$ models with $\eta=1$ that K-dJ and K-vL lead to final masses
of $18.51\,M_{\odot}$ and $22.53\,M_{\odot}$, respectively, but their
$\xi_{2.5}^{\mathrm{O\ depl}}$ are $0.244$ and $0.243$. The situation is more
complicated for $25\,M_{\odot}$ and $20\,M_{\odot}$ models that straddle the
$M_{\mathrm{ZAMS}}$ range where $\xi_{2.5}$ varies chaotically (Sukhbold &
Woosley, 2014). From Tab. 7 we see that for the $20\,M_{\odot}$, $\eta=1$
models all considered K- and V- combinations yield roughly the same
$\xi_{2.5}^{\mathrm{O\ depl}}$. In $25\,M_{\odot}$, $\eta=1$ models, on the
other hand, cool phase mass loss has the dominant impact on
$\xi_{2.5}^{\mathrm{O\ depl}}$. For example, K-dJ and K-vL have compactness of
$0.161$ and $0.200$, respectively, although their final masses differ by only
$\sim 1\,M_{\odot}$ due to the self-regulation of the vL algorithm in this
mass range.

####  3.6.4 Comparison with Sukhbold & Woosley (2014)

Sukhbold & Woosley (2014) employed the NJ mass loss algorithm without any
efficiency scaling factor (i.e. at efficiency $\eta=1.0$) throughout the
entire evolution of their models. Comparing our Tab. 7 with their Fig. 23, we
find that the compactness parameter values at oxygen depletion of our models
lie in the same range as theirs, with a tendency toward slightly higher
values.

Our $15\,M_{\odot}$ models produce values of $\xi_{2.5}^{\mathrm{O\
depl}}\simeq 0.15$, which is slightly higher than their value of $\sim
0.11-0.13$, especially for reduced wind mass loss rates (i.e. $\eta<1.0$). For
this initial mass, increasing $\eta$ decreases the compactness of the core and
reduces the difference between our models and those of Sukhbold & Woosley
(2014).

Most of our $20\,M_{\odot}$ models have $\xi_{2.5}^{\mathrm{O\ depl}}\simeq
0.18$, close to the corresponding models of Sukhbold & Woosley (2014). For
these models, the maximum difference varying the wind mass loss algorithm is
$\Delta\xi_{2.5}^{\mathrm{O\ depl}}\lesssim 0.05$. This is not surprising,
because large variations of $\xi_{2.5}$ are expected because of the transition
from convective to neutrino-cooled and radiative carbon shell burning, which
happens around $M_{\mathrm{ZAMS}}\simeq 20\,M_{\odot}$ (Sukhbold & Woosley,
2014). Therefore, in this mass range, changing the wind mass loss algorithm
can substantially change $\xi_{2.5}^{\mathrm{O\ depl}}$ by shifting the
evolutionary track of the star in the slightest way.

Our $25\,M_{\odot}$ models have values of $\xi_{2.5}^{\mathrm{O\ depl}}$
similar to those of the $20\,M_{\odot}$ models, in agreement with Sukhbold &
Woosley (2014). However, once again, we obtain a large variation of
$\xi_{2.5}^{\mathrm{O\ depl}}$ changing the treatment of mass loss
($\Delta\xi_{2.5}^{\mathrm{O\ depl}}\lesssim 0.05$).

Most of our $30\,M_{\odot}$ models have $\xi_{2.5}^{\mathrm{O\ depl}}\simeq
0.23$, which is significantly larger than the corresponding value of $\sim
0.16$ found by Sukhbold & Woosley (2014). The variations of the compactness
parameter with mass loss algorithm combination and efficiency are smaller for
the $30M_{\odot}$ models than for lower $M_{\mathrm{ZAMS}}$, except for the
model computed with $\eta=1.0$ and the V-dJ algorithm. This model has
$\xi_{2.5}^{\mathrm{O\ depl}}=0.179$, which is much closer to the values of
Sukhbold & Woosley (2014). We expect that the similar V-NJ algorithm
combination with $\eta=1.0$ would produce a structure close to this model at
oxygen depletion, based on the similarities between the models at the end of
the mass loss phase. The relatively low compactness of this model is therefore
likely a consequence of the V hot phase mass loss algorithm. For
$M_{\mathrm{ZAMS}}\gtrsim 30\,M_{\odot}$, V with $\eta=1.0$ produces
substantially more mass loss than K (cf. Tab. 6), and therefore, the
subsequent evolution is closer to the path of less massive stars – which are
also expected to reach oxygen depletion with a lower compactness.

We speculate that the quantitative differences between our findings for the
compactness parameter at oxygen depletion and those of Sukhbold & Woosley
(2014) are in fact due primarily to their choice of mass loss algorithm: they
employ the NJ algorithm in both the hot and the cold phase, which results in
overall greater early mass loss than the K algorithm (and the V algorithm
until the bistability jump). This notion is corroborated by our finding that
our $30\,M_{\odot}$ model closest to theirs is the one that loses the most
mass during the hot phase (using $\eta=1$ and the V-dJ combination; cf. Tab.
6).

###  3.7 Models at the onset of core collapse

We select a subset of six of our models at oxygen depletion for continuation
to the onset of core collapse. Reducing the model set is necessary to limit
the computational cost of this study. We choose two $15\,M_{\odot}$ models
with efficiency $\eta=1.0$ and mass loss algorithm combination V-NJ and K-VL,
two $20\,M_{\odot}$ models with $\eta=0.33$ computed using the V-vL and the
K-dJ combination, and two $30\,M_{\odot}$ models with $\eta=0.33$ and the V-dJ
or the K-NJ combination. When restarting our models from oxygen depletion, we
switch from the 45-isotopes nuclear reaction network used so far to a larger
customized network with 203-isotopes (see Appendix B). This is necessary to
capture core deleptonization due to electron capture during and after silicon
burning. Continuing the evolution with a larger nuclear reaction network also
requires reducing the number of computational mesh points (from $\sim 10^{4}$
to $\sim 10^{3}$) to run the simulations within the memory constraints of MESA
(see Appendix B for details). By the time oxygen depletion is reached, the
effect of wind mass loss on the core structure is already pronounced, and our
reduced-resolution models still have spatial resolution that is comparable to
that of published models (e.g., Woosley et al., 2002; Woosley & Heger, 2007).
Furthermore, resolution tests in Appendix B.1.3 give us confidence that the
presently unavoidable reduction in resolution does not affect our overall
results and conclusions.

Figure 6 shows the time evolution of $\xi_{2.5}$ (top panel) and the central
temperature–central density evolutionary tracks (bottom panel) from oxygen
depletion to the onset of core collapse for our pre-SN model set. The top
panel shows that $\xi_{2.5}$ settles onto its final value before the criterion
for the onset of core collapse is reached. The bottom panel clearly shows a
hook at $\log_{10}(T_{c}/\mathrm{[K]})\simeq 9.55$, where $\rho_{c}$ decreases
at roughly constant temperature, indicating the point of silicon core
ignition. This is also the more “noisy” part of these tracks, indicating that
this phase of nuclear burning with very high and nearly balancing reaction
rates is the most challenging to simulate (Hix & Thielemann, 1996; Hix et al.,
2007).

Table 8: Properties of the subset of models run to the onset of core collapse,
defined as the time when $\mathrm{max}\\{|v|\\}\geq 10^{3}\ \mathrm{km\
s^{-1}}$. $M_{4}$ and $\mu_{4}$ are the parameters used to predict the SN
outcome of a stellar model in Ertl et al. (2016), see also text.
$M_{\rho_{6}}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}M(\rho=10^{6}\mathrm{g\
cm^{-3})}$ is the mass enclosed in the location where the density drops below
$10^{6}\ \mathrm{g\ cm^{-3}}$, $M_{\mathrm{CO}}$ and $M_{\mathrm{Fe}}$ are the
carbon-oxygen and iron core masses, respectively.  Onset of core collapse:
$\mathrm{max}\\{|v|\\}\geq 10^{3}\ \mathrm{km\ s^{-1}}$
---
$M_{\mathrm{ZAMS}}\ [M_{\odot}]$ | $\eta$ | ID | $\xi_{2.5}^{\mathrm{pre-SN}}$ | $M_{4}\ [M_{\odot}]$ | $\mu_{4}$ | $M_{\rho_{6}}\ [M_{\odot}]$ | $M_{\mathrm{CO}}\ [M_{\odot}]$ | $M_{\mathrm{Fe}}\ [M_{\odot}]$
15 | 1.0 | V-NJ | 0.103 | 1.71 | 0.045 | 1.68 | 2.91 | 1.39
K-vL | 0.132 | 1.78 | 0.051 | 1.79 | 3.07 | 1.50
25 | 0.33 | V-vL | 0.227 | 1.73 | 0.084 | 1.84 | 6.38 | 1.51
K-dJ | 0.308 | 2.05 | 0.100 | 2.19 | 6.40 | 1.63
30 | 0.33 | V-dJ | 0.358 | 1.60 | 0.163 | 2.21 | 7.98 | 1.56
K-NJ | 0.276 | 1.82 | 0.100 | 1.98 | 7.90 | 1.58

Figure 6: Time evolution of the compactness parameter $\xi_{2.5}$ (top panel)
and evolution in $\rho_{c}-T_{c}$ space (bottom panel) from oxygen depletion
to the onset of core collapse. Solid curves correspond to $15\,M_{\odot}$
models computed with $\eta=1.0$, dashed and dot-dashed curves correspond to
$25\,M_{\odot}$ and $30\,M_{\odot}$ models, respectively, computed with
$\eta=0.33$. The vertical dot-dashed line in the top panel indicates roughly
the time of core silicon depletion ($X(^{28}\mathrm{Si})\leq 0.01$).

Figure 7: Neutrino luminosity from oxygen depletion to core collapse for the
$15\,M_{\odot}$ model run with the V-NJ scheme and $\eta=1.0$. This shows that
variations in the nuclear burning rate cause the oscillations in $\xi_{2.5}$.
The red curve shows only neutrinos coming from nuclear reactions, the cyan
curve shows thermal neutrinos, and the blue curve is the sum of the two. The
dotted line shows the compactness parameter (right hand side vertical axis).
The vertical dot-dashed line indicates the approximate time of core silicon
depletion. The $\xi_{2.5}$ oscillations are driven by variations in silicon
and oxygen shell burning.

The evolution of the compactness parameter shows maxima and minima, which are
related to oxygen shell ignition, around $\log_{10}(\\{t_{\mathrm{pre-
SN}}-t\\}/\mathrm{[yr]})\simeq-2$, silicon core ignition around
$\log_{10}(\\{t_{\mathrm{pre-SN}}-t\\}/\mathrm{[yr]})\simeq-3$, and silicon
shell ignition at about $\log_{10}(\\{t_{\mathrm{pre-
SN}}-t\\}/\mathrm{[yr]})\simeq-4.5$. However, note that the ignition times and
the durations of these burning phases are mass-dependent. Figure 7 shows the
corresponding increase in the neutrino emission from nuclear reactions for the
$15\,M_{\odot}$ model computed with the V-NJ combination and $\eta=1.0$.
Similar features are present in the $T_{c}(\rho_{c})$ evolutionary tracks.
However, while the $T_{c}(\rho_{c})$ track only probes the innermost part of
the stellar core, the $\xi_{2.5}$ evolution is determined by the interplay
between silicon, oxygen, and carbon burning shells, core contraction, and
onset of electron degeneracy.

Interestingly, the final compactness of the $15\,M_{\odot}$ models is lower at
the onset of core collapse than at oxygen depletion: for example, the
$15\,M_{\odot}$ model computed with K-VL and $\eta=1.0$ has
$\xi_{2.5}^{\mathrm{pre-SN}}=0.132<\xi_{2.5}^{\mathrm{O\ depl}}=0.152$. This
is because of the presence of nuclear burning shells within
$\mathcal{M}=2.5\,M_{\odot}$ whose energy generation tends to expand the
material in the layers above them (the same process happens during the
Hertzsprung gap for hydrogen-shell-burning stars). The location of the shells
can be estimated using the core masses listed in Tab. 8 (and Tab. 7): the
$15\,M_{\odot}$ models have two shells of nuclear burning (Si and O,
respectively) within the $2.5\,M_{\odot}$ mass coordinate, while only one
shell exists in this region at oxygen depletion. Models with
$M_{\mathrm{ZAMS}}>15\,M_{\odot}$ settle on a pre-SN compactness that is
higher than the corresponding $\xi_{2.5}^{\mathrm{O\ depl}}$, because only the
silicon burning shell is within $\mathcal{M}=2.5\,M_{\odot}$.

Table 8 lists the properties of the six models that we run to the onset of
core collapse. $\xi_{2.5}^{\mathrm{pre-SN}}$, and ($M_{4}$, $\mu_{4}$), where
$M_{4}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}M(s=4)$ is the mass location
where the specific entropy is $s=4\,\mathrm{k_{b}\,baryon^{-1}}$ and
$\mu_{4}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}dm/dr|_{s=4}$ is the mass
gradient at that location, offer two different ways to estimate how hard it
will be for the SN shock to unbind the stellar mantle and leave a NS remnant
(see Ertl et al., 2016). The total mass at density higher than
$10^{6}\,\mathrm{g\ cm^{-3}}$ ($M_{\rho 6}$), the carbon-oxygen core mass
($M_{\mathrm{CO}}$), and the iron core mass ($M_{\mathrm{Fe}}$), defined as
the location where $X(^{28}\mathrm{Si})<0.01$) can be used to estimate the
nickel yields of the possible SN explosion and the remnant mass (Fryer et al.,
2012; Sukhbold et al., 2016).

The final $M_{\mathrm{CO}}$ and $M_{\mathrm{Fe}}$ depend, although weakly, on
the mass loss algorithm adopted during the hot evolutionary phase (V or K).
The V algorithm yields slightly smaller cores (and total masses at the end of
the hot phase, cf. Sec. 3.2). The $15\,M_{\odot}$ models with $\eta=1.0$ reach
core collapse with $M_{\mathrm{CO}}=2.91\,(3.07)\,M_{\odot}$, and
$M_{\mathrm{Fe}}=1.39\,(1.50)\,M_{\odot}$ when using the combination V-NJ
(K-vL). For $25\,M_{\odot}$ models with $\eta=0.33$, we find
$M_{\mathrm{CO}}=6.38\,(6.40)\,M_{\odot}$, and
$M_{\mathrm{Fe}}=1.51\,(1.63)\,M_{\odot}$ for V-vL (K-dJ). Finally, the
$30\,M_{\odot}$ models with $\eta=0.33$ yield
$M_{\mathrm{CO}}=7.98\,(7.90)\,M_{\odot}$, and
$M_{\mathrm{Fe}}=1.56\,(1.58)\,M_{\odot}$ for the combination V-dJ (K-NJ). The
differences in $M_{\mathrm{CO}}$ (and to a lesser extent $M_{\mathrm{Fe}}$)
decrease with $M_{\mathrm{ZAMS}}$. However, note that the decreasing
difference is most likely caused by the lower wind efficiency $\eta=0.33$ for
the $25\,M_{\odot}$ and $30\,M_{\odot}$ models, while our $15M_{\odot}$ models
use full efficiency, i.e. $\eta=1.0$.

As anticipated in Sec. 3.5, the spread in $\xi_{2.5}$ increases between oxygen
depletion and the onset of core collapse: the final variations are about $\sim
30\%$ for models with the same initial mass (cf. Tab. 8). The two
$15\,M_{\odot}$ models have $\Delta\xi_{2.5}^{\mathrm{O\ depl}}=0.011$, and
$\Delta\xi_{2.5}^{\mathrm{pre-SN}}=0.029$. For the $25\,M_{\odot}$ models, the
spread at oxygen depletion is $\Delta\xi_{2.5}^{\mathrm{O\ depl}}=0.046$,
while it is $\Delta\xi_{2.5}^{\mathrm{pre-SN}}=0.081$ at the onset of core
collapse. The two $30\,M_{\odot}$ models go from $\Delta\xi_{2.5}^{\mathrm{O\
depl}}=0.015$ to $\Delta\xi_{2.5}^{\mathrm{pre-SN}}=0.082$.

Figure 8: Chemical composition profiles at the onset of core collapse. We show
all models computed to this point and listed in Tab. 8. Blue, cyan, yellow,
green, magenta, orange, and gray curves correspond to the mass fractions of
${}^{1}\mathrm{H}$, ${}^{4}\mathrm{He}$,${}^{12}\mathrm{C}$,
${}^{16}\mathrm{O}$, ${}^{20}\mathrm{Ne}$, ${}^{28}\mathrm{Si}$, and iron
group elements (i.e. with atomic number $50\leq A\leq 70$), respectively. Note
the different horizontal scale in each panel.

The pre-SN chemical abundances in the core are also affected by the choice of
the mass loss algorithm combination. Fig. 8 shows a comparison of the
distribution of the dominant isotopes for our pre-SN models. The composition
of the oxygen rich layer is sensitive to the early (hot phase) mass loss, with
the V scheme producing a lower ratio of
$X(^{16}\mathrm{O})/X(^{20}\mathrm{Ne})$, owing to the higher early mass loss
(cf. Sec. 3.2), and thus lower core temperature during the late phases (cf.
bottom panel of Fig. 6). The distribution of the chemical elements in mass
coordinate is also sensitive to the adopted mass loss algorithm combination.
However, it is likely to depend more strongly on the treatment of mixing
processes (which we do not vary here), mainly convection and overshooting,
which are the only processes fast enough to have an effect inside the CO core
of a star before the onset of collapse.

##  4 Discussion

##### Sensitivity of wind mass loss to evolving stellar properties.

For a given $M_{\mathrm{ZAMS}}$, the luminosity varies little between models
that experience different mass loss rates and the effective temperature (and
radius) varies much more. Once mass-loss induced differences in
$T_{\mathrm{eff}}$ appear, they feed back on mass loss to amplify these
differences. Hence, the $T_{\mathrm{eff}}$ evolution is most important for
governing the mass loss of models of a given $M_{\mathrm{ZAMS}}$. The
dependence of $T_{\mathrm{eff}}$ is particularly strong if the dust-driven vL
algorithm is included.

##### Role of dust in stellar winds.

The role of dust as a driver of RSG winds is a subject of debate (e.g., van
Loon et al., 2005; Ferrarotti & Gail, 2006; Bennett, 2010). Dust-driven mass
loss might occur in only a part of RSG evolution, and possibly even under only
rather specific conditions. Even if dust exists in the envelopes of cool,
evolved, massive stars, uncertainties in the grain properties result in highly
uncertain mass loss rates (van Loon et al., 2005). The extreme dust-driven vL
mass loss seen in our models is just one example (see Sec. 3.3). In this
context, it is important to mention that our models switch to the vL algorithm
already at $T_{\mathrm{eff}}\lesssim 15\,000\,\mathrm{K}$ ($\lesssim
11\,000\,\mathrm{K}$) when K (V) is used in the hot phase. These temperatures
are clearly far too hot for dust to form. However, the strong temperature
dependence $(\sim T^{-6.3})$ in combination with the extremely short time
spent in the Hertzsprung gap prevent this from having a substantial
consequence for the stellar mass. Hence, the extreme mass loss we find for
$M_{\mathrm{ZAMS}}\lesssim 20\,M_{\odot}$ with vL and $\eta=1$ is not an
artifact of how we use the vL algorithm. Within its framework the predicted
very low final masses for $15\,M_{\odot}$ and $20\,M_{\odot}$ stars,
$5.25-5.70\,M_{\odot}$ and $8.8\,M_{\odot}$, respectively, are robust. _It may
thus be possible to rule out the extreme case of vL with $\eta=1$ using pre-
explosion imaging and SN ejecta mass estimates._

We also note that in stars with $M_{\mathrm{ZAMS}}\gtrsim 20\,M_{\odot}$ the
vL algorithm self-regulates since its extreme mass loss uncovers hotter layers
of the stars. The steep temperature dependence of the vL algorithm then
results in lower mass loss rates than dJ and NJ.

##### The role of the efficiency parameter $\eta$.

The wind efficiency $\eta$ is the parameter that has the greatest impact on
the evolution of initially identical models (cf. Fig. 1 and Tab. 3). However,
it lacks an interpretation from first principles. Here, we investigate only
values $\eta\leq 1$, focusing on reduced wind mass loss motivated by possible
inhomogeneities in the wind structure (also called “clumpiness”, Smith, 2014;
Puls et al., 2008, and references therein). However, starting from first
principles, it cannot be excluded that clumps might actually enhance the mass
loss rate (Lucy & White, 1980). This might be the case if the overdense clumps
are efficiently pushed outward by impinging photons.

Furthermore, we assume a constant efficiency factor throughout the evolution,
but in principle the “clumpiness” of the wind might evolve (possibly even in a
stochastic way) and may require different $\eta$ in different evolutionary
phases. The mass loss routines in MESA333These are available at
https://stellarcollapse.org/renzo2017 are already adapted for varying this
parameter in different evolutionary phases. It is also possible that during
each single evolutionary phase we define, the efficiency of mass loss may vary
significantly, producing enhanced mass loss episodes separated by reduced mass
loss phases.

Although changing $\eta$ induces substantial changes in the mass loss rate and
final mass, the appearance (i.e., luminosity, effective temperature) is less
affected, making it difficult to use observed stellar populations to constrain
$\eta$ (Renzo, 2015).

##### Metallicity effects.

We do not investigate the effects of decreasing the metallicity in this study.
However, an approximate picture can be drawn by considering models with
reduced efficiency $\eta$ as proxy for low metallicity models. Most stellar
evolution codes implement the metallicity dependence of the wind mass loss by
just rescaling the mass loss rate at solar metallicity (cf. Eq. 1), which is
exactly the purpose of $\eta$. While the metallicity at the surface of the
star can change throughout the evolution, the main element driving a wind is
iron (Vink et al., 2001; Tramper et al., 2016), and its abundance is unlikely
to change because of upward mixing from the stellar interior. Therefore, it is
not unrealistic to consider a constant metallicity-related reduction factor
for the entire evolution. Nevertheless, the approach of considering reduced
$\eta$ as a proxy for lower metallicity does not take into account metallicity
effects on stellar radius and nuclear burning. These could indirectly affect
wind mass loss.

##### WR stars.

We emphasize in Sec. 2.2 the shortcomings of the computational definition of
WR stars adopted in stellar evolution codes ($X_{s}<0.4$). Although this
definition is artificial, we stress that the mass loss algorithms used during
this phase are derived from the observation of real WR stars. In the mass
range considered here, few WR stars are expected, and indeed our results show
that only $M_{\mathrm{ZAMS}}=35\,M_{\odot}$ stars with $\eta=1.0$ can become
sufficiently hydrogen-depleted at their surface to switch to a WR mass loss
algorithm. We compare only two WR mass loss algorithms (NL and H), and neglect
algorithms obtained by fitting either very luminous (i.e. more massive) or
hydrogen-free WR stars (e.g. Gräfener & Hamann, 2008; Tramper et al., 2016),
since none of our WR models reach the corresponding regions of the parameter
space. Therefore, in the framework of single nonrotating stars with wind
efficiency $\eta=1.0$, the minimum ZAMS mass to obtain a WR model is somewhere
between 30 and 35$\,M_{\odot}$. Although on the high end, this is in
relatively good agreement with the results obtained with other stellar
evolution codes (e.g. Woosley et al., 2002; Limongi & Chieffi, 2006; Eldridge
& Vink, 2006; Georgy et al., 2015). Even with full efficiency of the wind
before the WR phase, our models would likely underestimate the number of
observable WR stars.

Lowering $\eta$ has the obvious consequence of decreasing the mass loss rate,
and thus increasing the minimum ZAMS mass for single nonrotating WR stars.
However, the standard picture of a single nonrotating star misses pieces of
physics of great importance for the formation of WR stars (see, e.g., Maeder,
1996; Meynet & Maeder, 2003; Eldridge & Vink, 2006). These include, but are
not necessarily limited to, rotationally-enhanced mass loss, rotational mixing
processes (which can help depleting hydrogen from the surface, Meynet &
Maeder, 2003; Maeder, 1996), and binary interactions. Binarity can lead to the
formation of WR stars via envelope stripping in RLOF. Alternatively, accretion
or merger with a companion could increase the mass (and luminosity) of the
star sufficiently to enhance the wind mass loss and remove the hydrogen-rich
envelope. Also note that the envelope hydrogen depletion needed to switch to a
WR mass loss algorithm can be reached also because of upward mixing of
thermonuclearly processed material, e.g., because of efficient rotational
mixing, and not only because of mass loss. The choice for the algorithmic
representation of mixing processes (see Appendix B) influences directly the
surface mass fraction of hydrogen, but also the core mass, and consequently
the luminosity. Indirectly, these effects can change the mass loss rate and
consequently the fate of a model from/to WR.

##### Nucleosynthetic yields.

The nucleosynthetic yields of massive stars are mass loss (and angular
momentum loss) dependent (see, e.g., Maeder, 1992; Frischknecht et al., 2016).
Processes such as rotational mixing can bring thermonuclearly processed
material upward that is then lost through winds. This is especially relevant
for s-process elements which are synthesized during the hydrostatic lifetime
of massive stars and the ratio of carbon to oxygen abundance in the stellar
wind yields. On top of this, the success or failure of the SN explosion, and
the details of the explosive nucleosynthesis, depend on the interior structure
of the exploding star and thus on its mass loss history (e.g., Sukhbold et
al., 2016).

##### Consequences for SN explosions and compact remnants.

We find that mass loss affects the core structure and the burning shells
surrounding the core. Mass loss during the hot phase of the evolution (i.e.
the main sequence, roughly speaking) is important for the core structure,
because the core itself re-adjusts quasi-statically to the wind from the
stellar atmosphere. During this phase, different algorithms produce small
variations in the core, which are then amplified by the subsequent evolution
(cf. Sec. 3.2 and Fig. 4). The general trend is that a higher mass loss rate
during the hot phase produces structures with lower core compactness (cf. Sec.
3.6). The cool phase mass loss also impacts the core compactness, but more
indirectly, through its effect on the burning shells. Most of the mass is lost
during the cool phase of the evolution (cf. Sec. 3.3), and the vigor, extent
and type of mixing in the burning shells all depend on the amount and the
timing of mass loss. Since our present understanding of core-collapse SN
explosions strongly depends on the details of the input stellar models (e.g.,
Janka et al., 2012; Couch & Ott, 2015; Chatzopoulos et al., 2016), overlooking
the impact of wind mass loss on the core structure might bias detailed
hydrodynamical simulations of stellar explosions. This, in turn, can have
significant implications for the NS/BH ratio and mass distribution, and
consequently also for the inferred gravitational wave sources. We
provide444Data are available at https://zenodo.org/record/292924#.WK_eENWi60i
and input parameter files at https://stellarcollapse.org/renzo2017 stellar
models at the end of the main sequence, at the end of the hot phase of the
evolution, at the end of the mass loss phase ($T_{c}\geq
10^{9}\,\mathrm{[K]}$), at oxygen depletion, and at the onset of core
collapse. These models can be used as starting points for stellar experiments
during late evolutionary phases (see e.g. Couch & Ott, 2015; Chatzopoulos et
al., 2016).

##### Observational $M_{\mathrm{ZAMS}}$ estimates of SN progenitors.

The large vertical spread in final mass caused by differences in mass loss for
a given $M_{\mathrm{ZAMS}}$ (Fig. 1) suggests at first sight that these large
uncertainties may map to equally large uncertainties in $M_{\mathrm{ZAMS}}$
estimates from pre-explosion observations. However, the pre-SN mass is not a
direct observable. Rather, $M_{\mathrm{ZAMS}}$ estimates typically rely on
observational measurements of luminosity and effective temperature that are
then compared with stellar models (e.g., Smartt 2009). In Fig. 9, we plot the
final luminosity $L$ and effective temperature $T_{\mathrm{eff}}$ for our
entire model set. Table 4 summarizes the numerical results. These results
demonstrate that wind mass loss variations have very little effect on the
final luminosity for stars with $M_{\mathrm{ZAMS}}\lesssim 30\,M_{\odot}$ and
luminosity variations are smaller than typical observational uncertainties
(see, e.g., Smartt 2009). The similarity in luminosity of models of a given
$M_{\mathrm{ZAMS}}$ is a consequence of the rather small effect that mass loss
has on the core mass (cf. Tab. 3). Interestingly, $T_{\mathrm{eff}}$
variations are also small for $M_{\mathrm{ZAMS}}\lesssim 30\,M_{\odot}$. The
maximum variation from the average $T_{\mathrm{eff}}$ of a ZAMS mass is only
$0.08\,\mathrm{dex}$ and comes from the $\eta=1$ K-vL and V-vL $20\,M_{\odot}$
models that are YSGs at the end of their lives. Wind mass-loss dependent
variations in $L$ and $T_{\mathrm{eff}}$ are much larger for $35\,M_{\odot}$
models, some of which die as BSGs and some as WR stars.

Given the above results, it appears that for massive stars with
$M_{\mathrm{ZAMS}}\lesssim 30\,M_{\odot}$ wind mass loss uncertainties do not
increase the overall level of uncertainty with which the SN progenitor
$M_{\mathrm{ZAMS}}$ can be estimated from pre-explosion observations. However,
our results do not say anything about the other possibly existing mass loss
channels (e.g., binary interactions and/or impulsive phenomena, see, for
example, Smith, 2014; Smith & Arnett, 2014; Morozova et al., 2015; Margutti et
al., 2017) that are usually neglected in stellar evolution calculations, or
included in an very simplified way using enhanced winds (see, e.g., Meynet et
al., 2015).

Figure 9: Spread in pre-SN appearance due to the uncertainty in wind mass
loss. Each color represents an initial mass $M_{\mathrm{ZAMS}}$, each marker
corresponds to a different efficiency factor $\eta$. The vertical dot-dashed
lines indicate the BSG-YSG boundary, and YSG-RSG boundary according to Georgy
et al. (2016). $35\,M_{\odot}$ models exhibit a larger spread, because some
models (indicated by a thick black edge on their markers) develop into WR
stars.

##### The RSG problem.

The observed lack of RSG SN progenitors with $M_{\mathrm{ZAMS}}\gtrsim
16\,M_{\odot}$ (Smartt et al., 2009; Smartt, 2009) might be explained with the
effects of mass loss on the pre-SN stellar appearance (Georgy et al., 2016).
Table 4 and Fig. 9 show that different wind algorithms can change the pre-
explosion appearance of a massive star. However, this effect is relatively
small in luminosity and $T_{\mathrm{eff}}$, and the relative number of YSG to
RSG we find (cf. Sec. 3) largely depends on the adopted definition of YSG vs.
RSG. Therefore, the systematic uncertainty in the treatment of wind mass loss
for $16\lesssim M_{\mathrm{ZAMS}}/M_{\odot}\lesssim 30$ does not seem
sufficient to solve the RSG problem. It is likely that either (_i_) stars with
initial mass in the range $\sim 16-30\,M_{\odot}$ do not produce a bright
transient at their death, or (_ii_) mass loss phenomena other than winds,
happening at late stages in the evolution, change the pre-SN appearance of the
star. SN observations can constrain the total ejected mass and pre-explosion
imaging might also provide constraints (although subject to large
uncertainties) on the total mass lost and on the mass loss timing (if the
surface chemical composition can be inferred, e.g., Smartt et al., 2009;
Gordon et al., 2016).

##### Potential observational constraints.

A range of special systems, events, and phenomena offer alternative and
complimentary ways to constrain mass loss to the more traditional spectral
observations of massive stars. These special systems include bow shocks of
runaway stars (e.g., Gull & Sofia, 1979; Meyer et al., 2016), SN shocks
running into the circumstellar material (e.g., Maeda et al., 2015; Chakraborti
et al., 2016; Margutti et al., 2017), flash-spectroscopy of material ejected
shortly before core collapse (e.g., Khazov et al., 2016), accretion in wind-
fed high mass X-ray binaries, or binary wind collisions. Constraints from
special systems can be combined with those arising from observed populations
of stars and their compact remnants (including gravitational wave sources).
The different mass loss algorithms that we compare here give a range of mass
loss timing, which suggests that the chemical composition and dust properties
of the circumstellar material may also contain hints regarding the mass loss
of massive stars.

##### Other known unknowns and uncertainties.

In this study, we have only varied the mass loss algorithms used throughout
the evolution and their efficiency parameter $\eta$. However, other
uncertainties (with possible degeneracies) are known, for example in the
treatment of rotation, magnetic fields (e.g., Petit et al., 2017), convective
mixing, and overshooting (e.g., Arnett et al., 2015; Arnett, 2015; Arnett &
Meakin, 2016; Farmer et al., 2016). The coupling of these uncertainties with
the wind mass loss may modify the outcomes of our numerical experiments.

##  5 Conclusions

Massive star mass loss is a longstanding issue in stellar evolution. Despite
decades of observational and theoretical work it remains incompletely
understood. Mass loss influences the lifetime and appearance of massive stars,
their internal structure at the onset of core collapse, and their total
nucleosynthetic yields. Through its effect on the pre-supernova (pre-SN)
structure, wind mass loss can impact the outcome of core collapse and the
nature of the compact remant, with potential implications for gravitational
wave astronomy.

We studied the impact of a broad range of wind mass loss algorithm
combinations on the evolution and pre-collapse structure of nonrotating,
single, solar-metallicity stars with initial masses of $15,\,20,\,25,\,30,$
and $35\,M_{\odot}$. We compared 12 different mass loss algorithm
combinations, drawing from 2 algorithms for the hot phase of the evolution
(corresponding roughly to the main sequence), 3 for the cool phase of
evolution, and 2 for the Wolf-Rayet phase (if it is reached). We explored the
effects of reducing the mass loss rate with an efficiency scaling factor
$\eta=1.0,\,0.33$, or $0.1$ to crudely account for reduced stellar mass loss
that could be caused, for example, by inhomogeneities in the wind (e.g.,
clumpiness). The resulting differences in stellar structure and total mass at
various stages of the evolution are caused by the different algorithmic
representations of stellar winds.

The different mass loss efficiencies and algorithm combinations have profound
effects on the evolution and the pre-SN masses of massive stars. On the one
hand, this can be expected given the inherent differences of the various mass
loss algorithms and the various assumptions that enter them. On the other
hand, these algorithms all attempt to describe the same physical process –
steady wind mass loss – and less sensitivity to the theoretical/empirical
treatment of mass loss would in general be desirable.

We find that the choice of wind efficiency scaling factor $\eta$ has the
greatest impact on our stellar models. It affects their total mass loss, their
evolutionary path, and their pre-SN structure. $\eta$ is therefore the main
uncertainty and limiting factor for our present understanding of wind mass
loss and its effects. If the wind efficiency is low, the differences between
various mass loss algorithms are less important.

Considering the full range of wind efficiencies $0.1\leq\eta\leq 1.0$, we find
that there is a $\sim 50\%$ uncertainty in the pre-SN mass for a given
$M_{\mathrm{ZAMS}}$. For fixed efficiency $\eta=1.0$, the uncertainty varies
with initial mass and is $\sim 15-30\%$ in most cases. Impulsive mass loss
events (eruptions, pulsational instabilities, etc., all neglected here) could
only make the uncertainties in the initial to pre-SN mapping more severe.

Despite the large uncertainty in the pre-SN mass, we find that the key
observables from pre-SN imaging, the luminosity and effective temperature
before explosion, are only mildly affected by varying the mass loss algorithm
combination and wind efficiency. The uncertainties in $L$ and
$T_{\mathrm{eff}}$ from our models are within observational limits, suggesting
that wind mass-loss uncertainties do not affect observational estimates of SN
progenitor masses from pre-SN observations for most massive stars (assuming a
single star evolution scenario). Nevertheless, the impact of stellar winds on
the internal structure of the star can affect the mass and composition of the
SN ejecta.

Independent of the employed algorithm during the hot phase, the amount of mass
lost in this phase is only a small fraction ($\sim$ few percent) of the total
mass. However, since the core can respond directly to it, mass loss during the
hot phase creates seed differences that grow during the subsequent evolution,
leading to changes in the pre-SN structure and composition of the stellar
core. Later, burning shells re-adjust to the mass loss instead of the core
itself, and the effect of mass loss is more indirect.

Most of the mass is lost during the late and short cool phase of the
evolution. Wind mass loss during this phase is more uncertain. On the cool
side of the Hertzsprung-Russel diagram, two different wind driving mechanism
might exist: line-driving or dust-driving. The latter, assumed by the van Loon
et al. (2005) (vL) algorithm, has a much stronger temperature dependence and
produces much higher mass loss rates than algorithms describing line-driven
mass loss. vL mass loss can be so strong that it reveals the deep and hot
layers of the star. This results in self-damping of the wind itself for higher
mass progenitors and thus a higher value of the final mass. The vL algorithm
produces very different evolutionary tracks from those obtained with
algorithms assuming line driving, with the vL algorithm driving blueward
displacements on the Hertzsprung-Russel diagram. The more indirect effect of
late mass loss on the pre-SN core structure is difficult to pinpoint with our
limited model grid.

In our model grid, only models with $M_{\mathrm{ZAMS}}\geq 35\,M_{\odot}$ and
full wind efficiency ($\eta=1.0$) develop into Wolf-Rayet stars. The absence
of Wolf-Rayet models with reduced wind efficiency suggests that either (_i_)
our initial mass range is too small to produce Wolf-Rayet stars from single,
nonrotating stars, or (_ii_) other formation channels, such as binarity or
impulsive mass loss events, might be dominant. The standard picture for the
evolution of massive stars, the so-called “Conti scenario” (e.g., Maeder &
Conti, 1994; Lamers, 2013; Smith & Tombleson, 2015), predicts the formation of
Wolf-Rayet stars with $M_{\mathrm{ZAMS}}\lesssim 40\,M_{\odot}$, therefore, if
mass loss occurs in nature with reduced efficiency ($\eta<1$), then our
results with reduced wind efficiency disagree with this prediction.

From our limited sample of six models that we were able to evolve to the pre-
SN stage, we find that changing the wind mass loss algorithm combination can
lead to changes of the compactness parameter $\xi_{2.5}$ of up to $30\%$ for a
given $M_{\mathrm{ZAMS}}$. Moreover, the pre-SN models show a spread in terms
of composition profiles, density profiles, and core masses. These
uncertainties add to those arising from the incomplete understanding of mixing
processes (mainly convection and overshooting). They complicate the study of
the core-collapse SN explosion mechanism by adding uncertainty to the initial
conditions from which core-collapse SN simulations start. _This finding
underlines that systematic uncertainty in massive star mass loss can have
important implications for the relative number of BHs and NSs resulting from
core collapse events, and consequently for gravitational wave sources, and for
the nucleosynthetic yields from the explosions of massive stars._

Although the present study provides new insights into the effects of wind mass
loss on the evolution and pre-SN structure and appearance of massive stars, it
suffers from a number of important limitations that must be addressed by
future work.

Our model grid coarsely samples a limited mass range, includes only
nonrotating solar-metalicity models, and we could evolve only six models to
the onset of core collapse. Extending this grid to finer mass sampling, higher
masses, lower metallicity, and evolving all models to the pre-SN stage will be
important future work needed to infer more robustly how wind mass loss affects
pre-SN structure. Furthermore, we did not consider time-dependent wind