11institutetext: Zentrum für Astronomie der Universität Heidelberg,
Astronomisches Rechen-Institut, Mönchhofstr. 12-14, 69120 Heidelberg,
Germany
11email<EMAIL_ADDRESS>22institutetext: Department of
Astronomy, University of Geneva, Chemin Pegasi 51, 1290 Versoix, Switzerland

# Impact of main-sequence mass loss on the appearance, structure and evolution
of Wolf-Rayet stars

J. Josiek Impact of main-sequence mass loss on the appearance, structure and
evolution of Wolf-Rayet starsImpact of main-sequence mass loss on the
appearance, structure and evolution of Wolf-Rayet starsImpact of main-sequence
mass loss on the appearance, structure and evolution of Wolf-Rayet starsImpact
of main-sequence mass loss on the appearance, structure and evolution of Wolf-
Rayet stars S. Ekström Impact of main-sequence mass loss on the appearance,
structure and evolution of Wolf-Rayet starsImpact of main-sequence mass loss
on the appearance, structure and evolution of Wolf-Rayet stars A.A.C. Sander
Impact of main-sequence mass loss on the appearance, structure and evolution
of Wolf-Rayet starsImpact of main-sequence mass loss on the appearance,
structure and evolution of Wolf-Rayet stars

(Received 19 January 2024 / Accepted 22 April 2024)

###### Abstract

Context. Stellar winds are one of the most important drivers of massive star
evolution and a vital source of chemical, mechanical, and radiative feedback
on the galactic scale. Despite its significance, mass loss remains a major
uncertainty in stellar evolution models. In particular, the interdependencies
of different approaches with subsequent evolutionary stages and predicted
observable phenomena are far from being systematically understood.

Aims. In this study, we examine the impact of main sequence mass loss on the
structure of massive stars throughout their entire evolution. A particular
focus is placed on the consequences for entering the Wolf-Rayet (WR) regime
and the subsequent evolution.

Methods. Using the Geneva stellar evolution code, we compute grids of single,
non-rotating stellar models at solar and Large Magellanic Cloud (LMC)
metallicity of initial masses between 20 and 120 solar masses, with two
representative prescriptions for high and low main sequence mass loss.

Results. We obtain detailed numerical predictions regarding the structure and
evolution of massive stars, and infer the role of main sequence mass loss by
comparison of the mass-loss rate prescriptions. We present implications for
the overall evolutionary trajectory, including the evolution of WR stars, as
well as the effect on stellar yields and stellar populations.

Conclusions. Mass loss during the main sequence plays an important role due to
its ability to affect the sequence and duration of all subsequent phases. We
identify several distinct evolutionary paths for massive stars which are
signiﬁcantly influenced by the chosen main sequence mass-loss description. We
also discuss the impact of uncertainties other than mass loss on the
evolution, in particular those relating to convection. We further demonstrate
that not just the total mass loss, but the specific mass-loss history
throughout a star’s life is a crucial determinant of many aspects, such as the
resulting stellar yields.

###### Key Words.:

Stars: evolution – Stars: massive – Stars: mass-loss – Stars: Wolf-Rayet

##  1 Introduction

Through their powerful winds and high ionizing flux, massive stars play a
vital role in shaping the physical and chemical structure of their
surroundings (Krause et al. 2013). Unfortunately, massive stars are greatly
outnumbered by their low-mass counterparts, making them much more challenging
to study. Some phenomena related to massive stars are so rare that we only
have a handful of observations of them. For instance, the number of Wolf-Rayet
(WR) stars discovered to date is counted in the hundreds (Massey et al. 2014;
Mauerhan et al. 2011; Rosslowe & Crowther 2015; Neugent & Massey 2019), an
incredibly small number compared to the billions of stars we know in total.

One of the most important ingredients in stellar models is their mass-loss
rate. Within the course of their lives, single massive stars can eject more
than half of their own initial mass in the form of stellar winds, which not
only transfer energy and momentum into the interstellar medium (ISM), but also
feed it with new gas that has been chemically enriched with products of
nucleosynthesis. Through mass loss, stars also change their own appearance as
they progressively remove layers from their surface. When modeling the
evolution of stars, this mass loss cannot be predicted within the usual
physical framework of the structure models, and thus needs to be taken into
account in a prescribed way. Numerous different methods have been proposed to
measure, calculate and predict the mass-loss rate of stars, each with their
own caveats and uncertainties. For this reason, it is important to consider a
wide range of possible mass loss scenarios when attempting to explain stellar
evolution as completely as possible.

While on the main sequence (MS), massive stars appear as hot OB-type stars
with surface temperatures between 10 000–50 000 K. During this phase, stars
lose mass mainly as a result of line-driven winds, i.e. through the absorption
of photons in specific transitions of atoms in the atmosphere (see, e.g., Vink
2022, for a comprehensive review). In this work, we consider two recipes which
can be used to compute the mass-loss rates in this phase and study their
impact on the later phases of evolution.

The MS is arguably the most inconspicuous stage of evolution, marked by slow
processing of hydrogen in the core, and comparatively weak mass loss from the
surface. It is therefore often pushed into the background in favor of studying
phases of stronger and less certain mass-loss rates such as the red supergiant
(RSG) phase (e.g., Massey et al. 2023; Beasor & Davies 2018) or the WR phase
(e.g., Sander et al. 2020, 2022). However, there is up to one order of
magnitude of discrepancy between different mass-loss descriptions for OB stars
(e.g., Vink et al. 2001; Brands et al. 2022; Björklund et al. 2023), which is
large enough to introduce major uncertainties into their evolution, even
without considering the uncertainties of post-MS mass loss. The goal of this
paper is to perform a systematic study on the effects of MS mass loss on the
evolution of massive stars. For this, we compute grids of evolutionary models
using two different prescriptions for the MS mass-loss rate and compare the
outcome.

In Sect. 2, we review the two MS mass-loss prescriptions compared in this
study and briefly summarize the mass-loss prescriptions used in the post-MS
phases of evolution. We then present the parameters of the model grid in Sect.
3. In Sect. 4 we present a preliminary characterization of the mass-loss
rates, effects on the evolution during the MS, as well as consequences for
stellar yields. Sect. 5 is dedicated to the formation, structure, and
evolution of WR stars as a key phase of evolution for many massive stars. In
Sect. 6, we will present the impact of MS mass loss on the evolutionary
endpoints of the models as possible SN progenitors as well as populations of
massive stars as a whole. Finally, we will discuss the findings in Sect. 7 and
summarize the results in Sect. 8.

##  2 Mass loss descriptions

###  2.1 The main sequence

In this work, we employ two mass-loss prescriptions that will act as
prototypes for stronger and weaker mass loss on the MS. The de-facto standard
MS mass-loss prescription that is used in most of the current grids of
evolutionary models and population synthesis codes is the one introduced by
Vink et al. (2001, hereafter referred to as Vin01). This recipe is distilled
from a set of numerical stellar wind models and provides the mass-loss rate as
a function of the global stellar parameters, split over multiple effective
temperature domains. When the temperature transitions into a different domain,
the wind models predict the recombination or ionization of iron, causing a
significant change in the available opacities to drive the wind, which in the
models leads to a discontinuity in the trends for the predicted mass-loss
rates. Motivated by the bi-stable wind solution found by Pauldrach & Puls
(1990) for P Cygni and the jump in the observed ratio between terminal and
escape velocity by Lamers et al. (1995), these discontinuities are referred to
as bistability jumps. The existence of the jumps in $\dot{M}$ and their
magnitude are an active topic of theoretical and empirical research (e.g.,
Petrov et al. 2016; Björklund et al. 2023; Bernini-Peron et al. 2023; Krtička
et al. 2021, 2024).

In our models implementing the Vin01 mass-loss prescription, we account for
both the Fe IV/III as well as the Fe III/II jumps. First, the so-called
“characteristic density” $\langle\rho\rangle$, defined as the density of the
wind at 50% terminal velocity, is calculated via Eq. (23) in Vin01, i.e.,

$\log\langle\rho\rangle=-14.94+0.85\log(Z/Z_{\odot})+3.2\,\Gamma_{\mathrm{e}}$
(1)

with $\Gamma_{\mathrm{e}}$ denoting the electron-scattering Eddington
parameter. This parameter encodes the ratio between $L$ and $M$ and thus the
ratio between the outward force caused by radiation pressure and the inward
pull by gravity, which is inherently linked to radiation-driven mass loss.
From $\langle\rho\rangle$ we then obtain the location of the bistability jumps
with

$T_{\mathrm{eff,jump1}}/\mathrm{K}=61\,200+2\,590\,\log\langle\rho\rangle\text{,}$
(2)

and

$T_{\mathrm{eff,jump2}}/\mathrm{K}=100\,000+6\,000\,\log\langle\rho\rangle\text{.}$
(3)

For the hot domain, i.e. $T_{\mathrm{eff}}>T_{\mathrm{eff,jump1}}$, the mass-
loss rate is then determined as

$\begin{split}\log\dot{M}&=-6.697\\\
&+2.194\log\left(\frac{L}{10^{5}L_{\odot}}\right)\\\
&-1.313\log\left(\frac{M}{30\,M_{\odot}}\right)\\\
&-1.226\log\left(\frac{v_{\infty}/v_{\mathrm{esc}}}{2.0}\right)\\\
&+0.933\log\left(\frac{T_{\mathrm{eff}}}{40\,000\,\mathrm{K}}\right)\\\
&-10.92\left[\log\left(\frac{T_{\mathrm{eff}}}{40\,000\,\mathrm{K}}\right)\right]^{2}\\\
&+0.85\log\left(\frac{Z}{Z_{\odot}}\right)\text{,}\end{split}$ (4)

with $v_{\infty}/v_{\mathrm{esc}}=2.6$. For the intermediate domain, i.e.
$T_{\mathrm{eff,jump1}}>T_{\mathrm{eff}}>T_{\mathrm{eff,jump2}}$, we use

$\begin{split}\log\dot{M}&=-6.688\\\
&+2.21\log\left(\frac{L}{10^{5}L_{\odot}}\right)\\\
&-1.339\log\left(\frac{M}{30\,M_{\odot}}\right)\\\
&-1.601\log\left(\frac{v_{\infty}/v_{\mathrm{esc}}}{2.0}\right)\\\
&+1.07\log\left(\frac{T_{\mathrm{eff}}}{20\,000\,\mathrm{K}}\right)\\\
&+0.85\log\left(\frac{Z}{Z_{\odot}}\right)\text{,}\end{split}$ (5)

with $v_{\infty}/v_{\mathrm{esc}}=1.3$. For the cool temperature domain, i.e.
$T_{\mathrm{eff}}<T_{\mathrm{eff,jump2}}$, we also use Eq. (5), but increase
the leading constant to $-5.99$ and set $v_{\infty}/v_{\mathrm{esc}}=0.7$. The
unit of $\dot{M}$ is as usual implied to be $M_{\odot}\,\text{yr}^{-1}$, which
will be valid throughout this paper without explicitly stating so.

Alternatively to directly modeling the winds, there are also empirical and
semi-empirical descriptions of mass loss. Based on the wind theory by Castor
et al. (1975) combined with the Eddington stellar model, Bestenlehner (2020,
hereafter known as Bes20) suggested a description of the mass-loss rate
depending mainly on $\Gamma_{\mathrm{e}}$. Bes20 proposed a mass-loss recipe
in the form

$\log\dot{M}=\log\dot{M}_{0}+\left(\frac{1}{\alpha}+0.5\right)\log(\Gamma_{\text{e}})-\left(\frac{1}{\alpha}+1\right)\log(1-\Gamma_{\text{e}})\text{,}$
(6)

where $\alpha$ is a force-multiplier parameter used to account for the effect
of line opacities, and $\dot{M}_{0}$ is a global scale factor. This recipe
intends to connect the regime of stronger, optically thick winds of very
massive WNh stars (e.g., Vink et al. 2011; Bestenlehner et al. 2020; Sabhahit
et al. 2022) with the weaker optically thin winds of OB stars. Specifically,
for thin winds (where $\Gamma_{\mathrm{e}}\ll 1$), the
$\log(\Gamma_{\mathrm{e}})$-term dominates, whereas for thick winds (where
$\Gamma_{\mathrm{e}}\rightarrow 1$) the $\log(1-\Gamma_{\mathrm{e}})$-term is
dominant. At the transition point between the two regimes, where both
$\Gamma_{\mathrm{e}}$-terms cancel each other, the mass-loss rate is
$\dot{M}_{0}$. As a semi-empirical recipe, both of its parameters
$\dot{M}_{0}$ and $\alpha$ must be calibrated on external data, as was done
for example by Bestenlehner et al. (2020) and Brands et al. (2022) on the R136
cluster in the Large Magellanic Cloud (LMC). In this study, we do not aim for
a detailed analysis of individual mass-loss recipes, but employ Bes20 as a
prototype of considerably weaker MS mass loss compared to Vin01. Thus, we use
the values of $\log\dot{M}_{0}=-5.19$ and $\alpha=0.456$ from Brands et al.
(2022) for all models computed using the Bes20 recipe, including those at
Galactic metallicity, noting that the mass loss in these models will be
underestimated when using LMC-calibrated data.

###  2.2 Validity domains of the mass-loss prescriptions

Due to the changing physical conditions in a star’s atmosphere, the mechanisms
that drive its mass loss vary considerably throughout the course of its
evolution. Therefore, the implementation of a single mass loss description
must be restricted to the domain where it is applicable based on the
underlying physical assumptions, and multiple different mass loss descriptions
are used to cover the entire evolution of a star.

For the Vin01 models, we apply this mass-loss rate in the same domain defined
in previous GENEC models (Ekström et al. 2012), i.e., when
$\log(T_{\text{eff}}[\text{K}])>3.9$ and the surface hydrogen mass fraction is
above $30\%$. For models where we apply the Bes20 prescription, we use their
recommended domain of $T_{\text{eff}}>30\,000\,\text{K}$ and the surface
hydrogen mass fraction above $10^{-5}$, since this satisfies the assumptions
of the Eddington stellar model and includes WNh stars, which are also
adequately described by this prescription according to Bestenlehner (2020).
Although the domains of these two mass-loss prescriptions are not equivalent,
they both cover most of the MS, which is the longest phase and also the one
during which we are interested in comparing different mass-loss rates.

###  2.3 Post-main-sequence mass-loss rates

For parts of the evolution not covered by the aforementioned prescriptions, we
use the same mass-loss scheme as described in the GENEC grids by Ekström et
al. (2012) , which we briefly summarize here. RSG mass loss is calculated as
follows, based on a fit of Fig. 3 in Crowther (2001):

$\log\dot{M}=-13.83+1.7\log(L/L_{\odot})\text{.}$ (7)

The mass-loss rate of WR stars is computed with Nugis & Lamers (2000), with an
additional $Z$-scaling from Eldridge & Vink (2006) based on the calculations
of Vink & de Koter (2005). For WN stars, the following equation is used:

$\begin{split}\log\dot{M}&=-13.6+1.63\log(L/L_{\odot})+2.22\log(X_{\mathrm{He}})\\\
&+0.85\log(Z_{\mathrm{ini}}/Z_{\odot})\text{.}\end{split}$ (8)

For WC and WO stars, the mass-loss rate is as follows:

$\begin{split}\log\dot{M}&=-8.3+0.84\log(L/L_{\odot})+2.04\log(X_{\mathrm{He}})+1.04\log(Z)\\\
&+0.66\log(Z_{\mathrm{ini}}/Z_{\odot})\text{.}\end{split}$ (9)

In the small domain covered by Gräfener & Hamann (2008), this prescription is
used instead for WR stars:

$\begin{split}\log\dot{M}&=10.046+\beta\log(\Gamma_{\mathrm{e}}-\Gamma_{0})\\\
&-3.5\log(T_{\mathrm{eff}}[\mathrm{K}])+0.42\log(L/L_{\odot})-0.45X_{\mathrm{H}}\text{,}\end{split}$
(10)

with $\Gamma_{0}=0.326-0.301\log(Z/Z_{\odot})-0.045\log(Z/Z_{\odot})^{2}$,
and$\beta=1.727+0.25\log(Z/Z_{\odot})$.

BSG mass loss is computed in the same way as MS OB stars, since their stellar
parameters fall into the accepted validity domain of those prescriptions.

Outside the validity domain of all specialized mass-loss prescriptions, we
apply the prescription by de Jager et al. (1988):

$\begin{split}\log\dot{M}=\sum_{j=0}^{4}\sum_{i=0}^{5-j}a_{ij}T_{i}&\left(\frac{\log(T_{\mathrm{eff}}[\mathrm{K}])-4.05}{0.75}\right)\\\
&\times
T_{j}\left(\frac{\log(L/L_{\odot})-4.6}{2.1}\right)\text{.}\end{split}$ (11)

where $T_{j}(x)=\cos(j\arccos x)$ are Chebychev polynomials, and $a_{ij}$ are
fitted coefficients found in Table VI of de Jager et al. (1988). The
definition of each prescription’s validity domain results in the above formula
effectively only being used for YSGs.

##  3 Evolution modeling

In order to investigate the effect of mass loss on the evolution of massive
stars, we computed grids of evolutionary models using the Geneva stellar
evolution code (GENEC, Eggenberger et al. 2008). Specifically, we computed
models with the following parameter grid:

  * •

Initial mass: 20, 25, 30, 40, 50, 60, 66, 73, 80, 85, 95, 105, 120
$M_{\odot}$,

  * •

Metallicity: 0.014 (Solar), 0.006 (LMC),

  * •

OB mass-loss prescription: Vin01, Bes20.

The full grid contains one model for each possible combination of these three
parameters, resulting in a total of 52 models computed for this work. The set
of initial masses was chosen to be roughly equidistant in the HRD at the zero-
age MS (ZAMS), but having a slightly finer resolution in the higher mass
range. The exact initial elemental abundances for the models at solar and LMC
metallicity are described in Ekström et al. (2012) and Eggenberger et al.
(2021), respectively. All of the models are non-rotating in order to avoid
introducing an additional source of uncertainty as well as new effects that
may blur any conclusions pertaining to mass loss only. Convection was treated
using the Schwarzschild criterion and a step-overshoot scheme with
overshooting parameter $\alpha_{\mathrm{ov}}=0.1$. For all other physical
ingredients not explicitly mentioned here, we use the same setup as described
in Ekström et al. (2012). The models were computed until the core carbon mass
fraction drops below $10^{-5}$ at the end of core carbon burning, except both
solar $40\,M_{\odot}$ models, as well as the solar $25\,M_{\odot}$ model using
Vin01, for which convergence issues prevented the computation from finishing.
These models will be excluded from analysis at the relevant points.

##  4 Early evolution and yields

###  4.1 Characterization of mass-loss prescriptions

From the evolution models we obtain that, at solar metallicity, Vin01 produces
a mass-loss rate that is systematically around one order of magnitude above
Bes20 across the entire initial mass range. At LMC metallicity, the difference
between the two mass-loss rates is smaller, as the Vin01 recipe scales down
for lower metallicity, whereas our implementation of Bes20 does not.
Independent of their accuracy, the large contrast in the resulting mass-loss
rate is beneficial for the purpose of comparing the effect of weaker versus
stronger MS mass loss on stellar evolution.

The Vin01 mass-loss rate is much more variable during the MS for stars above
60 solar masses. This can be attributed to the fact that these more massive
stars expand more during their MS, and thus the cooler atmospheres cause the
mass loss to be boosted due to the bistability jumps. To illustrate the
evolution of the mass-loss rate during the MS using the different
prescriptions, we show an example in Fig. 1. As expected, Bes20 produces a
monotonic increase in the mass-loss rate of around one order of magnitude
during the MS as the star steadily increases its luminosity and therefore its
Eddington parameter. Close to the end of the MS, the star is too cool and thus
we switch to the stronger mass-loss rate by de Jager et al. (1988).

Figure 1: Evolution of the mass-loss rate during the MS of the $60\,M_{\odot}$
models at solar metallicity using the Vin01 and the Bes20 prescriptions.

###  4.2 Main-sequence evolution

Figure 2: Total mass of the models remaining at the terminal-age main sequence
(TAMS) as a function of initial mass at solar metallicity (left) and LMC
metallicity (right). The data is fitted with a piece-wise linear function to
indicate the two distinct mass-loss regimes. The gray dashed line indicates
the 1:1 relation in both plots.

Fig. 2 shows the mass remaining at the terminal-age main sequence (TAMS),
defined here as the moment in time where the central hydrogen mass fraction of
the star drops below $10^{-5}$. The results show a clear division into two
distinct regimes of MS mass loss for both of the considered OB mass-loss
prescriptions and at both metallicities. In each regime, the TAMS mass has a
linear dependence on the initial mass, parametrized as
$M_{\mathrm{TAMS}}=aM_{\mathrm{ini}}+b$, for which we determine the best fit
parameters $a$ and $b$ using the Python routine optimize.curvefit provided by
the scipy package (Virtanen et al. 2020). The results are shown in Table 1,
and drawn in Fig. 2.

$Z$ | OB $\dot{M}$ prescription | $a$ | $b/M_{\odot}$
---|---|---|---
$0.014$ | Bes20 | ($M_{\mathrm{ini}}\leq 66\,M_{\odot}$) | $0.869\pm 0.014$ | $3.1\pm 0.6$
|  | ($M_{\mathrm{ini}}\geq 80\,M_{\odot}$) | $0.601\pm 0.007$ | $4.2\pm 0.7$
| Vin01 | ($M_{\mathrm{ini}}\leq 50\,M_{\odot}$) | $0.76\pm 0.04$ | $5.3\pm 1.3$
|  | ($M_{\mathrm{ini}}\geq 60\,M_{\odot}$) | $0.458\pm 0.019$ | $10.0\pm 1.7$
$0.006$ | Bes20 | ($M_{\mathrm{ini}}\leq 85\,M_{\odot}$) | $0.866\pm 0.012$ | $3.6\pm 0.7$
|  | ($M_{\mathrm{ini}}\geq 93\,M_{\odot}$) | $0.56\pm 0.08$ | $19\pm 8$
| Vin01 | ($M_{\mathrm{ini}}\leq 73\,M_{\odot}$) | $0.844\pm 0.014$ | $3.6\pm 0.7$
|  | ($M_{\mathrm{ini}}\geq 85\,M_{\odot}$) | $0.53\pm 0.05$ | $15\pm 5$
Table 1: Best fit parameters for the piece-wise linear relation between TAMS
mass and initial mass of the models at solar and LMC metallicity using the
different OB mass-loss prescriptions. The parameters fulfill the relation
$M_{\mathrm{TAMS}}=aM_{\mathrm{ini}}+b$.

The two regimes apparent in Fig. 2 arise from intrinsic differences in the
mass-loss prescriptions and their validity limits. More massive stars have a
higher mass loss and thus a larger core-to-envelope mass ratio. Thus further
leads to a larger expansion during the MS and thus a cooler
$T_{\mathrm{eff}}$. Since there are multiple temperature thresholds which
increase the applied mass-loss rate during cooling, such as the bistability
jumps for Vin01 or the limit of the validity domain for Bes20, MS expansion
causes additional boosts in mass-loss rate for stars above the initial mass
corresponding to the regime break in Fig. 2. This means that the stars in the
high MS mass-loss regime undergo most of their MS mass loss close to the end
of the MS. For example, 64% of the mass loss experienced during the MS by the
solar $60\,M_{\odot}$ model using Vin01 (plotted in Fig. 1) happens after the
first bistability jump in the last 8% of the MS evolution time. Stars at lower
metallicity are generally hotter and thus these stars must have a higher
initial mass to reach the high mass-loss regime compared to solar metallicity.
It is worth noting that for stars at the highest end of the mass range, the
surface-temperature trend reverses because mass loss is so strong that it
actually begins to expose hotter inner layers of the star. A star with an
initial mass of more than $100\,M_{\odot}$ already experiences intrinsically
high mass loss for either of the OB mass-loss prescriptions without any boost
from cooler effective temperatures.

Finally, it is worth drawing attention to another important effect of mass
loss during the MS, namely the reduction of the central convective region.
Stars with stronger mass loss have less massive helium cores by the end of
their MS, implying that MS winds do not only remove hydrogen directly from the
surface, but also prevent hydrogen within the star from being converted to
helium. In other words, mass loss is effectively a damper of core
nucleosynthesis despite being a surface process.

###  4.3 Stellar Yields

There is much interest in exploring the chemical output of massive stars,
given their unique ability to synthesize elements heavier than oxygen and
expel them into the interstellar medium through their winds and their final SN
explosion. While the aim of this paper is not to conduct a full study on the
effect of mass loss on stellar yields, we present of some important
implications of our results for yield calculations.

Figure 3: Total yields of carbon (left), nitrogen (center), and oxygen (right)
for models using different MS mass-loss prescriptions (colors and shapes) at
solar metallicity (large) and LMC metallicity (small). Figure 4: Relative
yield (i.e. the ratio of total yield and total mass loss) of oxygen versus
total mass loss (left), MS mass loss (center) and post-MS mass loss (right)
for models using different MS mass-loss prescriptions (colors and shapes) at
solar metallicity (large) and LMC metallicity (small).

Fig. 3 shows the total yield of CNO elements for each of our models,
calculated via

$\mathrm{Yield}_{X}=\sum
X_{\mathrm{surf,i}}\,(M_{\mathrm{i+1}}-M_{i})\text{,}$ (12)

where $M_{i}$ is the total mass of the star and $X_{\mathrm{surf,i}}$ is the
surface mass fraction of element $X$ at evolution timestep $i$.

As expected, the results show a significant impact of metallicity on the
yields of these elements. The separation by metallicity is especially apparent
for the nitrogen yield, which is mainly dependent on the initial abundance of
CNO elements and enhanced by the CNO cycle. For carbon and oxygen yields, we
see a more prominent influence of the chosen MS mass-loss prescription than
for the nitrogen yields. Since carbon and oxygen are expelled primarily during
the late evolution of the star when the CO core becomes exposed, this
illustrates how mass loss during the MS is able to impact the chemical
structure of a star much later during the evolution. This chain of causality
is indeed quite complex as there are multiple factors to consider. Firstly,
the yields of carbon and oxygen are affected by the efficiency of the
${}^{12}\text{C}(\alpha,\gamma)\,^{16}\text{O}$ reaction during helium
burning, which is set by the core temperature during helium burning, which in
turn is mainly a result of the mass of the star at the end of the MS.
Secondly, the transport of carbon and oxygen towards the surface depends on
the timing of mass loss with respect to the interior evolution. If mass loss
is weaker, core material has more time to be processed, but might not reach
the surface before the end of the evolution and therefore not contribute to
yields by winds. Conversely, strong mass loss might succeed in stripping the
star down further, but the material will then not be enriched as much with
carbon and/or oxygen. One may be tempted to try to relate the CO yields
directly either to the strong post-MS mass loss (which strips the star), or to
the weaker but longer MS mass loss (which sets the TAMS mass). However, as it
turns out, the total yield of carbon and oxygen correlates neither with the
amount of mass lost in the MS, nor with the mass lost after the MS, nor with
the total mass lost, as shown in Fig. 4 for oxygen. Instead, the detailed mass
loss history throughout all phases of evolution must be considered in order to
understand pre-collapse yields of massive stars. Mass loss during the MS plays
a pivotal role through its ability to determine the overarching mass loss
history of a star by influencing the duration and sequence of its post-MS
evolutionary phases. This will be explored further in Sect. 6.

##  5 Implications for the evolution of Wolf-Rayet stars

###  5.1 Spectral classification in evolution models

A WR star is defined as a spectral type characterized by strong and broad
emission lines, which is an indicator of a strong, optically thick wind that
obscures the hydrostatic layers of the star (Abbott & Conti 1987). Stellar
evolution models include only a simple approximated treatment of the stellar
atmosphere, and thus offer no way of determining a spectral type directly.
Therefore, stellar parameters such as surface abundances and temperature serve
as proxies to spectroscopic definitions of the various classes of stars,
although it is often unclear how closely this overlaps with observations.

In this work, we use these common criteria to study the consequences on
stellar evolution and the perceived populations. Our definitions of the
various types of stars are given in Table 2. It should be noted that the WC/WO
transition has been redefined since the analysis of WR stars by Georgy et al.
(2012) from an abundance to a temperature criterion, due to new results from
spectral modeling on evolutionary tracks by Groh et al. (2014). This has
further been supported by recent analyses of WC and WO stars by Aadland et al.
(2022). Moreover, we have renamed the subtypes WNL and WNE in Georgy et al.
(2012) to WNH and WN, respectively, in order to convey more directly that the
differentiating parameter between these two subtypes is the surface hydrogen
abundance, which does not necessarily correlate with the star’s position in
the HRD.

| Type | $\log(T_{\mathrm{eff}}\,[\mathrm{K}])$ | $X_{\mathrm{H}}$ |
---|---|---|---|---
| OB | $\geq 4$ | $>0.3$ | MS
| BSG | Post-MS
| YSG |  $3.7$ – $4$ |  |
| RSG | $\leq 3.7$ |  |
WR | WNH | $\geq 4$ |  $10^{-5}$ – $0.3$ |
WN | $\leq 10^{-5}$ | $\mathrm{C}\leq\mathrm{N}$
WC |  $4$ – $5.25$ | $\mathrm{C}>\mathrm{N}$
WO | $\geq 5.25$
Table 2: Definitions of the types of stars used in this work.

###  5.2 Surface hydrogen depletion

Figure 5: Evolutionary tracks of the $60\,M_{\odot}$ models at solar
metallicity, showing the surface hydrogen mass fraction as a function of
current total mass (hydrogen depletion curve). The horizontal line indicates
the threshold of surface hydrogen below which a hot star becomes classified as
a WR in our models. The “/” markers on the tracks separate the curves into the
four characteristic sections (I)–(IV), explained in the main text.

The WR definition based on low surface hydrogen abundance and high surface
temperature intrinsically links them to a formation scenario involving
significant mass loss, since hot and hydrogen-poor layers tend to be exposed
by stripping. Let us therefore consider more closely the evolution of the
surface hydrogen abundance $X_{\mathrm{s}}$ as a function of total mass loss
$\Delta M$, and for simplicity refer to this relation as the hydrogen
depletion curve of a given stellar model.

In Fig. 5, we show this curve for our $60\,M_{\odot}$ solar metallicity
models. As introduced in Sect. 4, the Bes20 model of this mass falls into the
low MS mass-loss regime, whereas the Vin01 model corresponds to high MS mass
loss (See Fig. 2). The curves can be divided into four phases, labeled as (I)
to (IV). During phase (I), the star is stripping its chemically homogeneous
envelope, and thus the surface hydrogen abundance remains constant. Phase (II)
marks a gradual decline of surface hydrogen abundance, caused either by
convective dredge-up after hydrogen-burning, or – if the mass loss happens
quickly enough – by the exposure of the ZAMS convective core that has since
shrunk but previously burned some of this region’s hydrogen. After the end of
the MS, when hydrogen shell-burning ignites around the core, a convective
region is established above this shell which removes the chemical gradient in
these layers. This region is sometimes referred to as the “intermediate
convective zone” (ICZ) in the literature, in reference to its position above
the convective core but below the layers that eventually become the convective
envelope in supergiants. When the ICZ is exposed, the hydrogen depletion curve
enters phase (III). Finally, when the hydrogen-burning shell is lost, the
sharp chemical gradient below this shell causes the surface hydrogen abundance
to drop to zero almost instantaneously and the star becomes a stripped helium
star (phase IV).

The fact that the two models in Fig. 5 yield different curves shows that there
must be some feedback mechanism linking MS mass loss to the chemical structure
of the interior layers before they are brought to the surface. This affects in
particular the resulting length and hydrogen abundance of phase (III), i.e.,
the ICZ. If MS mass loss is high, the gas pressure and temperature in this
region are low when hydrogen ignites in the shell since the envelope above it
is smaller at that moment. Therefore, the ICZ will be smaller and more
hydrogen-poor due to its proximity to the core (cf. Vin01 in Fig. 5).
Conversely, if the mass-loss rate during the MS is low, a large ICZ is
established and is able to transfer hydrogen-rich material from higher layers
into the interior of the star (cf. Bes20 in Fig. 5).

###  5.3 Evolution of WR stars

Figure 6: Effective temperature and evolution time as a function of initial
mass at the point of WR formation, i.e. at the first computational time-step
where each model fulfills our evolutionary criteria for WR stars. The colors
and shapes indicate the chosen OB mass-loss prescription as either Vin01
(yellow squares) or Bes20 (red circles). The results are shown in separate
panels depending on metallicity: solar (left) and LMC (right). The y-axis of
the bottom plots refers to the time-span elapsed since the ZAMS at the moment
of WR formation, normalized by the time spent on the MS by each model. The
purple dotted line represents the separation of the models into two evolution
scenarios, described in the text: the “early-formed” and “late-formed” WR
stars are below and above this line, respectively. Figure 7: Evolutionary
tracks in the Hertzsprung-Russell diagram (HRD) for the $30\,M_{\odot}$,
$60\,M_{\odot}$, and $85\,M_{\odot}$ models (from top to bottom) at solar
metallicity (left) and LMC metallicity (right). The shaded red area indicates
the domain of the red supergiant (RSG) phase
($\log(T_{\mathrm{eff}}\,[\mathrm{K}])<3.7$), and the shaded gray area
indicates the domain of the yellow supergiant (YSG) phase
($3.7\leq\log(T_{\mathrm{eff}}\,[\mathrm{K}])<4.0$), as defined in GENEC for
the application of specific mass-loss rates. The triangles mark time-steps of
$20\,000\,\mathrm{yr}$ starting from core helium burning. The onset of the WR
classification is marked by dotted squares and circles, for models using Vin01
and Bes20, respectively. The blue dashed line indicates the Humphreys-Davidson
limit (Davies et al. 2018; Lamers & Fitzpatrick 1988)

We now characterize the properties of our models at the moment they enter the
WR classification according to our defined criteria, to understand how stars
in single-star evolution and population synthesis models transition to the WR
phase. From the computations, all models at solar metallicity with an initial
mass of at least $30\,M_{\odot}$ reach the WR classification, regardless of
the applied MS mass-loss prescription. For LMC metallicity, the minimum
initial mass for WR formation is $40\,M_{\odot}$ using the Vin01 prescription,
and $50\,M_{\odot}$ with Bes20. In Fig. 6 an overview of the effective
temperature and the evolution time of the models at the entry into WR
classification is given as a function of initial mass, separated by
metallicity. The $y$-axis on the second row of this figure is the evolution
time (since the ZAMS) normalized by each model’s MS timescale. The MS
therefore ends at $t/\tau_{\mathrm{MS}}=1.0$, and helium burning (lasting
$\sim$$10\%$ of the MS timescale) ends at around $t/\tau_{\mathrm{MS}}=1.1$.

Fig. 6 essentially reveals two types of evolution towards WR classification,
which directly correspond to the two MS mass-loss regimes established in Sect.
4. The complete evolution of these two different classes of stars in the HRD
is shown for a few examples in Fig. 7, which also shows the observational
Humphreys-Davidson (HD) limit taken from Lamers & Fitzpatrick (1988), with the
cool LMC component taken as $\log L=5.5$ from Davies et al. (2018).

####  5.3.1 High MS mass-loss regime: “early-formed” WR stars

Stars in the high MS mass-loss regime are stripped very rapidly during the
cooler YSG phase, having already lost most of their envelope during the MS.
Since these stars have no hydrogen-rich convective zone, complete envelope
stripping is not required in order for the stellar surface to become hydrogen-
poor, and thus become classified as WR as soon as they cross the
$\log(T_{\mathrm{eff}})$ threshold of $4.0$. Since this happens very shortly
after the end of the MS (cf. Fig. 6), we refer to these WR stars as “early-
formed”. As they still retain a small amount of hydrogen on their surface,
they are first categorized as subtype WNH.

Of the models depicted in the HRD (Fig. 7), the “early-formed” (i.e., high MS
mass-loss regime) WR stars are the $85\,M_{\odot}$ models at solar metallicity
(both prescriptions) and LMC metallicity (Vin01 only), as well as the
$60\,M_{\odot}$ model at solar metallicity using Vin01.

In the HRD, such a star enters the WR phase while still having the radius and
luminosity of a supergiant, since it still has an inflated envelope at this
stage. As it loses more mass, the star gradually transitions to hotter
effective temperatures while reducing its luminosity. Once it has completely
lost its hydrogen, it becomes a stripped WN star with an effective temperature
hotter than $100\,000\,\mathrm{K}$. At this point, the star continues to strip
down completely, losing the helium envelope above the convective helium
burning core (which briefly expands just before being ejected causing a small
redward spike in the HRD), and part of the carbon-enriched interior (at which
point it would be classified as a WO/WC star). During this process, the star
remains at a relatively constant high temperature, but significantly reduces
its luminosity as it progressively reduces in size to less than one solar
radius. The dimming is halted by the ignition of helium shell burning, which
gives the star a burst of luminosity, further increasing its surface
temperature. This creates the hook-like feature at the end of the HRD tracks
which is characteristic of “early-formed” WR stars as well as some of the
earliest (most massive) late-formed WR stars (e.g., $60\,M_{\odot}$ LMC
models). The condition for this WR hook to appear is that the star needs to be
stripped of hydrogen before helium shell burning begins, i.e. the mass loss
needs to happen fast enough relative to the nuclear timescale. The more of the
helium envelope is then stripped before helium shell burning begins, the lower
the luminosity at the WR hook onset will be (e.g., compare both solar
$60\,M_{\odot}$ models).

To summarize, stars in the high MS mass-loss regime reach the WR stage early
and follow the evolutionary trajectory O–[YSG]–WNH–WN–WC/WO, where the YSG
phase is a very short transition period, if it exists at all. Since for these
stars, the onset of the WR stage (i.e., their hydrogen-depletion) coincides
with the beginning of core helium burning, they are close to the “Helium main
sequence”. These “Helium stars” are used as a proxy to study WR evolution
independent of the formation scenario (e.g., McClelland & Eldridge 2016; Yoon
2017). As we will discuss below, this proxy is only valid for the “early-
formed” subclass of WR stars.

####  5.3.2 Low MS mass-loss regime: “late-formed” WR stars

The stars in the low MS mass-loss regime spend at least the first half of
their helium burning time with a considerable envelope. This happens mostly in
the cool region of the HRD, where strong RSG/YSG mass loss gradually removes
the envelope. Once their envelope is thin enough, these stars then transition
to the blue region of the HRD. Since this envelope is hydrogen-rich due to the
large ICZ that formed after the TAMS, the stars continue their evolution being
classified as BSGs until their envelope is completely stripped, the surface
hydrogen abundance drops to zero, and the stars fulfill the criteria to be
classified as WN.

Of the depicted models in Fig. 7, the “late-formed” (low MS mass-loss regime)
WR stars are: the $85\,M_{\odot}$ LMC model using Bes20, both $60\,M_{\odot}$
LMC models, the $60\,M_{\odot}$ solar model using Bes20, and both
$30\,M_{\odot}$ solar models. The $30\,M_{\odot}$ models at LMC metallicity as
well as all $20\,M_{\odot}$ and $25\,M_{\odot}$ models lose too little mass to
ever be considered WR and these finish their lives as RSGs.

The “late-formed” WR stars follow the evolution path O–YSG/RSG–BSG–WN(–WC/WO),
formally “skipping” the WNH stage and instead evolving directly as a hydrogen-
free classical WR star with a more evolved core. Some of these stars can also
have a pronounced BSG phase just after the MS before expanding to become
YSGs/RSGs, which could be a result of our treatment of convection (See also
Sibony et al. 2023). Although these stars spend less time in the WR phase than
the “early-formed” WR stars, they are more numerous due to their lower initial
mass and therefore may also contribute significantly to observed and simulated
populations of WR stars.

##  6 Core-collapse progenitors and stellar populations

###  6.1 Evolutionary endpoints

Figure 8: HRD showing the evolutionary endpoints of the models that were
computed until the end of core carbon burning at solar metallicity (large
markers) and LMC metallicity (small markers). The labeled ellipses show
different populations of SN progenitors. Dashed black lines show equal radius
of several different values. The blue line indicates the HD limit. Figure 9:
Core mass at the evolutionary endpoints of the models that were computed until
end of core carbon burning at solar metallicity (large markers) and LMC
metallicity (small markers). The core region is defined as the layers with a
helium abundance below 1%.

We will now investigate the nature of the core-collapse progenitors by
examining the physical characteristics of the models at their endpoint.
Massive stars are commonly thought to end their lives as SNe, but recent
studies now show that a considerable number of massive stars do not produce
SNe but collapse silently or with a “failed” SN (Smartt 2009; Reynolds et al.
2015; Sukhbold et al. 2018; Beasor et al. 2023). We computed (most of) the
models until the end of core carbon burning, which we define as the point at
which the central carbon mass fraction drops below $10^{-5}$. This happens at
most a few years before core collapse. Considering the location of our models
in the HRD, as shown in Fig. 8, we see that our models split clearly into
several populations of possible core-collapse progenitors:

$M_{\mathrm{ini}}/M_{\odot}$ | Endpoint
---|---
| $Z=0.014$ | $Z=0.006$
| Vin01 | Bes20 | Vin01 | Bes20
20 | RSG
25 | – | YSG | RSG
30 | WN | WNH | RSG
40 | – | WNH | RSG
50 | WO | WN
60 | WO
66 | WO
73 | WO
80 | WO | WRSG | WC
85 | WO | WRSG | WC
95 | WO | WC | WRSG
105 | WO | WRSG
120 | WO | WRSG
Table 3: Estimated spectral type of the models at their evolutionary endpoint
based on stellar parameters.

  * •

Compact WR stars: These are the endpoints of most of the models with strong
enough mass loss. By the time of their death, most these stars have undergone
so much mass loss that their surfaces are mainly composed of carbon and oxygen
with only $\sim$20% helium. These WC/WO stars have removed almost all of the
envelope above their helium-burning shell, making them extremely hot and
bright. The exception applies to the $50\,M_{\odot}$ LMC models, which are
still formally classified as WN at the end of their evolution, but are also
extremely hot ($\log(T_{\mathrm{eff}}\,[\mathrm{K}])>5.25$).

  * •

“WR supergiants” (WRSG): This group contains the stars with the highest masses
but comparatively weak mass loss, i.e. solar metallicity $105\,M_{\odot}$ and
$120\,M_{\odot}$ using Bes20, as well as several high-mass LMC models. These
stars have very massive cores which enable them to inflate their helium
envelope after the onset of helium shell burning. They have the size and
luminosity of supergiants while also having characteristics of WR stars such
as the absence of hydrogen and enrichment of carbon and oxygen on the surface.
They populate a region considerably beyond the Humphreys-Davidson limit in the
HR diagram, and move there only within the last few years of their life,
making their observation very unlikely if they do exist. However, one notable
example of a star in this region of the HRD is $\eta$ Car (Hillier et al.
2001). It is presently unclear whether the derived HRD positions can be
physically reached through the WR phase or whether this is just a consequence
of the atmosphere treatment in GENEC.

  * •

WN stars: These are the least massive WR stars in our sample (of those models
that could be computed until the end), i.e. the $30\,M_{\odot}$ solar models
as well as the $40\,M_{\odot}$ LMC model with Vin01. These stars have just
lost their hydrogen envelope, but still have an envelope with a surface made
of the products of core hydrogen-burning, composed of mostly helium ($\sim
98$–$99$%) and some nitrogen. Their helium-burning shell is still
substantially covered, resulting in a cooler surface than the compact WR
stars.

  * •

RSG stars: These are the endpoints of the least massive stars, whose mass loss
is too weak to significantly expose the hydrogen-burning shell and evolve back
to the blue. They die as cool supergiants with only a slightly hydrogen-
depleted surface which still consists of $50$–$65$% hydrogen. The upper
luminosity limit of RSGs in our models lies at around $\log(L/L_{\odot})=5.8$,
significantly higher than the upper limit observed among RSG progenitors of
SNe, which lies closer to 5.2–5.3 (Davies & Beasor 2020).

The presumed spectral types for the endpoint of each model are summarized in
Table 3.

It is also interesting to consider the mass of the core at the end of the
computed evolution, since this quantity can be an important factor in driving
the final SN explosion and determine the type and mass of the remnant left
behind. The final CO core mass (defined as the mass of the region with a
combined carbon and oxygen mass fraction above $95\%$) is shown in Fig. 9 for
the models of different metallicities and using different MS mass-loss
prescriptions. It is evident that both metallicity and the MS mass-loss rate
have a great influence on the final core mass of the star, especially for
stars initially more massive than $50\,M_{\odot}$. The trend illustrates that
stars at lower metallicity end with more massive cores (due to their higher
compactness and lower mass loss). At any metallicity, stars with a weaker MS
mass-loss rate also end with larger cores (with the exception of the
$80\,M_{\odot}$ and $85\,M_{\odot}$ LMC models).

###  6.2 Timescales and populations

Figure 10: Post-MS mass loss by the end of core carbon burning of the models
at solar metallicity (large markers) and LMC metallicity (small markers). The
labeled regions denote which phase of mass loss is the dominant contributor to
the total mass lost after the MS. The division of models into RSG/YSG and WR
dominant post-MS mass loss corresponds exactly to the division into the low
and high MS mass-loss regime, respectively. Figure 11: Mass lost after the MS
versus during the MS for models at solar metallicity (large markers) and LMC
metallicity (small markers). The data is fitted with a $n=5$ polynomial fit,
which agrees with the data within 21.2% (root mean square relative deviation).
Visually, this illustrates that the mass loss experienced during the MS of
massive star evolution is a reasonable predictor of the mass loss that will be
experienced in the later stages of evolution, regardless of metallicity or MS
mass-loss prescription. In particular, the shape of this curve illustrates the
“YSG bump” for stars that lose around $10\,M_{\odot}$ during the MS, which
spend the most time in the YSG phase out of any of the models. Figure 12:
Duration of the various post-MS phases (colors), for each grid of models using
the MS mass-loss prescriptions by Vin01 (top row) and Bes20 (bottom row), at
solar metallicity (left column) and LMC metallicity (right column). The dashed
white vertical lines represent the approximate location of the break between
the low and high MS mass-loss regimes. The durations are shown cumulatively,
so there is no information about the order in which the different phases are
traversed. Furthermore, the $x$-axis is a power scale (where $x\sim
M_{\mathrm{ini}}^{-1.35}$), designed in such a way that the area of a region
corresponding to one phase is proportional to the number of stars currently in
that phase, assuming a population of constant star formation rate using the
standard initial mass function by Salpeter (1955).

As we discussed in Sect. 5, our evolution models follow a different trajectory
based on which MS mass-loss regime they are in. The sequence and duration of
the post-MS evolutionary phases is then a main determiner of the overall mass
loss history of a star, since each evolutionary phase is characterized by a
specific mass-loss rate. The results regarding the stellar yields (see Sect.
4.3) revealed the importance of the detailed mass loss history versus simple
total mass loss estimates. In Fig. 10, we consider the mass lost over all
post-MS phases of the models. Our models can clearly be divided into groups
depending on which post-MS phase is the most dominant contributor to the total
mass loss of the star. As it turns out, this division into WR-dominated and
RSG/YSG-dominated mass loss corresponds exactly to the division into high and
low MS mass-loss regimes, respectively, which establishes an interesting
causal link between MS mass loss and post-MS mass loss. Indeed, Fig. 11 shows
that within our framework, the mass loss experienced during the MS is a good
predictor of the mass loss that will be experienced after the MS, as
underlined by the goodness of the $n=5$ polynomial fit to the data:

$\begin{split}\Delta
M_{\mathrm{postMS}}/10M_{\odot}&=0.039x^{5}-0.62x^{4}+3.6x^{3}\\\
&-9.1x^{2}+9.1x+0.76\text{,}\end{split}$ (13)

where $x=\Delta M_{\mathrm{MS}}/10M_{\odot}$.

Expectedly, RSG mass loss is most significant for the stars of the lowest mass
in our sample. Stars with masses between $50\,M_{\odot}$ and $85\,M_{\odot}$
are dominated by YSG mass loss, if they are in the low MS mass-loss regime.
Otherwise, the YSG phase is too short, and the post-MS evolution is almost
exclusively spent in the WR phase. Fig. 10 also shows that the stars with the
largest post-MS mass loss are not the most massive stars, but the stars with
the longest YSG phase, since in our models, this phase has the strongest mass
loss (de Jager et al. 1988). This effect can also be seen as a “YSG bump” in
Fig. 11. An interesting consequence of this is that at the same initial mass
(e.g., $80\,M_{\odot}$), a star at lower metallicity may lose more mass in
total due to it spending more time in the YSG phase than a star which has a
higher metallicity and consequently loses more mass while still on the MS, but
less mass afterwards. Unfortunately, mass loss in the YSG region of the HRD is
still among the most poorly understood.

To understand the effect of MS mass loss on the timescales of various phases
more closely, consider Fig. 12 which shows stacked line graphs of the duration
of the post-MS phases obtained from the models. The results nicely summarize
some conclusions of the previous analysis. Firstly, regardless of the MS mass-
loss prescription or metallicity, stars in the low MS mass-loss regime spend
more time as YSGs, during which time they are removing their large hydrogen-
rich envelope. Secondly, these stars also spend less time as WNH, since they
transition quickly from BSG to hydrogen-poor WN when they remove their
hydrogen-rich ICZ. Comparing the different mass-loss prescriptions with each
other, we notice that the main difference lies in the mass limit between the
high- and the low MS mass-loss regime, while the overall trends within the
regimes remain the same. Finally, the figure also shows that higher
metallicity environments seem to favor RSGs, which is partly due to the fact
that the higher core-to-envelope ratios at the end of the MS cause these stars
to expand faster through the BSG and YSG stages and reach the RSG stage
earlier in their evolution.

From the timescales of the various phases, one can further infer the
distribution of massive stars in a population, assuming a constant star
formation rate as well as some initial mass function. The number of stars
$N_{i}$ in the population predicted to be currently in a certain evolutionary
phase $i$ is then directly proportional to the integral of the product of the
initial mass function $\mathrm{IMF}(M_{\mathrm{ini}})$ and the duration of
that phase $\tau_{i}(M_{\mathrm{ini}})$ over the initial mass range of the
population:

$N_{i}=K\int\tau_{i}(M_{\mathrm{ini}})\;\mathrm{IMF}(M_{\mathrm{ini}})\;\mathrm{d}M_{\mathrm{ini}}\text{,}$
(14)

where $K$ is some proportionality constant that normalizes the IMF to a given
star formation rate. Using
$\mathrm{IMF}(M_{\mathrm{ini}})=M_{\mathrm{ini}}^{-2.35}$ from Salpeter
(1955), we can absorb this term into a new $x$ coordinate with
$x=M_{\mathrm{ini}}^{-1.35}$, and the result is:

$N_{i}\propto\int\tau(M_{\mathrm{ini}})\;\mathrm{d}x\text{,}$ (15)

i.e., the number of stars $N_{i}$ is proportional to the area of the region
associated with phase $i$ in Fig. 12, creating a nice visual representation of
populations of massive stars with different MS mass-loss rates and different
metallicities.

| Relative abundance [%]
---|---
Phase | $Z=0.014$ | $Z=0.006$
| Vin01 | Bes20 | Vin01 | Bes20
MS | 89.9 | 89.8 | 90.1 | 90.2
RSG | 3.2 | 3.0 | 1.0 | 0.9
YSG | 1.8 | 3.5 | 2.0 | 2.1
BSG | 3.5 | 2.4 | 6.1 | 6.4
WR | 1.6 | 1.3 | 0.8 | 0.4
WNH | 0.3 | 0.5 | 0.4 | 0.3
WN | 0.2 | 0.2 | 0.1 | 0.1
WC | 1.1 | 0.6 | 0.3 | 0.1
WO | $<0.05$
Table 4: Relative abundances of the different phases in populations of massive
stars between $20\,M_{\odot}$ and $120\,M_{\odot}$, assuming constant star
formation rate, a power-law IMF from Salpeter (1955), and different MS mass-
loss rates and different metallicities. The values are numerical results
rounded to one decimal place.

Computing this integral numerically from the models (using trapezoid
integration), and normalizing to $100\%$, we obtain the distributions given in
Table 4. The results show a few interesting properties, which we list here in
no particular order:

  * •

There are fewer BSGs produced at high metallicities, due to the tendency of
hot evolved stars to be classified as WR instead. In the lower mass range,
higher metallicity stars expand faster than at lower metallicity,
transitioning through more quickly to the RSG phase after the MS.

  * •

There are fewer RSGs in lower metallicity environments, because these stars
expand more slowly off the MS. Additionally, RSGs at lower metallicity cannot
be progenitors of WR stars.

  * •

There are fewer WR stars in lower metallicity environments, with a higher
proportion of them being WNH.

##  7 Discussion

The results obtained in this study are based on evolutionary models of single
massive stars at solar and LMC metallicity and without rotation. We now
discuss the applicability of the results on domains that were not directly
explored in the simulations, and evaluate the limitations and the impact of
this study and its conclusions in the context of stellar evolution.

###  7.1 Metallicity

The presence of metals in the atmospheres of stars has the general effect of
increasing the opacity. This induces larger radiation pressure and therefore
larger mass loss. Often, a metallicity dependence is included in mass-loss
prescriptions in the form of a power law $\dot{M}\sim Z^{a}$, as is the case
for Vin01 in our models, but also others such as Björklund et al. (2023). The
metallicity dependence of Bes20 is less direct. In this prescription, the mass
loss depends on the electron scattering Eddington parameter, with a force-
multiplier $\alpha$ that supposedly corrects for the contribution of metallic
lines to the total opacity. For implementation in GENEC, we have used this
prescription with parameters calibrated by Brands et al. (2022) on LMC stars.
These stars have a metallicity of around half of the Sun’s and therefore
exhibit a mass-loss rate that is lower than galactic OB stars by a factor of
around 1.4–1.8, assuming a metallicity scaling with an exponent of around
0.5–0.85 (Kudritzki & Puls 2000; Vink et al. 2001). Since we have directly
applied this prescription at solar metallicity without adjustment, we expect
the resulting mass-loss rate to be underestimated by this factor. However, we
do not expect this to impact the results significantly, since the transition
between the low and high MS mass-loss regime does not depend on the absolute
scaling of the mass-loss rate.

In general, we find the studied effects of mass loss to be qualitatively
similar at LMC and solar metallicities, with the mass limits between the low
and high MS mass-loss regimes simply being shifted at lower metallicity. It is
reasonable to suppose that for a further metallicity reduction, single star
evolution could no longer produce “early-formed” WR stars (due to the lack of
a high MS mass-loss regime). If the metallicity is reduced even further, WR
stars would no longer be a possible outcome of single star evolution at all.
At this point, binary interactions are expected to be the only remaining path
towards WR formation (e.g., Shenar et al. 2020b).

###  7.2 Binarity and rotation

Binarity has long been studied as an important aspect of massive star
evolution (Vanbeveren 1991; Podsiadlowski et al. 1992; Wellstein & Langer
1999; Mason et al. 2009). Most massive stars interact with a binary companion
at some point during their lives (Sana et al. 2012, 2013), and as many as 30%
of massive MS stars are the product of binary mergers (e.g., de Mink et al.
2014). The mass flow induced by binary interactions may dominate (by far) the
effect of radiatively driven winds in several regimes, complicating the
comparison of the results from the present study to real populations of
massive stars. However, the effect of mass loss or accretion via binary
interaction on the evolution of stars may be similar to those caused by
stellar winds and could in theory be interpreted in the same way. For example,
a massive star that experiences a binary stripping episode during the MS (case
A mass transfer) may experience a similar subsequent evolution as a star that
has lost the equivalent amount of mass through MS winds. Moreover, the mass
removal in cool regime can be treated as a proxy for the primaries in systems
that undergo case-B mass transfer, i.e., those transferring mass when the
primary expands towards the RSG regime.

In an effort to produce a synthetic WR population for the LMC, Pauli et al.
(2022) calculated both single and binary evolution models. In their
illustrated binary models with case-B mass transfer, the mass transfer
essentially removes a large part of the hydrogen-rich envelope. After the
significant mass removal, the models predict that the primary will spend most
of its lifetime as a hot WR, hydrogen-poor star, very similar to our cases
with high MS mass loss. Hence, the fate of such a primary marks a case of what
we call an “early formed” WR in our single-star models. Yet, the time scales
are not expected to be similar by default. Compared to the short time scale of
the binary mass transfer, the inner structure of the single star models might
be more evolved as even our single star models in the high MS mass-loss regime
will usually spend some time in the RSG and YSG stage, albeit this still being
very short compared to the total He-burning time (around $0.3$–$2.5$% of core
He-depletion). One aspect not mapped in this analogy is the possibility to
regain hydrogen, which exists in binary systems in the case of mass transfer
back from the secondary onto the primary. Hence, the level of the plateau in
the hydrogen depletion curve can also be higher for a post-interaction primary
than for single-star evolution. This would also contribute to a variety of
hydrogen surface abundances for partially stripped stars, such as those
observed by Ramachandran et al. (2023); Götberg et al. (2023), where some seem
to spend a considerable amount of their post-MS evolution at temperatures and
luminosities similar to main-sequence OB stars, albeit with much lower current
mass.

Binary interaction can also inﬂuence the rotational behavior of stars,
especially in the event of a merger. The surface anisotropy of rotating stars
introduces an entirely new layer of complexity into the treatment of mass
loss, even for single stars (e.g., Poe 1987; Bjorkman & Cassinelli 1993;
Müller & Vink 2014), which ideally requires modeling in (at least)
2-dimensions. Rotation also aﬀects other aspects of stellar evolution such as
convective mixing within the star (Maeder & Meynet 2000), which we identified
to be affected by prior mass loss with crucial implications for the subsequent
evolution.

As covering the whole implications from binarity and rotation is far beyond
the scope of the present study, the results from this work should therefore be
regarded as an exploration of the physical mechanisms relating mass loss and
stellar evolution.

###  7.3 Unexplored uncertainties in single star evolution models

With this study being focused on mass loss, various parameters describing
convection in 1D stellar models were not explored in this study. Yet, many of
our results follow from properties of convective zones in the star, in
particular the ICZ. ICZs form upon ignition of the hydrogen shell after the MS
in all of our models, albeit they are more prominent in stars that have
previously experienced lower mass loss (cf. Fig. 5). The appearance of ICZs is
also a consequence of our choice of convection parameters. For example, if one
varies the overshooting parameter or the semiconvective efficiency, one can
reasonably obtain models that do not form ICZs (Schootemeijer et al. 2019).
Additionally, our use of the Schwarzschild criterion for convection favors the
formation of larger convective regions as opposed to models which implement
the more restrictive Ledoux criterion (Sibony et al. 2023).

If we chose convection parameters which inhibit the formation of large ICZs,
more models would instead follow the “early-formed” WR formation scenario,
where WR classification would be reached before total hydrogen depletion. As a
result, models would predict a greater number of “cool” WNH stars. Moreover,
stars would spend more time in the RSG phase and perhaps even stars with
initial masses above $40\,M_{\odot}$ could become RSGs, since less efficient
convection favors redward evolution (Schootemeijer et al. 2019; Sibony et al.
2023).

Secondly, it is known that there are large uncertainties in the mass-loss
rates of various phases of post-MS evolution. These types of mass loss and
their impact on massive stars has been the subject of previous investigations,
for example Yang et al. (2023); Beasor et al. (2021); Massey et al. (2023) for
RSGs, Owocki et al. (2019); Smith & Owocki (2006) for eruptive LBV mass loss,
and Higgins et al. (2021); Sander et al. (2023) for WRs.

The uncertainties in convective mixing and supergiant/WR mass loss most
importantly affect the evolution after the MS. Since the post-MS treatment in
GENEC is fixed within the framework of our study, our results highlight the
ways that varying the mass-loss rate only during the MS can affect stars in
every stage of evolution and show that the uncertainties in MS mass loss are
indeed not negligible. MS mass loss must therefore be viewed as an important
initial condition in studies concentrating on post-MS massive stars.

###  7.4 Impact of MS mass loss and existence of the mass-loss regimes

Our results show a clear separation of the models into two MS mass-loss
regimes for both the Vin01 and the Bes20 mass-loss prescriptions, and many
subsequent conclusions follow this divide. These results are expected to
propagate into the predictions of stellar population synthesis codes, which
widely rely on prototypical stellar evolution models such as these ones. Based
on the evolution models, the choice of MS prescription mainly has an impact on
the initial mass at which the MS mass-loss regime break is located, and thus
most importantly affects stars between $55$–$75$ $M_{\odot}$ at solar
metallicity, and between $75$–$95$ $M_{\odot}$ at LMC metallicity. It is also
clear that without any MS mass loss, all stars would behave as if they were in
the low MS mass-loss regime (see 1:1 relation in Fig. 2) and thus would end in
RSGs or late-formed WRs. MS mass loss is therefore crucial for the evolution
of early-formed WR stars, and in particular can explain the occurrence of WR
stars with hydrogen (WNH).

One may question whether the clustering into the two MS mass-loss regimes
presented in Fig. 2 is real. These regimes are the result of boosts in the
mass-loss rate as a star expands while on the MS, yet it debated whether MS
stars truly behave this way. Björklund et al. (2023) argue with their
hydrodynamically consistent wind models for the absence of the iron
bistability jump found by Vink et al. (1999). However, due to the constraints
of their parameter space to stars smaller than $80\,M_{\odot}$, it is unclear
whether their models would also show a mass-loss regime change at higher
masses, in particular as the results are in conflict with theoretical
predictions by Krtička et al. (2021) and Bernini-Peron et al. (2023). The
$\dot{M}(\Gamma_{\mathrm{{e}}})$ dependencies in the Bes20 prescription are
motivated by observational studies such as Bestenlehner et al. (2014) and
Brands et al. (2022). While this recipe does not include a discontinuity in
itself, there is still a sharp regime change for the Bes20 models caused by
the star exiting the domain of validity with $T_{\mathrm{eff}}$ below 30 000 K
and entering the region of the HRD where we apply stronger YSG winds (see Fig.
1). The high mass-loss rate of YSGs in our models is a consequence of applying
the prescription by de Jager et al. (1988), which is a generic prescription
not designed specifically for YSGs and merely used as a “fall-back” since no
better prescription exists for this regime. However, YSGs are indeed
associated with significant mass loss (Koumpia et al. 2020; Bonanos et al.
2023; Kraus et al. 2023; Humphreys et al. 2023), which motivates the use of
this prescription.

###  7.5 Observability

In the study of massive stars, observational constraints of evolutionary
models are often provided in the form of population surveys, with key
indicators being the abundance ratios between different types of stars, e.g.
WR/O, BSG/RSG, and so on (see, e.g., Eggenberger et al. 2002; Neugent & Massey
2019; Higgins & Vink 2020; Wagle et al. 2020). We present a visual
representation of these abundance ratios in populations with constant star
formation rate in Fig. 12, with data in Table 4. Assuming no WR stars form
below the minimum mass of our model grid ($20\,M_{\odot}$), our predictions
can be compared fairly to properties of observed WR populations. For example,
we expect the WC/WN ratio to be approximately $2.2$ (MW, Vin01), $0.85$ (MW,
Bes20), $0.6$ (LMC, Vin01) and $0.25$ (LMC, Bes20) for our grids. The Bes20
populations fit remarkably well with the observations of Neugent & Massey
(2019), who find a WC/WN ratio of $0.83$ for the MW and $0.23$ for the LMC.
This indicates that weaker MS mass loss, as tested by our implementation of
Bes20, may be a viable approach in order to correctly reproduce the
distribution of WR stars in stellar populations. Of course, the other
uncertainties discussed in this paper may also have an impact on the WC/WN
ratio.

Moreover, we note the discrepancy between our modeled evolutionary tracks (see
Fig. 7) and the observational Humphreys-Davidson (HD) limit in the Galaxy and
the LMC (Lamers & Fitzpatrick 1988; Humphreys & Davidson 1984). Especially the
models using Bes20 spend much more time in the “forbidden” cool area of the
HRD, which suggests that this mass-loss rate is lower than reality, or that
there are other mechanisms involved in avoiding too much inflation, in
particular for the more massive models. To cite an extreme example, the $85\
M_{\odot}$ LMC model in Fig. 7 spends $\sim$ $250\,000$ years beyond the HD
limit, which represents most of its helium-burning lifetime. However, models
with initial masses as low as $40\,M_{\odot}$, including those from previous
work (e.g., Ekström et al. 2012), do also cross the HD limit. Gilkis et al.
(2021) has suggested that an appropriate choice of rotational velocities and
convective overshooting parameters can satisfy the HD limit, but not without
in turn creating an excess of hot supergiants. In general, the tendency of
evolutionary tracks to go beyond the HD limit is a common problem in evolution
modeling that is still a subject of active investigation (e.g., Higgins & Vink
2020).

As mentioned previously, the conclusions of this work should be taken as an
analysis of the physical link between mass loss and evolution. The effects
explored here could potentially be used as indicators for the validity of
certain mass-loss prescriptions, but this cannot be done without more complete
model grids, nor without accounting for other uncertainties in stellar models.
Observations could also be used to compare and calibrate the mass-loss
prescriptions directly so that their implementation in stellar evolution codes
is more robust.

###  7.6 Theoretical domain limits

For this work, we have defined different types of stars (e.g. RSG, YSG, WR,
etc.) via strict criteria on stellar parameters that are easily provided by a
stellar evolution code. Specifically, we have defined WR stars by a threshold
of surface temperature and hydrogen mass fraction. Based on this definition,
we have found that the classification of a WR star as “early-formed” or “late-
formed” depends directly on the hydrogen content of its ICZ with respect to
the threshold value of 0.3. If this shell is hydrogen-rich, the WR star forms
“late”, whereas if it is hydrogen-poor, the WR star forms “early”. In general,
we would reach similar conclusions regarding the formation and evolution of WR
stars if the threshold were chosen differently. Shifting the WR surface
hydrogen threshold to a higher value (e.g., 0.4 as it is sometimes defined,
Schootemeijer & Langer 2018) would simply lower the mass limit between early-
formed and late-formed WR stars. Early-formed WR stars would form earlier,
while the remaining late-formed WR stars would form at about the same moment,
since the drop from hydrogen-rich to hydrogen-depleted happens almost
instantaneously when the hydrogen-burning shell is removed and thus the exact
hydrogen threshold value does not matter.

Yet, it is questionable whether a strict hydrogen and temperature threshold is
a valid assumption at all. In reality, WR stars are a spectroscopic class,
independent of surface hydrogen abundance, and evolutionary criteria are a
simply proxy defined to match the spectroscopy as closely as possible.
Problematic aspects of this classification have already been identified, for
example the existence of hot hydrogen-poor stars with optically thin winds
whose mass loss is significantly lower than classical WN stars (Vink 2017;
Shenar et al. 2020a). These kinds of stars would be classified as WR by our
evolutionary criteria, although they show spectra that are distinct from WR
stars. It is part of future work to investigate the connections between
spectroscopic and evolutionary classification and establish a more consistent
picture of massive star evolution from both perspectives.

##  8 Summary and conclusions

In this paper we performed a systematic study of the effects of MS mass loss
on the evolution of massive stars in the range of $20\,M_{\odot}$ to
$120\,M_{\odot}$ at two different metallicities. Using the stellar evolution
code GENEC, we compared the effect of applying two different mass-loss
prescriptions for OB-type stars. For the direct impact on the MS evolution
itself, including the total MS mass loss, we obtain two distinct regimes of
low or high MS mass loss.

Mass loss during the MS acts early and over long timescales. The results show
that it affects the evolution in two ways: First, it leads directly to the
removal of the outer layers of the star, beginning with the hydrogen envelope.
This has the general consequence of exposing metals on the surface and
increasing the effective temperature once enough material has been removed. In
this regard, mass loss acts in the same way on the MS as during any other
phase. However, MS mass loss has shown here to have a second effect which is
unique to this phase and acts more indirectly, namely through its connection
with the formation and evolution of interior convective zones. Our results
show that mass loss not only affects the size of the convective core during
the MS, but more importantly, MS mass loss is almost the sole determinant of
the mass a star has when it inflates to the RSG stage, which arguably marks
the most dramatic structure change in the life of a massive star. It is right
before the inflation that our models form a hydrogen-burning shell with a
convective zone above it (known as the ICZ), and that the core contracts and
begins to burn helium in the center. The size and composition of the ICZ is
directly dependent on MS mass loss, and will go on to play an important role
in determining the formation and appearance of WR stars. Their evolution
across the HRD, which is mostly independent of our WR classification scheme,
is also influenced by MS mass loss.

After the inflation, mass loss loses its connection to the evolution of the
stellar interior and its effect becomes more straightforward, simply stripping
the star progressively. A useful visual representation of how MS mass loss
affects the interior structure of stars is the hydrogen depletion curve (e.g.,
Fig. 5), which shows a “fossilized” record of the hydrogen abundance of
interior layers as they were brought to the surface.

Our study shows that even during the MS, the effects of mass loss cannot be
ignored and must be accounted for in the uncertainties of any stellar model.
One important conclusion of this paper is the importance of mass loss history
for the evolution of massive stars. This is shown in particular by the results
regarding stellar yields, which do not seem to be correlated simply with
either total mass loss or mass loss during just one phase of evolution. The MS
plays a special role in determining the subsequent mass loss history of a
star, since the mass loss during this phase determines the sequence and
duration of the post-MS phases of evolution, each of which has its own
characteristic mass-loss rate. The models considered here, which are
prototypical of current population synthesis models, show the following
possible evolutionary paths, which are significantly influenced by MS mass
loss:

  * (a)

OB $\longrightarrow$ RSG ($\longrightarrow$ blue loop)

  * (b)

OB $\longrightarrow$ RSG/YSG $\longrightarrow$ BSG $\longrightarrow$ WN
($\longrightarrow$ WC/WO)

  * (c)

OB $\longrightarrow$ WNH $\longrightarrow$ WN $\longrightarrow$ WC/WO

Stars with the lowest initial masses in our grid experience type (a) evolution
due to their mass loss being too low to significantly strip the hydrogen
envelope. Some of these stars may attempt to cross the HRD again, shortly
becoming BSGs in a blue loop, but none will reach the WR phase. Type (b)
evolution corresponds to the models that have strong enough mass loss to form
WR stars, but are still located in the low MS mass-loss regime. As such, they
spend a considerable amount of time as cool RSGs/YSGs and then evolve back to
the blue and become hydrogen-depleted WN stars, and eventually nitrogen-
depleted WC/WO stars if the mass loss is strong enough. These stars only enter
the WR phase in the latter half of core helium burning, which we have referred
to as “late-formed” WR stars. Finally, models with type (c) tracks evolve to
the WR phase directly from the MS, only briefly passing through the YSG phase
for a few thousand years. These models become WR stars close to the beginning
of core helium-burning and we have therefore called them “early-formed” WR
stars, representing WNH stars while hydrogen still remains on the surface.
These evolutionary tracks correspond largely to the so-called “Conti scenario”
of massive star evolution (Conti 1975; Maeder 1996), but the mass limits are
very sensitive to the chosen MS mass-loss rate and metallicity.

Next to our detailed study of MS mass loss, we have also addressed other
uncertainties of the models in Sect. 7.3, namely the parameters describing
convection, which are critical for explaining the appearance of the ICZ.

Although we only model the evolution of individual stars, the obtained paths
towards WR classification can also be applied to binary systems with mass
transfer. In systems where mass transfer occurs when stars enter the
Hertzsprung gap, the mass loss onto the companion fulfills a similar role as
the high assumed mass loss in the RSG and YSG stages, albeit not necessarily
in a similar amount of time. Thus, the general distinction between stars
entering the WR stage early, i.e., at the beginning of He burning or even
sooner, or late remains valid even when accounting for binaries. To better
predict the individual populations, more detailed studies for specific
metallicities will be required.

###### Acknowledgements.

The authors thank the anonymous referee whose detailed and constructive
comments helped improve the quality of this paper considerably. The authors
acknowledge fruitful discussions made possible by the International Space
Science Institute (ISSI) in Bern, through ISSI International Team project 512
(Multiwavelength View on Massive Stars in the Era of Multimessenger Astronomy,
PI Oskinova). JJ acknowledges funding from the Deutsche Forschungsgemeinschaft
(DFG, German Research Foundation) Project-ID 496854903 (SA4064/2-1, PI Sander)
and is a member of the International Max Planck Research School for Astronomy
and Cosmic Physics at the University of Heidelberg (IMPRS-HD). SE acknowledges
the STAREX grant from the ERC Horizon 2020 research and innovation programme
(grant agreement No. 833925). AACS is supported by the Deutsche
Forschungsgemeinschaft (DFG - German Research Foundation) in the form of an
Emmy Noether Research Group – Project-ID 445674056 (SA4064/1-1, PI Sander) and
acknowledges funding from the Federal Ministry of Education and Research
(BMBF) and the Baden-Württemberg Ministry of Science as part of the Excellence
Strategy of the German Federal and State Governments.

## References

  * Aadland et al. (2022) Aadland, E., Massey, P., Hillier, D. J., et al. 2022, ApJ, 931, 157
  * Abbott & Conti (1987) Abbott, D. C. & Conti, P. S. 1987, ARA&A, 25, 113
  * Beasor & Davies (2018) Beasor, E. R. & Davies, B. 2018, MNRAS, 475, 55
  * Beasor et al. (2021) Beasor, E. R., Davies, B., & Smith, N. 2021, ApJ, 922, 55
  * Beasor et al. (2023) Beasor, E. R., Hosseinzadeh, G., Smith, N., et al. 2023, JWST reveals a luminous infrared source at the position of the failed supernova candidate N6946-BH1
  * Bernini-Peron et al. (2023) Bernini-Peron, M., Marcolino, W. L. F., Sander, A. A. C., et al. 2023, A&A, 677, A50
  * Bestenlehner (2020) Bestenlehner, J. M. 2020, MNRAS, 493, 3938
  * Bestenlehner et al. (2020) Bestenlehner, J. M., Crowther, P. A., Caballero-Nieves, S. M., et al. 2020, MNRAS, 499, 1918
  * Bestenlehner et al. (2014) Bestenlehner, J. M., Gräfener, G., Vink, J. S., et al. 2014, A&A, 570, A38
  * Björklund et al. (2023) Björklund, R., Sundqvist, J. O., Singh, S. M., Puls, J., & Najarro, F. 2023, A&A, 676, A109
  * Bjorkman & Cassinelli (1993) Bjorkman, J. E. & Cassinelli, J. P. 1993, ApJ, 409, 429
  * Bonanos et al. (2023) Bonanos, A. Z., Tramper, F., de Wit, S., et al. 2023, arXiv e-prints, arXiv:2312.04626
  * Brands et al. (2022) Brands, S. A., de Koter, A., Bestenlehner, J. M., et al. 2022, A&A, 663, A36
  * Castor et al. (1975) Castor, J. I., Abbott, D. C., & Klein, R. I. 1975, ApJ, 195, 157
  * Conti (1975) Conti, P. S. 1975, Memoires of the Societe Royale des Sciences de Liege, 9, 193
  * Crowther (2001) Crowther, P. A. 2001, in Astrophysics and Space Science Library, Vol. 264, The Influence of Binaries on Stellar Population Studies, ed. D. Vanbeveren, 215
  * Davies & Beasor (2020) Davies, B. & Beasor, E. R. 2020, MNRAS, 493, 468
  * Davies et al. (2018) Davies, B., Crowther, P. A., & Beasor, E. R. 2018, MNRAS, 478, 3138
  * de Jager et al. (1988) de Jager, C., Nieuwenhuijzen, H., & van der Hucht, K. A. 1988, A&AS, 72, 259
  * de Mink et al. (2014) de Mink, S. E., Sana, H., Langer, N., Izzard, R. G., & Schneider, F. R. N. 2014, ApJ, 782, 7
  * Eggenberger et al. (2021) Eggenberger, P., Ekström, S., Georgy, C., et al. 2021, A&A, 652, A137
  * Eggenberger et al. (2002) Eggenberger, P., Meynet, G., & Maeder, A. 2002, A&A, 386, 576
  * Eggenberger et al. (2008) Eggenberger, P., Meynet, G., Maeder, A., et al. 2008, Ap&SS, 316, 43
  * Ekström et al. (2012) Ekström, S., Georgy, C., Eggenberger, P., et al. 2012, A&A, 537, A146
  * Eldridge & Vink (2006) Eldridge, J. J. & Vink, J. S. 2006, A&A, 452, 295
  * Georgy et al. (2012) Georgy, C., Ekström, S., Meynet, G., et al. 2012, A&A, 542, A29
  * Gilkis et al. (2021) Gilkis, A., Shenar, T., Ramachandran, V., et al. 2021, MNRAS, 503, 1884
  * Götberg et al. (2023) Götberg, Y., Drout, M. R., Ji, A. P., et al. 2023, ApJ, 959, 125
  * Gräfener & Hamann (2008) Gräfener, G. & Hamann, W. R. 2008, A&A, 482, 945
  * Groh et al. (2014) Groh, J. H., Meynet, G., Ekström, S., & Georgy, C. 2014, A&A, 564, A30
  * Higgins et al. (2021) Higgins, E. R., Sander, A. A. C., Vink, J. S., & Hirschi, R. 2021, MNRAS, 505, 4874
  * Higgins & Vink (2020) Higgins, E. R. & Vink, J. S. 2020, A&A, 635, A175
  * Hillier et al. (2001) Hillier, D. J., Davidson, K., Ishibashi, K., & Gull, T. 2001, ApJ, 553, 837
  * Humphreys & Davidson (1984) Humphreys, R. M. & Davidson, K. 1984, Science, 223, 243
  * Humphreys et al. (2023) Humphreys, R. M., Jones, T. J., & Martin, J. C. 2023, AJ, 166, 50
  * Koumpia et al. (2020) Koumpia, E., Oudmaijer, R. D., Graham, V., et al. 2020, A&A, 635, A183
  * Kraus et al. (2023) Kraus, M., Kourniotis, M., Arias, M. L., Torres, A. F., & Nickeler, D. H. 2023, Galaxies, 11, 76
  * Krause et al. (2013) Krause, M., Fierlinger, K., Diehl, R., et al. 2013, A&A, 550, A49
  * Krtička et al. (2021) Krtička, J., Kubát, J., & Krtičková, I. 2021, A&A, 647, A28
  * Krtička et al. (2024) Krtička, J., Kubát, J., & Krtičková, I. 2024, A&A, 681, A29
  * Kudritzki & Puls (2000) Kudritzki, R.-P. & Puls, J. 2000, ARA&A, 38, 613
  * Lamers & Fitzpatrick (1988) Lamers, H. J. G. L. M. & Fitzpatrick, E. L. 1988, ApJ, 324, 279
  * Lamers et al. (1995) Lamers, H. J. G. L. M., Snow, T. P., & Lindholm, D. M. 1995, ApJ, 455, 269
  * Maeder (1996) Maeder, A. 1996, in Liege International Astrophysical Colloquia, Vol. 33, Liege International Astrophysical Colloquia, ed. J. M. Vreux, A. Detal, D. Fraipont-Caro, E. Gosset, & G. Rauw, 39
  * Maeder & Meynet (2000) Maeder, A. & Meynet, G. 2000, ARA&A, 38, 143
  * Mason et al. (2009) Mason, B. D., Hartkopf, W. I., Gies, D. R., Henry, T. J., & Helsel, J. W. 2009, AJ, 137, 3358
  * Massey et al. (2023) Massey, P., Neugent, K. F., Ekström, S., Georgy, C., & Meynet, G. 2023, ApJ, 942, 69
  * Massey et al. (2014) Massey, P., Neugent, K. F., Morrell, N., & Hillier, D. J. 2014, ApJ, 788, 83
  * Mauerhan et al. (2011) Mauerhan, J. C., Van Dyk, S. D., & Morris, P. W. 2011, AJ, 142, 40
  * McClelland & Eldridge (2016) McClelland, L. A. S. & Eldridge, J. J. 2016, Monthly Notices of the Royal Astronomical Society, 459, 1505
  * Müller & Vink (2014) Müller, P. E. & Vink, J. S. 2014, A&A, 564, A57
  * Neugent & Massey (2019) Neugent, K. & Massey, P. 2019, Galaxies, 7, 74
  * Nugis & Lamers (2000) Nugis, T. & Lamers, H. J. G. L. M. 2000, A&A, 360, 227
  * Owocki et al. (2019) Owocki, S. P., Hirai, R., Podsiadlowski, P., & Schneider, F. R. N. 2019, MNRAS, 485, 988
  * Pauldrach & Puls (1990) Pauldrach, A. W. A. & Puls, J. 1990, A&A, 237, 409
  * Pauli et al. (2022) Pauli, D., Langer, N., Aguilera-Dena, D. R., Wang, C., & Marchant, P. 2022, A&A, 667, A58
  * Petrov et al. (2016) Petrov, B., Vink, J. S., & Gräfener, G. 2016, MNRAS, 458, 1999
  * Podsiadlowski et al. (1992) Podsiadlowski, P., Joss, P. C., & Hsu, J. J. L. 1992, ApJ, 391, 246
  * Poe (1987) Poe, C. H. 1987, PhD thesis, University of Wisconsin, Madison
  * Ramachandran et al. (2023) Ramachandran, V., Klencki, J., Sander, A. A. C., et al. 2023, A&A, 674, L12
  * Reynolds et al. (2015) Reynolds, T. M., Fraser, M., & Gilmore, G. 2015, MNRAS, 453, 2885
  * Rosslowe & Crowther (2015) Rosslowe, C. K. & Crowther, P. A. 2015, Monthly Notices of the Royal Astronomical Society, 447, 2322
  * Sabhahit et al. (2022) Sabhahit, G. N., Vink, J. S., Higgins, E. R., & Sander, A. A. C. 2022, MNRAS, 514, 3736
  * Salpeter (1955) Salpeter, E. E. 1955, ApJ, 121, 161
  * Sana et al. (2013) Sana, H., de Koter, A., de Mink, S. E., et al. 2013, A&A, 550, A107
  * Sana et al. (2012) Sana, H., de Mink, S. E., de Koter, A., et al. 2012, Science, 337, 444
  * Sander et al. (2023) Sander, A. A. C., Lefever, R. R., Poniatowski, L. G., et al. 2023, A&A, 670, A83
  * Sander et al. (2020) Sander, A. A. C., Vink, J. S., & Hamann, W. R. 2020, MNRAS, 491, 4406
  * Sander et al. (2022) Sander, A. A. C., Vink, J. S., Higgins, E. R., et al. 2022, in The Origin of Outflows in Evolved Stars, ed. L. Decin, A. Zijlstra, & C. Gielen, Vol. 366, 21–26
  * Schootemeijer & Langer (2018) Schootemeijer, A. & Langer, N. 2018, A&A, 611, A75
  * Schootemeijer et al. (2019) Schootemeijer, A., Langer, N., Grin, N. J., & Wang, C. 2019, A&A, 625, A132
  * Shenar et al. (2020a) Shenar, T., Bodensteiner, J., Abdul-Masih, M., et al. 2020a, A&A, 639, L6
  * Shenar et al. (2020b) Shenar, T., Sablowski, D. P., Hainich, R., et al. 2020b, A&A, 641, C2
  * Sibony et al. (2023) Sibony, Y., Georgy, C., Ekström, S., & Meynet, G. 2023, A&A, 680, A101
  * Smartt (2009) Smartt, S. J. 2009, ARA&A, 47, 63
  * Smith & Owocki (2006) Smith, N. & Owocki, S. P. 2006, ApJ, 645, L45
  * Sukhbold et al. (2018) Sukhbold, T., Woosley, S. E., & Heger, A. 2018, ApJ, 860, 93
  * Vanbeveren (1991) Vanbeveren, D. 1991, Space Sci. Rev., 56, 249
  * Vink (2017) Vink, J. S. 2017, A&A, 607, L8
  * Vink (2022) Vink, J. S. 2022, ARA&A, 60, 203
  * Vink & de Koter (2005) Vink, J. S. & de Koter, A. 2005, A&A, 442, 587
  * Vink et al. (1999) Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 1999, A&A, 350, 181
  * Vink et al. (2001) Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2001, A&A, 369, 574
  * Vink et al. (2011) Vink, J. S., Muijres, L. E., Anthonisse, B., et al. 2011, A&A, 531, A132
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261
  * Wagle et al. (2020) Wagle, G. A., Ray, A., & Raghu, A. 2020, ApJ, 894, 118
  * Wellstein & Langer (1999) Wellstein, S. & Langer, N. 1999, A&A, 350, 148
  * Yang et al. (2023) Yang, M., Bonanos, A. Z., Jiang, B., et al. 2023, A&A, 676, A84
  * Yoon (2017) Yoon, S.-C. 2017, MNRAS, 470, 3970