# The cosmic carbon footprint of massive stars stripped in binary systems

R. Farmer<EMAIL_ADDRESS>Max-Planck-Institut für Astrophysik, Karl-
Schwarzschild-Straße 1, 85741 Garching, Germany  Anton Pannekoek Institute for
Astronomy and GRAPPA, University of Amsterdam, NL-1090 GE Amsterdam, The
Netherlands  E. Laplace Anton Pannekoek Institute for Astronomy and GRAPPA,
University of Amsterdam, NL-1090 GE Amsterdam, The Netherlands  S. E. de Mink
Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Straße 1, 85741
Garching, Germany  Anton Pannekoek Institute for Astronomy and GRAPPA,
University of Amsterdam, NL-1090 GE Amsterdam, The Netherlands  Center for
Astrophysics — Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138,
USA  S. Justham School of Astronomy & Space Science, University of the Chinese
Academy of Sciences, Beijing 100012, China  National Astronomical
Observatories, Chinese Academy of Sciences, Beijing 100012, China  Anton
Pannekoek Institute for Astronomy and GRAPPA, University of Amsterdam, NL-1090
GE Amsterdam, The Netherlands  Max-Planck-Institut für Astrophysik, Karl-
Schwarzschild-Straße 1, 85741 Garching, Germany

###### Abstract

The cosmic origin of carbon, a fundamental building block of life, is still
uncertain. Yield predictions for massive stars are almost exclusively based on
single star models, even though a large fraction interact with a binary
companion. Using the MESA stellar evolution code, we predict the carbon
ejected in the winds and supernovae of single and binary-stripped stars at
solar metallicity. We find that binary-stripped stars are twice as efficient
at producing carbon (1.5–2.6 times, depending on choices on the slope of the
initial mass function and black hole formation). We confirm that this is
because the convective helium core recedes in stars that have lost their
hydrogen envelope, as noted previously. The shrinking of the core disconnects
the outermost carbon-rich layers created during the early phase of helium
burning from the more central burning regions. The same effect prevents carbon
destruction, even when the supernova shock wave passes. The yields are
sensitive to the treatment of mixing at convective boundaries, specifically
during carbon-shell burning (variations up to 40%) and improving upon this
should be a central priority for more reliable yield predictions. The yields
are robust (variations less than 0.5%) across our range of explosion
assumptions. Black hole formation assumptions are also important, implying
that the stellar graveyard now explored by gravitational-wave detections may
yield clues to better understand the cosmic carbon production. Our findings
also highlight the importance of accounting for binary-stripped stars in
chemical yield predictions and motivates further studies of other products of
binary interactions.

††software:  mesaPlot (Farmer, 2018), mesaSDK (Townsend, 2019),
ipython/jupyter (Pérez & Granger, 2007; Kluyver et al., 2016), matplotlib
(Hunter, 2007), NumPy (van der Walt et al., 2011), MESA(Paxton et al., 2011,
2013, 2015, 2018, 2019), and pyMesa (Farmer & Bauer, 2018).

##  1 Introduction

Understanding the Cosmic production of the elements that form the building
blocks of life is still one of the main quests for modern astronomy. Massive
stars are known to play a critical role in the synthesis of heavy elements
over cosmic time, but many questions remain open (e.g. Burbidge et al., 1957;
Cameron, 1959; Woosley & Weaver, 1995; Nomoto et al., 2006). Our current
understanding of the nucleosynthesis products of massive stars is still almost
exclusively based on single star progenitor models (Maeder, 1992; Woosley et
al., 1993a; Kobayashi et al., 2006; Curtis et al., 2019). However,
observational studies of young massive stars have indicated that massive stars
are nearly always born in multiple systems (Abt, 1983; Mason et al., 2009),
with at least one companion star nearby enough for binary interaction (Sana et
al., 2012; Moe & Di Stefano, 2017). Such interactions can dramatically change
the final fate of massive stars (Podsiadlowski et al., 1992; Wellstein &
Langer, 1999; Langer, 2012). How such interactions may change the chemical
yields of massive stars is still not fully clear (see, however, De Donder &
Vanbeveren 2004; Izzard 2004; Izzard et al. 2006; Woosley 2019). Models of
${}^{26}\rm{Al}$ nucleosynthesis find that binary stars can produce
significantly different yields from single stars (Braun & Langer, 1995;
Brinkman et al., 2019).

The source of carbon in the Universe is still uncertain (Bensby & Feltzing,
2006; Romano et al., 2017). Observations, and theoretical modelling, suggest
contributions from the winds of asymptotic giant branch (AGB) stars (Nissen et
al., 2014), massive star winds and their core collapse supernovae (Franchini
et al., 2020), and type Ia supernovae (Leung & Nomoto, 2020). The source of
carbon matters not just for the amount of carbon expected in the Universe, but
also for understanding when and where it is formed (Carigi et al., 2005;
Cescutti et al., 2009). Which can then be used to understand the star
formation history of a galaxy (Carilli & Walter, 2013; Romano et al., 2020).

Massive stars are able to eject carbon a few million years after formation
(Woosley et al., 1993b), while AGB winds and type Ia supernovae require much
longer timespans before releasing their carbon (Henry et al., 2000; Akerman et
al., 2004). The relative contribution from each source may vary over time as
the metallicity, and thus wind mass loss, increases (Dray & Tout, 2003; Lau et
al., 2020).

Carbon plays a crucial role in the interstellar medium (ISM) through its
complex chemistry and its ability to form a wide range of carbon-rich
molecules (Herbst & van Dishoeck, 2009) and carbonaceous dust (Li & Draine,
2001; Weingartner & Draine, 2001). Atomic carbon plays key roles in heating
and cooling interstellar gas (Wolfire et al., 1995) as well as in tracing the
properties of the ISM (Wolfire et al., 2003). The presence of CO is an
important observational tracer of molecular gas (Frerking et al., 1982;
Solomon et al., 1987). Dust formation from supernovae is also governed by the
presence of carbon (Bevan et al., 2017; Sarangi et al., 2018; Lau et al.,
2020; Brooker et al., 2021). Thus understanding the formation of carbon and
its distribution is key to understanding the ISM (Burton & Gordon, 1978;
Gullberg et al., 2018).

Here we study the effect of binary evolution on the carbon yields ejected by
massive stars (Langer, 1991). Previous work by Laplace et al. (2020, 2021)
explored the evolution of binary-stripped stars up to core collapse, and
showed how the mass loss during a binaries evolution alters the final
structure of a star. These structural differences in binary-stripped stars, as
compared to single stars, is expected to lead to differences in the final
supernova and the yields (Woosley, 2019; Schneider et al., 2021).

We take this work further by exploring the nucleosynthetic yield of carbon
($\mathrm{{}^{12}C}$) before and after core collapse and over a larger range
of initial masses. We consider here the fate of binaries that are stripped by
their companion in case B mass transfer, i.e they lose mass during their
evolution after core hydrogen depletion but before core helium ignition. An
analysis of all nucleosynthetic yields is deferred to later work.

Our paper is structured as follows, in Section 2 we describe our method for
following the evolution of single and binary stars, as well as their supernova
explosions. In Section 3 we compute the carbon yields for single and binary
stars. We discuss the uncertainties in our pre-supernova evolution and
supernova explosions in Section 4. In Section 5 discusses the initial mass
function (IMF) weighted yields. Finally, we discuss our results in Section 6
and conclude in Section 7.

##  2 Method

###  2.1 Pre-supernova evolution

We use the MESA stellar evolution code (version 12115, Paxton et al., 2011,
2013, 2015, 2018, 2019) to evolve massive single and binary stars from the
zero-age main sequence to core collapse. Our single stars and the primary
(initially most massive star) in the binary have initial masses between
$\rm{M_{init}}=11$ – $45\,\rm{M}_{\odot}{}$. For binary stars we set the
initial period to be between 38–300 days. This period range ensures that all
binary stars undergo case B mass transfer (Paczyński, 1967). We set the
secondary star’s mass such that the mass ratio $\rm{M_{2}/M_{1}}=0.8$. All
models are computed with a initial solar metallicity of Z=0.0142 and are non-
rotating. Single stars and binary stars are evolved with initial Y=0.2684
($\rm{Y=2Z+0.24}$ Pols et al., 1995a; Tout et al., 1996). Inlists with all
input parameters and models are made available online at
https://doi.org/10.5281/zenodo.4545836.

To evolve the systems we build upon the method in Laplace et al. (2020, 2021).
We follow in detail the structure of the primary star and the period evolution
of the system during Roche lobe overflow. We take the secondary star in the
binary to be a point mass and do not follow its evolution. Mass transfer is
assumed to occur conservatively such that no mass is lost from the system.
With the initial periods chosen we do not expect further Roche lobe overflow
(RLOF) to occur during the systems lifetime (Laplace et al., 2020). After core
helium burning ceases we evolve only the initial primary star of the binary.

Wind driven mass loss follows the prescriptions of Vink et al. (2001) for
stars with $\rm{T_{eff}}>10^{4}K$ and surface hydrogen mass fraction
$\rm{X_{H}}>0.4$, Nugis & Lamers (2000) for $\rm{T_{eff}}>10^{4}K$ and surface
$\rm{X_{H}}<0.4$, and de Jager et al. (1988) at all other times, with wind-
scaling factors of 1.0. In binary systems we define all of the mass lost when
the radius of the primary is greater than its Roche lobe radius to be RLOF,
even though mass loss via winds will still occur. This is a reasonable
assumption as the mass loss via winds during Roche lobe overflow is small due
to the short timescale over which RLOF occurs.

Convective overshoot is calibrated to that of Brott et al. (2011), with a step
overshoot value of $f=0.385$ and $f_{0}=0.05$. In MESA overshoot starts inside
a convection zone at a distance of $f_{0}$ (in pressure scale heights), and
extends from this point a distance of $f$ (in pressure scale heights).
Therefore, overshoot will extend a distance $f-f_{0}$ from the edge of the
convective boundary. We also apply the same amount of overshooting above the
metal burning zones during the late stage evolution of the models. We add a
small amount of overshoot ($f=0.05$, $f_{0}=0.01$) below metal burning shells
to improve numerical stability. We use MLT++ for all models to improve the
numerical stability of the low density envelopes (Paxton et al., 2013). We
include semiconvection with a mixing efficiency of $\alpha_{\rm{semi}}=1.0$
and we do not include thermohaline mixing. Additional physics choices are
specified in Appendix A.

We evolve our stars with MESA’s approx21.net which contains 21 isotopes,
following the alpha-chain up to iron. Farmer et al. (2016) has shown the need
to use larger nuclear networks when evolving stars to core-collapse, to
compute the core structure accurately. However, we show in Appendix B.1 that
models computed with approx21.net predict similar $\mathrm{{}^{12}C}$ yields
to models using such a larger network (mesa_128.net).

We define the helium core mass of the star as the first point in time (and
space) when the helium mass fraction, $\rm{X_{He}}>0.1$ and the hydrogen mass
fraction (at the same mass coordinate) is $\rm{X_{H}}<0.01$. Core oxygen
depletion is defined when the oxygen mass fraction at the center of the star
drops below $\rm{X_{O}}<10^{-4}$. Finally, we define core collapse to occur
when the inner regions of the star infalls at $300{\mathrm{km\ s^{-1}}}$.

###  2.2 Core-collapse supernovae

To model the supernova explosion, its shock, and the resulting nucleosynthesis
we place a “thermal bomb” at the center of our model (Aufderheide et al.,
1991; Sawada & Maeda, 2019). First we excise a portion of the star’s core, the
material that will form a compact object, by placing the inner boundary of our
model at the point where the entropy per baryon $S=4$ (Brown & Woosley, 2013).
We then inject energy into the base of the material outside this boundary over
a mass range of $0.01\,\rm{M}_{\odot}$, over 0.0045 seconds. We inject
sufficient energy to bring the total energy of the star (the sum of the
kinetic plus thermal energy minus the gravitational binding energy) to
$10^{51}\rm{erg/s}$. These values specify our default model assumptions.

This injection of energy then generates a hydrodynamic shock which travels
from the inner boundary of the star to the surface. As it passes through the
star it shock heats material and begins nuclear burning. This nucleosynthesis
is computed with MESA’s mesa_128.net which contains 128 isotopes up to
$\mathrm{{}^{64}Zn}$. In Section 4.2 we discuss how our choice of explosion
parameters affects the resulting nucleosynthesis. Also in Section 4.2 we
discuss our choice of temporal and spatial resolution during the explosion.
The star is then evolved until the shock reaches a location
$0.1\,\rm{M}_{\odot}$ below the surface, by which time the shock has cooled to
the point of no further nucleosynthesis occurring, except for beta decays. We
do not add by hand any $\mathrm{{}^{56}Ni}$ to the stars.

During the shock propagation through the star we track the energy change due
to photo-disintegrations and nuclear burning. At shock breakout our models
will have a different final energy as compared to the amount of energy we
injected during the explosion. The total energy is 1.05–1.20$\times 10^{51}$
$\rm{erg/s}$ at shock breakout while the kinetic energy at shock breakout is
between 0.5–1.3$\times 10^{51}$ $\rm{erg/s}$.

We define the yield of an isotope as (Karakas & Lugaro, 2016):

$\rm{Yield}=\sum_{T}\Delta M_{T}\times\left(X_{j}-X_{j,int}\right)$ (1)

Where $\Delta\rm{M_{T}}$ is the mass lost over the time interval $\rm{T}$,
$\rm{X_{j}}$ is the surface mass fraction of isotope $\rm{j}$, and
$\rm{X_{j,int}}$ is the initial mass fraction of isotope $\rm{j}$. With this
definition negative yields will occur in cases of net destruction of an
isotope. In this case, the mass fraction in the ejected material will be lower
than in the initial composition of the ejected material. We use the solar
composition of Grevesse & Sauval (1998) which sets
$\rm{X_{C_{12},int}=0.00244}$.

Figure 1: Top panel: The total $\mathrm{{}^{12}C}$ yield from all ejection
sources, Left Panel: The $\mathrm{{}^{12}C}$ yield ejected during wind mass
loss, Right Panel: The $\mathrm{{}^{12}C}$ yield ejected during core collapse.
The core-collapse yields assume all stars eject their envelope. In all panels
red regions denote binary models while blue regions denote single star models.
Open circles mark models which show anomalous carbon-burning behaviour, see
section 3.3. Dashed lines in panels (a) and (c) denote extrapolations over the
anomalous carbon-burning behaviour and models which do not reach core
collapse. The black arrows show the approximate location where each type of
mass loss dominates the $\mathrm{{}^{12}C}$ yield, taking into account
reasonable assumptions for which stars eject their envelopes (Sukhbold et al.,
2016; Zapartas et al., 2021).

##  3 Total carbon

Figure 2: The composition profile of a 15$\,\rm{M}_{\odot}$ binary-stripped
star at the start and end of the shock propagation phase of core collapse. The
colors denote the mass fraction of $\mathrm{{}^{12}C}$. The animated version
shows the shock propagation and resulting nucleosynthesis with time.

Figure 1(a) shows the total $\mathrm{{}^{12}C}$ yields from our single and
binary stars, from all sources of mass loss. Figure 1(b) shows the
$\mathrm{{}^{12}C}$ yield from winds, while Figure 1(c) shows the
$\mathrm{{}^{12}C}$ yield from core-collapse ejecta. In appendix C we include
table 2 that breaks down the total carbon yield by its source.

The $\mathrm{{}^{12}C}$ yield from the mass loss during RLOF is negative and
small $\approx-0.01\,\rm{M}_{\odot}$ and approximately independent of the
initial mass of the primary star. This is due to the envelope being
unprocessed and its carbon content reflects this initial abundance. The deeper
layers of the envelope have been processed by CN and CNO-cycling and show
depleted carbon abundances (Maeder, 1983). The most massive stars show a
slight decrease in $\mathrm{{}^{12}C}$ ejected during RLOF as a larger
fraction of the envelope is processed by CN(O)-cycling.

Figure 1(b) shows the $\mathrm{{}^{12}C}$ yields due to wind mass loss. The
mass loss due to winds (for all stars) can be broken into two groups, for
stars with $\rm{M_{init}}\lessapprox 35\,\rm{M}_{\odot}$ their winds are not
$\mathrm{{}^{12}C}$ enriched as compared to their initial composition, and
thus not visible in Figure 1(b). Stars with $\rm{M_{init}}\gtrapprox
35\,\rm{M}_{\odot}$ have $\mathrm{{}^{12}C}$ enriched winds. This is due to
both the single and binary stars becoming fully stripped, removing both the
hydrogen and helium layers of the star. In this initial mass range carbon rich
material has been mixed out of the core and is then ejected. This occurs at
$\rm{M_{init}}=37\,\rm{M}_{\odot}$ for single stars and
$\rm{M_{init}}=35\,\rm{M}_{\odot}$ for the binaries. The transition occurs at
a lower initial mass for the binaries due to RLOF removing some of the
envelope. Lower mass objects which are not fully stripped do not eject more
$\mathrm{{}^{12}C}$ in their stellar winds than they started with.

Figure 1(c) shows the $\mathrm{{}^{12}C}$ yields for core-collapse ejecta,
assuming that all stars eject their envelopes. The $\mathrm{{}^{12}C}$ yield
is relatively flat as a function of initial mass for stars with
$\rm{M_{init}}\lessapprox 27\,\rm{M}_{\odot}$, at $\approx
0.2\,\rm{M}_{\odot}$. Above this initial mass the yield rapidly increases up
to $\approx 1.25\,\rm{M}_{\odot}$. This transition between low and high carbon
yields occurs between $27\lessapprox\rm{M_{init}}/\,\rm{M}_{\odot}\lessapprox
35$. The increased $\mathrm{{}^{12}C}$ yields are due to the wind mass loss
removing the hydrogen envelope but not all of the helium envelope from the
stars. Thus there is enough mass loss to alter the core structure of the star,
leaving behind a carbon layer but not enough wind mass loss to then expose
that layer (see Section 3.2).

###  3.1 Shock nucleosynthesis

Figure 2 shows the distribution of $\mathrm{{}^{12}C}$ inside a binary-
stripped star with $\rm{M_{init}}=15\,\rm{M}_{\odot}$ at the start and end of
the shock propagation phase during core collapse. The online animation shows
the propagation of the shock and the explosive nucleosynthesis this generates.

At the start of the core-collapse phase, we can see that the
$\mathrm{{}^{12}C}$ is concentrated in the helium shell, with only small
amount of material near the core. As the shock propagates it photo-
disintegrates some of the $\mathrm{{}^{12}C}$ near the inner boundary, and
cools from a peak temperature of $\approx 10^{10}\rm{K}$, to $\approx 3\times
10^{8}\rm{K}$ when it reaches the carbon-rich helium shell, and to $\approx
10^{7}\rm{K}$ near shock breakout. By the time the shock reaches the carbon-
rich layers it has cooled sufficiently that it can no longer burn those
layers. The difference in pre- and post- supernova nucleosynthesis is $\approx
1\%$ for $\mathrm{{}^{12}C}$. The main effect on the $\mathrm{{}^{12}C}$ yield
is to smear the $\mathrm{{}^{12}C}$ distribution over a slightly larger range
of mass coordinates. Thus for future studies it is not necessary to model the
supernova explosion itself to predict $\mathrm{{}^{12}C}$ yields. This was
already found for single-star progenitors (e.g Thielemann et al., 1996; Young
& Fryer, 2007), and we now confirm it for self-consistent binary-stripped
structures.

Figure 3: Kippenhahn diagrams for the inner regions of single and binary-
stripped stars with $\rm{M_{init}}=19$, $37$, and $45\,\rm{M}_{\odot}$, during
core hydrogen and core helium burning. The left column shows single star
models, while the right column shows binary-stripped models. The x-axis shows
the time until core-collapse. Colors show the mass fraction of
$\mathrm{{}^{12}C}$ at each mass co-ordinate. Hatching shows mixing regions
due to convection and overshoot. The red horizontal line shows the mass
coordinate for what will become the CO core at the end of core helium burning.

###  3.2 Differences in core structure

To understand the differences between single stars and stripped binaries in
Figure 1 it is instructive to examine the evolution of two stars (a single and
a binary-stripped) of the same initial mass. Figure 3 shows the time evolution
of a single and a binary-stripped star with $\rm{M_{init}}=19,\ 37,\rm{and}\
45\,\rm{M}_{\odot}$. Considering the $\rm{M_{init}}=19\,\rm{M}_{\odot}$ case,
both stars start on the main sequence and lose only a very small amount of
mass before they evolve into Hertzsprung gap stars. At this point the donor
star in the binary exceeds its Roche lobe radius and begins losing mass, via
RLOF. In Figure 3b this occurs at
$\log_{10}\;\left(\rm{\tau_{cc}-\tau}\right)/\rm{yr}\approx 5.9$. Most of the
binary star’s envelope is lost at this point, bringing its mass down to
$\approx 8.0\,\rm{M}_{\odot}$, but the RLOF does not completely remove the
hydrogen envelope (Götberg et al., 2017; Yoon et al., 2017).

However, the mass that was lost was only comprised of the stellar envelope
(light blue region), which has a composition similar to the star’s birth
composition. Thus the mass loss from the binary during RLOF does not enrich
the Universe in $\mathrm{{}^{12}C}$. In fact, the material is slightly carbon
poor, due to some of the material having been CNO processed before the outer
edge of the convective hydrogen-burning core receded (dark blue region).

As the binary loses mass from its outer layers the core structure of the
binary is altered. The binary-stripped star forms a smaller helium core, and
the edge of the convective core recedes during core helium burning (Langer,
1989; Woosley et al., 1993b), while in the $\rm{M_{init}}=19\,\rm{M}_{\odot}$
single star the mass in the convective helium-burning core stays constant.
This receding convective core leaves behind $\mathrm{{}^{12}C}$ that was
produced by the $3\alpha$ process but had not yet been converted into oxygen
in ${}^{12}\rm{C}\left(\alpha,\gamma\right)^{16}\\!\rm{O}$ (Langer, 1991).
This left-over $\mathrm{{}^{12}C}$ is outside of what will become the CO core
of the star (denoted by the red lines).

Any $\mathrm{{}^{12}C}$ that is produced, and stays, inside the helium core
will either be burnt into oxygen at the end of core helium burning, destroyed
during $\mathrm{{}^{12}C}{}+\mathrm{{}^{12}C}{}$ burning, or will be accreted
onto the compact object during core-collapse. Thus only the
$\mathrm{{}^{12}C}$ that is mixed out of the core has a chance to survive
until core-collapse. Some of this outwardly mixed $\mathrm{{}^{12}C}$ will
not, however, survive until core collapse as it will be converted into oxygen
during helium shell burning, or be mixed back into the core during carbon
burning via the carbon convection zone intersecting the helium shell (see
Section 3.3).

Figure 3(c,d) shows a single and binary-stripped star with
$\rm{M_{init}}=35\,\rm{M}_{\odot}$. Here the wind mass loss is sufficient to
remove the remaining helium layers above the core. This exposes
$\mathrm{{}^{12}C}$ rich material, which is then ejected in a wind, similar to
the process in the binary-stripped star with
$\rm{M_{init}}=19\,\rm{M}_{\odot}$. Thus the single star has a net positive
$\mathrm{{}^{12}C}$ yield for its winds. The mass loss due to winds is also
now strong enough in the single-star case to force the helium core to recede.

Figure 3(e,f) shows stars with $\rm{M_{init}}=45\,\rm{M}_{\odot}$. At this
initial mass the winds of the single star are sufficient both to cause the
helium core to recede and to expose the $\mathrm{{}^{12}C}$ rich layers to be
ejected in a wind. Though this occurs at a slightly later time than in the
binary star case (as the binary-stripped star lost some mass in RLOF), thus
there is less time for the star to eject this carbon rich material leading to
a lower total carbon wind yield for the single star.

###  3.3 Carbon shell burning

Figure 4: Kippenhahn diagram for a 23$\,\rm{M}_{\odot}$ single star during
carbon and oxygen burning. Colors show the mass fraction of
$\mathrm{{}^{12}C}$ at each mass co-ordinate. Hatching shows mixing regions
due to convection and overshoot.

In Figure 1 there is significant noise in the $\mathrm{{}^{12}C}$ yields from
core-collapse ejecta, especially at the higher masses. There are also several
models marked with open symbols in the right panel of Figure 1. These are
$\rm{M_{init}}=24$, $28$, and $33\,\rm{M}_{\odot}$ for the single stars and
$\rm{M_{init}}=31\,\rm{M}_{\odot}$ for the binary-stripped stars. These
variations occur due to changes in the behaviour of carbon burning shells
during the star’s carbon and oxygen burning phases.

During the carbon shell burning phase, carbon initially ignites either at the
center of the star or off-center and burns outwards (Arnett & Truran, 1969;
Sukhbold & Adams, 2020). Figure 4 show the ignition of carbon in a
23$\,\rm{M}_{\odot}$ single star. This star ignites $\mathrm{{}^{12}C}$ at the
center at $\log_{10}\;\left(\rm{\tau_{cc}-\tau}\right)/\rm{yr}\approx 1.2$,
radiatively, this burning then moves outwards before it begins driving an off-
center convection zone. This initial burning phase stops, before an additional
$\mathrm{{}^{12}C}$ burning zone ignites off-center, at
$\log_{10}\;\left(\rm{\tau_{cc}-\tau}\right)/\rm{yr}\approx 0.0$ at the same
time as the core ignites $\mathrm{{}^{20}Ne}$. It is this burning zone that
causes the variability in the $\mathrm{{}^{12}C}$ yields.

As this convection zone extends outwards, it mixes $\mathrm{{}^{12}C}$ from
the outer layers of the core inwards where the $\mathrm{{}^{12}C}$ is then
burnt. Thus the maximal extent of this zone and the amount of time that it has
to mix $\mathrm{{}^{12}C}$ downwards sets the final $\mathrm{{}^{12}C}$
yields. Figure 4 shows that in this model a small pocket of
$\mathrm{{}^{12}C}$ survives between the outer edge of the carbon convection
zone and the lower edge of the helium burning shell (Laplace et al., 2021).
This pocket is a mix of $\mathrm{{}^{12}C}$ and $\mathrm{{}^{4}He}$ (left over
when the helium core receded at the end of core helium burning). At core
collapse, we find that for most models this pocket will contribute $\approx
50\%$ of the final $\mathrm{{}^{12}C}$ yield, with the remaining
$\mathrm{{}^{12}C}$ yield coming from $\mathrm{{}^{12}C}$ produced in the
helium shell.

We expect that differences in the treatment of mixing boundaries during
carbon-shell burning are also important for understanding the differences
between the results obtained with different codes, see Appendix B.2.

##  4 Physics variations

###  4.1 Sensitivity to physics choices

The sensitivity to the size and timing of the carbon-burning shells may
suggest that our models are under resolved, or that our choice of convective
overshoot above the carbon shell is significantly expanding the size of the
convection zone. To test this we ran two grids of additional models. Firstly
we ran models for with a $23\,\rm{M}_{\odot}$ single star model varying our
resolution controls. Next we randomly varied our choice of overshoot controls
($f$ and $f_{0}$) above the carbon shell for the same $23\,\rm{M}_{\odot}$
single star model (with our default resolution controls). We varied $f$ and
$f_{0}$ between 0.0 and 0.05, with $f_{0}<f$. For both sets of models we keep
our default model assumptions until the end of core helium burning, where we
then change our model assumptions and evolve the models until core oxygen
depletion.

We randomly picked values for each of the following spatial resolution
controls, mesh_delta_coeff (0.5–1.0), mesh_delta_coeff_for_highT (0.5–1.0),
and max_dq ($10^{-4}$–$10^{-3}$). For our temporal controls we varied
dX_nuc_drop_limit ($10^{-4.5}$–$10^{-3.0}$) and varcontrol_target
($10^{-4.5}$–$10^{-3.0}$). With our choices we increase the spatial resolution
by up to a factor of 5 and increase (decrease) the temporal resolution by a
factor of 3 (10). This leads to the highest resolution models having $\approx
100,000$ time steps (for just carbon and oxygen burning) and $\approx 15,000$
mesh points. Overall, we are primarily sensitive to the choice of
dX_nuc_drop_limit, which limits the timestep based on the rate of change of
the most abundant isotopes, primarily variations in the $\mathrm{{}^{12}C}$
and $\mathrm{{}^{16}O}$ abundances at the center of the star.

(a) Temporal resolution

(b) Overshoot

Figure 5: The relative change in the final total mass of $\mathrm{{}^{12}C}$
in the star measured at the point of core oxygen depletion, as a function of
the temporal resolution and strength of convective overshoot for a
$23\,\rm{M}_{\odot}$ single star. In panel (a) we vary a range of MESA’s
temporal and spatial resolution controls. In panel (b) we vary both of the
overshoot parameters $f$ and $f_{0}$. The x-axis in panel (a) is the average
timestep relative to our default model, where the resolution increases towards
the left. The x-axis in panel (b) is the physical extent of convective
overshoot beyond the convective boundary above carbon burning shells. In both
panels the orange star denotes our default model assumptions, and the black
plus symbol denotes a model which was evolved with mesa_128.net from the ZAMS
(with otherwise default assumptions).

Figure 5 shows the relative change in the total mass of $\mathrm{{}^{12}C}$ in
our model measured at the end of the core oxygen depletion for both variations
in the resolution and the amount of convective overshoot above the carbon
burning shell. The black plus symbols in Figure 5 show a comparison with a
model ran with mesa_128.net from the ZAMS, which has only a $\approx 2\%$
difference in the total mass of $\mathrm{{}^{12}C}$ (measured at core oxygen
depletion).

Figure 5(a) shows the relative change in the total mass of $\mathrm{{}^{12}C}$
compared to the relative change in the average time step taken. As the time
resolution increases the total mass of $\mathrm{{}^{12}C}$ at the end of core
oxygen burning decreases, relative to our default choices. However, even with
at our highest resolutions, there is still a large spread possible in the
final mass of $\mathrm{{}^{12}C}$ ($\approx 40\%$). We do not plot the spatial
resolution variations as they show no correlation with the carbon yields.

Figure 5(b) shows the relative change in the $\mathrm{{}^{12}C}$ yields as a
function of the amount of overshoot beyond the top of the carbon burning
shells convective boundary (with all other overshoot regions keeping their
same values). There is a slight trend for both the final total mass of
$\mathrm{{}^{12}C}$ to decrease, and the spread in $\mathrm{{}^{12}C}$ values
to decrease, as the amount of overshoot increases. However, it is possible to
achieve the same final total mass of $\mathrm{{}^{12}C}$ as our default model
with the full range of overshoot values considered here. While this shows we
may not be sensitive to the overall amount of overshoot (above the carbon
burning shells), there is considerable scatter in the final total mass of
$\mathrm{{}^{12}C}$ with the the same choice of the amount of overshoot. Thus
MESA users should consider carefully their physical choices for the amount of
overshoot ($f-f_{0}$) and the individual numerical choices ($f$ and $f_{0}$)
needed to achieve this value.

More work is needed to understand convection during this burning phase
(Cristini et al., 2017), as well as the role of convective overshoot that can
enhance this effect by extending the region over which the carbon convection
zone can mix. We may be seeing a similar effect as shown in Paxton et al.
(2018) and Paxton et al. (2019) with the improper placement of the convective
boundary. We tested with both the predictive mixing (Paxton et al., 2018) and
convective pre-mixing schemes (Paxton et al., 2019). Both schemes still show
scatter in how large the carbon convective shell grows, and for how long it is
able to mix $\mathrm{{}^{12}C}$ into the core. We note however that both
predictive mixing and convective pre-mixing assume that the mixing timescale
is shorter than the computational timestep, which breaks down during carbon
burning.

###  4.2 Sensitivity to the explosion properties

(a) Injected energy

(b) Mass cut

(c) Injection time

(d) Injection mass

(e) Spatial resolution

(f) Temporal resolution

Figure 6: The relative difference in the amount of $\mathrm{{}^{12}C}$ ejected
by supernovae of our binary-stripped star model with
$\rm{M_{init}}=16\,\rm{M}_{\odot}$, compared to our default ($\rm{def}$)
model. Panel (a) shows the variation with respect to the injected energy,
Panel (b) shows the variation with respect to the mass cut, Panel (c) shows
variations with respect to the time the energy is injected over, Panel (d)
shows variations with respect to the range of masses that the energy is
injected over, Panel (e) shows the variation with respect to changes in the
spatial resolution, and Panel (f) shows the variation with respect to changes
in the temporal resolution. $\langle dm\rangle$ is the time-weighted average
cell mass in $\,\rm{M}_{\odot}$, $\langle dt\rangle$ is the average timestep
in seconds. The spatial and temporal resolution increases to the left. For
each row, we simultaneously Monte Carlo sampled both the parameters shown
while keeping the other explosion parameters set to our default assumptions.
The orange star denotes our default assumptions.

To explore the sensitivity of our results to the assumptions made we perform
two sets of tests on a $\rm{M_{init}}=16\,\rm{M}_{\odot}$ binary-stripped
star. We do not extensively test the pre-supernova physics variations as they
have previously been explored in Farmer et al. (2016); Renzo et al. (2017);
Laplace et al. (2021). Figure 6 shows the effect of varying the physics
assumptions made during the core-collapse supernovae and the effects of
varying the numerical assumptions made during the core-collapse supernova.
Each row of Figure 6 tests two physics/numerical assumptions at a time. The
vertical spread in the distribution of points at a fixed x-coordinate shows
how sensitive the $\mathrm{{}^{12}C}$ yield is to the _other_ parameter shown
in the same row. Thus a tight correlation indicates that the _other_ parameter
in the same row only minimally affects the predicted $\mathrm{{}^{12}C}$
yield. A horizontal line implies that the $\mathrm{{}^{12}C}$ yield is
insensitive to _that_ parameter. Negative
$\rm{\left(M_{c_{12}}-M_{c_{12},def}\right)/M_{c_{12},def}}$ values indicate
that more $\mathrm{{}^{12}C}$ was destroyed than in our default model, while
positive values indicate that less $\mathrm{{}^{12}C}$ was destroyed.

In Figure 6(a) and 6(b) we randomly sampled $\approx 200$ times both the
injected energy (between $0.5$ – $5\times 10^{51}\rm{erg/s}$) and the mass cut
(between $1.4$ – $2.0\,\rm{M}_{\odot}$). Figure 6(a) shows the effect of
varying the injected energy on the total mass of $\mathrm{{}^{12}C}$ ejected
during the supernovae. We can see a strong correlation between injected energy
and $\mathrm{{}^{12}C}$ ejected, increasing the energy decreases the
$\mathrm{{}^{12}C}$ mass ejected. However, the change is $<1\%$ in the total
mass of $\mathrm{{}^{12}C}$ and in the context of Figure 1 the change would
only be of the order of the size of the symbols.

Figure 6(b) shows the effect as we vary the mass cut. We find no correlation
with the mass of ejected $\mathrm{{}^{12}C}$ for this range of mass cuts.
These trends can be explained as follows: most of the $\mathrm{{}^{12}C}$
resides in the helium shell and only a small amount has been mixed down (or
produced during late stage burning) to near where the compact object will
form. Thus changing the mass cut has little effect, as the total mass of
$\mathrm{{}^{12}C}$ ejected is much greater than the change possible by moving
the inner boundary.

The injected energy has a correlated (but small) effect on the
$\mathrm{{}^{12}C}$ yield due to the shock processing the small amount of
$\mathrm{{}^{12}C}$ that exists near the inner boundary (see Figure 2).
However by the time the shock reaches the bulk of the $\mathrm{{}^{12}C}$ in
the helium shell the shock has cooled below $10^{8}$K, thus it can no longer
burn the $\mathrm{{}^{12}C}$ present.

In Figure 6(c) and 6(d) we randomly sampled $\approx 200$ times both the time
over which we inject energy during the explosion
($10^{-3}<\rm{T_{inj}}/s<10^{0}$) and the mass over which we inject energy
($10^{-2}<\rm{\Delta M_{inj}}/\,\rm{M}_{\odot}<5\times 10^{-1}$). When either
the $\rm{T_{inj}}$ is small ($\rm{T_{inj}}<0.1$s) or $\rm{\Delta M_{inj}}$ is
small ($\rm{M_{inj}}<0.1\,\rm{M}_{\odot}$) their effects on the
$\mathrm{{}^{12}C}$ yield is small. Only once $\rm{M_{inj}}$ exceeds
$0.1\,\rm{M}_{\odot}$ does it begin to dominate the $\mathrm{{}^{12}C}$ yield.
However, once $\rm{T_{inj}}>0.3$s it begins to dominate the
$\mathrm{{}^{12}C}$ yield instead. When either $\rm{T_{inj}}$ or $\rm{\Delta
M_{inj}}$ increases the generated shock is weaker as either its power is lower
(as it spreads over more time) or its deposited energy per unit mass is lower
(as its spread over more mass) . This weakening of the shock lowers the amount
of $\mathrm{{}^{12}C}$ destroyed in the initial burning phases of the
explosion (Sawada & Maeda, 2019). The variations here are comparable to the
uncertainty in the injected energy (though with an opposite sign).

For Figure 6(e) and 6(f) we randomly sampled $\approx 100$ times both the
temporal controls and the mesh controls, simultaneously, to probe the
sensitivity of our predictions to the numerical resolution of our models. We
varied a series of mesh controls (split_merge_amr_nz_baseline (500–8000),
split_merge_amr_MaxLong (1.1–2.5), and split_merge_amr_MaxShort (1.1–2.5)),
during the explosion, as seen in Figure 6(e). These controls force MESA to
distribute its mesh according to the radius of each zone (See section 4 of
Paxton et al. (2018)). We varied dt_div_min_dr_div_cs_limit, which sets the
timestep based on the sound-crossing timescale for spatial zones near the
shock front, between $0.1$ – $10$. While MESA is an implicit code and thus not
limited by the sound-crossing timescale, it provides a physical and convenient
timescale over which to probe the numerical sensitivity of our models.

As can be seen in Figure 6(e), increasing the mesh resolution by a factor of 5
shows changes smaller than those which result from the variations we tested in
the injected energy or choice of mass cut. Increasing (or decreasing) the
temporal resolution also shows changes in ejected $\mathrm{{}^{12}C}$ mass of
$\lessapprox 0.05\%$ (Figure 6(f)). Thus our numerical uncertainties during
the supernova are much smaller than our uncertainties due to either the
stellar models or the physical explosion parameters.

These results suggest our $\mathrm{{}^{12}C}$ estimates are therefore not
sensitive to the uncertain parameters assumed for the explosion, however other
isotopes which are formed deeper into the star are more affected by both the
injected energy and the mass cut (Young & Fryer, 2007; Suwa et al., 2019;
Sawada & Maeda, 2019). While our results show that the total amount of
$\mathrm{{}^{12}C}$ ejected is insensitive to the physical explosion
parameters considered here, the amount of $\mathrm{{}^{12}C}$ that is
observable (via the production of carbon-rich dust at late times) is sensitive
to the explosion physics (Brooker et al., 2021).

##  5 IMF weighted yields

Table 1: Ratio of the IMF weighted yields between an equal number of massive binary-stripped stars and single stars for different assumptions about the ejection of the envelope during core collapse and the IMF power-law $\alpha$, see also Equation 2. In case of BH formation we assume that the carbon-rich layers fall back onto the BH. BH formation assumption | Ratio of $\mathrm{{}^{12}C}$ yields (binary-stripped/single)
---|---
| $\alpha=-1.9$ | $\alpha=-2.3$ | $\alpha=-2.7$
All stars explode successfully | 1.51 | 1.44 | 1.36
BHs form from stars with initial masses $\rm{\rm{M_{init}}>22}\,\rm{M}_{\odot}$ | 1.79 | 1.55 | 1.35
BHs form from stars final core masses $\rm{M_{He,final}}>7\,\rm{M}_{\odot}$ | 2.59 | 2.22 | 1.91
Schneider et al. 2021 like formation | 2.57 | 2.21 | 1.92

The limit between stars which form neutron stars (and are assumed to eject
their envelopes) and those that form black holes (which may or may not eject
their envelopes) is uncertain (O’Connor & Ott, 2011; Brown & Woosley, 2013).
Whether a star ejects its envelope or not will strongly affect the final
yields from that star. To probe this uncertainty in whether a star ejects its
envelope during core collapse we test several different ejection assumptions.
Table 1 show the IMF weighted ratio of the total $\mathrm{{}^{12}C}$ yields
for binary-stripped massive stars and single massive stars at solar
metallicity. We define the ratio as:

$\rm{Ratio=\frac{\int_{11}^{45}\left(Y_{b,ccsn}f_{b}+Y_{b,winds}+Y_{b,rlof}\right)M_{init,b}^{\alpha}\,dM}{\int_{11}^{45}\left(Y_{s,ccsn}f_{s}+Y_{s,winds}\right)M_{init,s}^{\alpha}\,dM}}$
(2)

Where we assumed a Salpeter-like IMF with different values of $\alpha$
(Schneider et al., 2018), $\rm{M_{init,b}}$ and $\rm{M_{init,s}}$ are the
initial mass of the primary star in the binary and the mass of the single star
in solar masses, $\rm{Y_{b}}$ and $\rm{Y_{s}}$ are the yields of
$\mathrm{{}^{12}C}$ in solar masses for the respective stars and for each type
of mass loss, and $\rm{f_{b}}$ and $\rm{f_{s}}$ are filter functions that are
either 0 or 1 depending on whether we assume that the star ejects its
envelope. The integration limits are taken over the entire range of initial
masses considered in this work.

Firstly assuming that all stars eject their envelopes at core-collapse then we
find binary-stripped stars contribute $\approx 40\%$ more $\mathrm{{}^{12}C}$
compared to the same initial mass of single massive stars, assuming a standard
Salpeter IMF $\alpha=-2.3$. This is due to extra mass loss in binary systems,
leading both to a higher $\mathrm{{}^{12}C}$ yield in the winds and in the
final supernovae for $\rm{M_{init}}\approx 30$ – $40\,\rm{M}_{\odot}$ binary-
striped stars as compared to single stars.

As we do not expect all stars to eject their envelopes at core-collapse we can
also filter out systems which we do not expect to eject their envelope.
Sukhbold & Woosley (2014) found that single stars with
$\rm{M_{init}}\lessapprox 22\,\rm{M}_{\odot}$ had a low compactness (O’Connor
& Ott, 2011), suggesting they would be likely to successfully eject their
envelope, this would set the maximum helium core mass as
$\rm{M_{He}}\lessapprox 7\,\rm{M}_{\odot}$.

Assuming all stars eject their envelopes with
$\rm{M_{init}}<22\,\rm{M}_{\odot}$, then stripped primary stars in massive
binaries contribute even more to the $\mathrm{{}^{12}C}$ in the Universe than
single massive stars.

Next, taking the helium core mass at $\rm{M_{He,final}}<7\,\rm{M}_{\odot}$,
measured at core collapse, binary-stripped massive stars eject approximately
twice as much carbon than single massive stars (for $\alpha=-2.3$). This is
due to the extra mass loss binaries undergo, lowering the final helium core
masses. Thus for a given helium core mass at core-collapse binary-stripped
systems have a higher initial mass (Kippenhahn & Weigert, 1967; Habets, 1986),
expanding the range in initial masses over which we assume a successful
explosion occurs (Vartanyan et al., 2021). For our systems, a cut of
$\rm{M_{He,final}}<7\,\rm{M}_{\odot}$ is equivalent to a cut of
$\rm{M_{init}}<19\,\rm{M}_{\odot}$ for single stars and
$\rm{M_{init}}<22\,\rm{M}_{\odot}$ for binaries.

Finally, we use the limits for BH formation found in Table 1 of Schneider et
al. (2021). They found that NSs would form from two populations, firstly for
single stars with $\rm{M_{init}}\leq 21.5\,\rm{M}_{\odot}$ and for stars which
were stripped as case B binaries with $\rm{M_{init}}\leq
31.5\,\rm{M}_{\odot}$. Then a second population of NSs are formed from more
massive objects, for single stars in the range
$23.5\leq\rm{M_{init}}/\,\rm{M}_{\odot}\leq 34.0$ and for case B binaries with
$34.0\leq\rm{M_{init}}/\,\rm{M}_{\odot}\leq 67.5$. This second population
arises from changes in the core carbon and neon burning, leading to
differences in whether the burning is convective or radiative, which leads to
differences in the final core compactness and structure (see also Brown et
al., 1999, 2001).

As the IMF $\alpha$ increases, thus favouring the production of more massive
stars, the contribution that binary-stripped makes to the $\mathrm{{}^{12}C}$
yields increases. This is due the greater weight now given to the
$\mathrm{{}^{12}C}$ yields from wind mass loss (as only the most massive stars
in our grid contribute to wind mass loss yields), which is where the
difference between the binary-stripped and single star yields is greatest (See
Figure 1).

For all envelope assumptions the contribution from the wind mass loss is
greater for stripped binaries then for the single stars. This is due to both
the binary-stripped stars having $\mathrm{{}^{12}C}$ positive yields at lower
initial masses, and due to the higher $\mathrm{{}^{12}C}$ yield at the
equivalent initial masses due to the extra mass loss from RLOF altering the
core structure (Section 3.2). Table 3 in Appendix C shows how the
$\mathrm{{}^{12}C}$ yields from only core collapse vary as a function of the
envelope ejection assumptions. Adopting the prescription of Schneider et al.
(2021) does not lead to a significantly different ratio of carbon yields from
when using the simple helium core mass cut. Firstly, this is due to the ratio
of the yields being dominated by the wind yield rather than the core-collapse
yield, and secondly because increasing the initial mass range for successful
envelope ejections increases both binary-stripped and single star core-
collapse yields by approximately the same relative amount.

The actual contribution to carbon enrichment in the Universe also depends on
an additional scaling factor weighting the fraction of stars that are single
against those that are binary-stripped (Sana et al., 2012), as well as the
fraction of binary systems which do not self-strip. Note also that we are
comparing equivalent total initial masses for ensembles of massive single
stars and the one star in each binary that we model as stripped, not the full
initial stellar mass of the binaries. We have also not included the
$\mathrm{{}^{12}C}$ yields from the secondary stars in the binaries. If the
secondary gains significant mass by accretion, we anticipate that this would
further increase the relative efficiency of massive binary systems in
producing $\mathrm{{}^{12}C}$.

We can combine the result shown in Figure 1(b), indicating that the
$\mathrm{{}^{12}C}$ yields from winds of massive stars are not positive until
$\rm{M_{init}}\approx 35\,\rm{M}_{\odot}$, with the expectation that only
massive stars with $\rm{M_{init}}\lessapprox 22\,\rm{M}_{\odot}$ eject their
envelopes at core collapse. This combination suggests that stars with initial
masses $22\lessapprox\rm{M_{init}}/\,\rm{M}_{\odot}\lessapprox 35$ do not
contribute to the net production of $\mathrm{{}^{12}C}$ in the Universe. The
contribution of massive stars with $\rm{M_{init}}<22\,\rm{M}_{\odot}$ can be
well described by only their $\mathrm{{}^{12}C}$ supernovae yield, while stars
with $\rm{M_{init}}>35\,\rm{M}_{\odot}$ are well described by only their wind
yields.

##  6 Discussion

This work has only modelled solar-metallicity stars, has neglected the
potential importance of rotational mixing, and our binary-stripped models all
underwent early case B mass transfer. We have also tested only one set of
physics assumptions for stellar evolution, even though there are many
uncertainties in the evolution of massive stars. Here we discuss potential
limitations from those approximations.

Wind mass loss prescriptions for massive stars are uncertain (Renzo et al.,
2017). Sander & Vink (2020) found that we would expect weaker Wolf-Rayet
winds, at solar metallicity, for stars in the $\rm{M=}10$ –
$15\,\rm{M}_{\odot}$ range than for our choice of Nugis & Lamers (2000).
Lowering the wind mass-loss rate in this mass range would likely lead to lower
supernova $\mathrm{{}^{12}C}$ yields from binaries (Langer & Henkel, 1995;
Eldridge & Tout, 2004). As the stripped star would lose less mass during core
helium burning, the helium core would recede less, leaving behind less
$\mathrm{{}^{12}C}$ in the helium shell. This would decrease the significance
of binary-stripped stars on the total $\mathrm{{}^{12}C}$ yields. Single stars
would be unaffected as they do not self-strip until much higher masses, where
the models of Sander & Vink (2020) agree with those of Nugis & Lamers (2000)
(until the initial masses go above those considered here, where Sander & Vink
(2020) would predict higher mass loss rates). See also Dray et al. (2003) for
discussion of the effects of the choice of wind mass-loss prescription on
carbon yields from Wolf-Rayet stars.

Correlated with the wind mass loss is the modelled metallicity. Lower
metallicity models would show a similar effect as having a lower assumed wind
mass loss rate by changing the fraction of stars that become stripped
(Eldridge & Vink, 2006). Also at low metallicities, RLOF in binaries may not
fully strip the star due to changes in the envelope physics (Götberg et al.,
2017; Yoon et al., 2017; Laplace et al., 2020). Without the envelope becoming
fully stripped we would not see the helium core recede which would lower both
the wind and core-collapse $\mathrm{{}^{12}C}$ yields.

Our binary star models were each given an initial period such that they would
undergo early case B mass transfer, before the onset of core helium burning.
However, binaries may also interact earlier during the main sequence (case A),
or after core helium burning has started (case C). Mass lost during case A
will primarily lead to the star having a smaller core mass (Wellstein &
Langer, 1999; Yoon et al., 2010). If the mass loss ceases before the star
undergoes core helium burning, and the star retains its hydrogen envelope,
then the helium core will not recede. As the core did not recede it will not
leave behind a pocket of $\mathrm{{}^{12}C}$ above it. Thus the winds from the
binary-stripped star will not be enhanced in $\mathrm{{}^{12}C}$ yields
relative to single stars. Countering this, smaller cores would lead to more
successful explosions, ejecting the envelope, which acts to increase the total
$\mathrm{{}^{12}C}$ yields. More work is needed to understand which effect
dominates. Mass lost during case C will occur after core helium burning has
started. By this time the core mass of the star is set, thus we would be
unlikely to see an enrichment of the helium shell with $\mathrm{{}^{12}C}$
from the helium core.

Rotation plays a key role in the evolution of stars (Meynet & Maeder, 2000)
and the observed chemical composition of massive stars (Hunter et al., 2007,
2008, 2009). The extra chemical mixing rotation generates inside the star can
lead to larger core masses (Ekström et al., 2012; Murphy et al., 2021), as
well as mixing elements to the surface of the star (Meynet et al., 2006;
Maeder et al., 2014; Groh et al., 2019). Prantzos et al. (2018) found
enhancements in $\mathrm{{}^{12}C}$ of a approximately a factor of three,
between initial rotation rates of $0$ – $300{\mathrm{km\ s^{-1}}}$ for a
$20\,\rm{M}_{\odot}$ star. Hirschi et al. (2004, 2005) found that the
$\mathrm{{}^{12}C}$ yields increase by a factor $1.5$ – $2.5$ for stars
$\rm{M_{init}}<30\,\rm{M}_{\odot}$, and with initial rotation rates of $0$ –
$300{\mathrm{km\ s^{-1}}}$ due to the more massive core that the rotating
models produce.

Fields et al. (2018) explore the sensitivity of a 15$\,\rm{M}_{\odot}$ solar
metallicity star to variations when all reaction rates in the model were
varied within their measured uncertainties. They found that at core helium
depletion the central carbon fraction can vary $\pm 80\%$ compared to the
model which adopts median reaction rates, driven by variations in
${}^{12}\rm{C}\left(\alpha,\gamma\right)^{16}\\!\rm{O}$ and the
triple-$\alpha$ rates. Sukhbold & Adams (2020) showed how the type of carbon
burning, radiative or convective, is altered by changing the
${}^{12}\rm{C}\left(\alpha,\gamma\right)^{16}\\!\rm{O}$ rate, in 14 –
26$\,\rm{M}_{\odot}$ stars. The location of carbon ignition plays a key role
in determining how compact the core will become and thus how likely it is to
successfully explode (Brown et al., 1999, 2001), and thus whether it will
eject its envelope (Weaver & Woosley, 1993; Timmes et al., 1996). Sukhbold &
Adams (2020) found that variations in the $\approx 1\sigma$ uncertainty in the
${}^{12}\rm{C}\left(\alpha,\gamma\right)^{16}\\!\rm{O}$ can move the location
between convective and non-convection carbon ignition by $\pm
2\,\rm{M}_{\odot}$ in initial mass. That change in initial masses is similar
to the change in initial masses seen in Section 5 when assuming stars with a
fixed final helium core mass eject their envelopes.

##  7 Conclusion

Motivated by the high binary fraction inferred for young massive stars and the
importance in improving chemical yield predictions, we have started a
systematic investigation of the impact of binarity on the chemical yields of
massive stars. In this paper, intended as the first in a series, we focus on
$\mathrm{{}^{12}C}$ yields. We present a systematic comparison of the
differences between massive stars stripped by a binary companion and single
stars with the same initial mass at solar metallicity.

To achieve this, we modelled the evolution of from the onset of hydrogen
burning until core collapse. We then followed the nucleosynthesis through the
resulting supernova shock and compute the ejected $\mathrm{{}^{12}C}$ yield
for different mass loss processes.

Our results can be summarised as follows:

  1. 1.

We find that massive stars stripped in binaries during hydrogen-shell burning
are nearly twice as efficient at contributing to the cosmic
$\mathrm{{}^{12}C}$ production as a similar number of massive single stars (a
factor 1.4–2.6 depending on the assumptions for black hole formation and the
slope of the initial mass function, see Equation 2 and Table 1).

  2. 2.

We confirm that the difference in yields between binary-stripped and single
massive stars can be explained by considering the outer-most
$\mathrm{{}^{12}C}$-rich layers produced in the early phases of central helium
burning. In single stars these layers tend to get mixed into the growing
convective core, leading to destruction of $\mathrm{{}^{12}C}$ by alpha
captures and later carbon burning. In binary-stripped stars (and in single
stars stripped by stellar winds) the convective helium-burning core cannot
grow and may even retreat. The outer-most $\mathrm{{}^{12}C}$-rich layers
disconnect and form a pocket where the temperature never becomes high enough
for $\mathrm{{}^{12}C}$ destruction (Figure 3), not even when a supernova
shock wave passes through and ejects these layers (Figure 2) (cf. Langer,
1991; Woosley, 2019; Laplace et al., 2021).

  3. 3.

Stellar winds also eject $\mathrm{{}^{12}C}$ once the carbon-rich layers are
exposed to the surface. This only happens for the most massive stars in our
grid (initial masses $>36{\mathrm{M}_{\odot}}$ and $>38{\mathrm{M}_{\odot}}$
for our solar-metallicity binary-stripped and single stars respectively).
Comparing stars of the same initial mass, we find that the wind yields for
binary-stripped stars are higher because the carbon-rich layers are more
massive and appear earlier at the surface (see Figure 1 and 1b).

  4. 4.

Mass loss from binary systems during conservative mass transfer does not
contribute to $\mathrm{{}^{12}C}$ yields, in our models. These layers are
either pristine or contain CNO-processed material which is $\mathrm{{}^{12}C}$
poor (c.f. de Mink et al., 2009).

  5. 5.

Our yield predictions are remarkably robust with respect to choices in the
treatment of the explosion, such as the explosion energy, how the energy is
inserted, and the mass that promptly forms a compact object. These lead to
variations smaller than $\approx 0.5\%$ (See Figure 6).

  6. 6.

We show that the $\mathrm{{}^{12}C}$ yield predictions of both single and
binary-stripped stars are sensitive to the treatment of overshooting,
specifically during carbon-shell burning. We find variations up to $\approx
40\%$ in the final carbon yields for individual models (see Figure 5),
depending on the size and lifetime of the shells. We identify this as a
primary source of uncertainty and cause for noise in our predictions of the
supernova yields.

  7. 7.

We note that the variations between the predictions of $\mathrm{{}^{12}C}$
yields for massive stars presented in other recent studies (Pignatari et al.,
2016; Limongi & Chieffi, 2018; Griffith et al., 2021) is large and of a
similar order as the differences we find between binary-stripped and single
stars.

We conclude that the yields of massive binary-stripped stars are
systematically different from those of massive single stars. The effects of
binarity should therefore not be ignored if we want to obtain more reliable
yield predictions for carbon and better understand the relative contribution
of massive stars with respect to the other proposed sources of carbon, namely
Asymptotic Giant Branch stars and type Ia supernova.

However, the main priority would be to better understand how to treat boundary
mixing, especially during the carbon shell burning phase, as this seems to be
the primary uncertainty in both single and binary star models. Further
aspects, which deserve further investigation, include the effects of
metallicity and how these interplay with uncertainties in the stellar wind
mass loss, and the still unknown conditions for which stars lead to successful
supernova explosions. Advancing in these areas would require efforts on both
the theory and observation side.

Finally, we emphasize that the binary-stripped stars studied here represent
only one of the many possible final fates for a star born in the close
vicinity of a companion. Binary interaction can affect massive stars in many
other ways. This study is thus only a first step towards answering the bigger
question of how binarity impacts the chemical yields of stars.

##  Appendix A Other physics choices

The EOS in MESA is blended from several sources; OPAL (Rogers & Nayfonov,
2002), SCVH (Saumon et al., 1995), PTEH (Pols et al., 1995b), HELM (Timmes &
Swesty, 2000), and PC (Potekhin & Chabrier, 2010). Opacities are primarily
from OPAL (Iglesias & Rogers, 1993, 1996) with additional data from Buchler &
Yueh (1976); Ferguson et al. (2005); Cassisi et al. (2007)

Nuclear reaction rates are a combination of rates from NACRE and JINA’s
REACLIB (Angulo et al., 1999; Cyburt et al., 2010), with additional weak
reaction rates from Fuller et al. (1985); Oda et al. (1994); Langanke &
Martínez-Pinedo (2000). Nuclear screening is computed with the prescription of
Chugunov et al. (2007). Thermal neutrino loss rates are from Itoh et al.
(1996). We compute the Roche lobe radii in binary systems using the fit from
Eggleton (1983). The mass transfer rate in our Roche lobe overflowing binary
systems is computed from the prescription of Ritter (1988).

##  Appendix B Comparison to other works

###  B.1 Laplace et al. 2021

Figure 7: Total $\mathrm{{}^{12}C}$ mass ejected during supernovae, comparing
models from this work (dashed) with those from Laplace et al. (2021) (solid).
In blue are plotted single-star models, while red represents binary-stripped
models. The models of Laplace et al. (2021) were evolved with mesa_128.net
after core oxygen depletion, and then exploded with the same method as in
Section 2.

In Figure 7 we compare the $\mathrm{{}^{12}C}$ mass ejected in the supernova
calculations presented in Figure 1 to the $\mathrm{{}^{12}C}$ ejected when
using progenitor models computed in Laplace et al. (2021). We took the models
from Laplace et al. (2021) and exploded them with the same method as in
Section 2. There are two primary differences between these sets of models;
Laplace et al. (2021) was computed with MESA r10398, while this work uses MESA
r11215; and Laplace et al. (2021) uses mesa_128.net after core oxygen
depletion, while this work continues to use approx21.net after core oxygen
depletion. For both sets of models we explode them with MESA r11215.

We can see that for both single and binary-stripped star models the amount of
$\mathrm{{}^{12}C}$ ejected in the supernovae is similar between both sets.
The largest disagreement occurs for the binary-stripped stars at
$\rm{M_{init}}=21\,\rm{M}_{\odot}$. This is due to the carbon shell destroying
more $\mathrm{{}^{12}C}$ in the pre-supernova phase of the evolution. As this
carbon shell acts before the switch to the larger nuclear network this
indicates the differences are due to changes in MESA between r10398 and
r11215, rather than to differences in the nuclear network.

This indicates that it is reasonable, when considering the $\mathrm{{}^{12}C}$
yields, to use a smaller but more computational efficient, nuclear network
during the pre-supernova evolution. This is due to two reasons, first most of
the $\mathrm{{}^{12}C}$ that will be ejected is set by the evolution up to the
end of core-carbon burning, which does not require large nuclear networks to
compute. Secondly, the $\mathrm{{}^{12}C}$ that is ejected in the supernova
resides in/near the helium shell. Thus the total $\mathrm{{}^{12}C}$ yield is
insensitive to the structure of the inner core, which is sensitive to the
choice of nuclear network (Farmer et al., 2016). This however may not be true
for other isotopes, or for inferring neutron star masses which depend
sensitively on pre-supernova structures.

###  B.2 Other core-collapse works

Figure 8: Comparison of nucleosynthetic $\mathrm{{}^{12}C}$ yields for solar
metallicity non-rotating single stars. In blue this work, magneta line
(Pignatari et al., 2016) (delayed model), orange line (Limongi & Chieffi,
2018), and teal line (Griffith et al., 2021). Dashed lines and empty circles
have the same meaning as in Figure 1

Figure 8 shows a comparison between our core collapse yields and the results
of Pignatari et al. (2016); Limongi & Chieffi (2018) and Griffith et al.
(2021) for solar metallicity non-rotating single stars. Below
$\rm{M_{init}}\approx 20\,\rm{M}_{\odot}$ there is a reasonable agreement
between the different models, and software instruments, for the final core-
collapse $\mathrm{{}^{12}C}$ yield. Above this point the different models
begin to diverge, most likely due to changes in the behaviour of carbon shells
in each set of models and choice of nuclear reaction rates. Differences in the
choices for the treatment of stellar wind mass loss likely also plays a role.
We note that the spread in $\mathrm{{}^{12}C}$ yields is of the order of the
spread we find due to the differences in carbon shell burning (see section
4.1). Note the two open symbols, which indicate models in our grid where the
mixing due to convective carbon shell burning destroyed less carbon. These two
models are more consistent with the higher yields found by Limongi & Chieffi
(2018) and Griffith et al. (2021).

##  Appendix C Table of yields

In Table 2 we show the total $\mathrm{{}^{12}C}$ yields based on the source of
the mass loss. Several of the models show anomalous carbon burning (see
section 3.3), however the single $\rm{M_{init}}=33\,\rm{M}_{\odot}$ star shows
a slightly different behaviour. Here the convection zone above the carbon-
burning shell extended sufficiently far that the mixing region dredged down
all the $\mathrm{{}^{12}C}$ that existed in/near the helium shell. Thus by the
time of core-collapse the star becomes depleted in $\mathrm{{}^{12}C}$.

We also provide the integrated yields for single and binary-stripped stars for
different assumption with respect which stars explode successfully in Table 3.
The results are normalized to the assumption that all stars explode
successfully. We can see that changing the assumption for the successful
explosion has a large impact on the final $\mathrm{{}^{12}C}$ yields, with the
yields decreasing by 60–70%. The binary-stripped models show no differences
between the different envelope assumptions as the final helium core mass cut
and cut on initial masses are equivalent.

Table 2: $\mathrm{{}^{12}C}$ yields in solar masses broken down by mass-loss type for both single stars and the primary star in the binary. $\rm{M_{He,final}}$ is the final helium core mass at core-collapse .  Initial mass [$\,\rm{M}_{\odot}$] | Single [$\,\rm{M}_{\odot}$] | Binary [$\,\rm{M}_{\odot}$]
---|---|---
| $\rm{M_{He,final}}$ | $\rm{M_{Winds}}$ | $\rm{M_{CC}}$ | $\rm{M_{Total}}$ | $\rm{M_{He,final}}$ | $\rm{M_{RLOF}}$ | $\rm{M_{Winds}}$ | $\rm{M_{CC}}$ | $\rm{M_{Total}}$
11 | 3.798 | -0.001 | 0.088 | 0.087 | 3.191 | -0.006 | -0.001 | 0.055 | 0.051
12 | 4.297 | -0.002 | 0.100 | 0.098 | 3.602 | -0.007 | -0.001 | 0.073 | 0.069
13 | 4.800 | -0.002 | 0.104 | 0.102 | 4.014 | -0.007 | -0.002 | 0.074 | 0.071
14 | 5.308 | -0.003 | 0.114 | 0.111 | 4.422 | -0.007 | -0.002 | 0.084 | 0.081
15 | 5.829 | -0.004 | 0.135 | 0.132 | 4.832 | -0.007 | -0.002 | 0.134 | 0.132
16 | 6.359 | -0.004 | 0.142 | 0.137 | 5.222 | -0.008 | -0.003 | 0.157 | 0.156
17 | 6.896 | -0.005 | 0.174 | 0.169 | 5.608 | -0.008 | -0.003 | 0.198 | 0.198
18 | 7.414 | -0.006 | 0.191 | 0.185 | 5.980 | -0.008 | -0.004 | 0.235 | 0.235
19 | 7.817 | -0.007 | 0.190 | 0.184 | 6.243 | -0.008 | -0.004 | 0.226 | 0.227
20 | 8.275 | -0.007 | 0.164 | 0.157 | 6.579 | -0.009 | -0.005 | 0.290 | 0.291
21 | 8.382 | -0.007 | 0.197 | 0.190 | 6.793 | -0.010 | -0.005 | 0.249 | 0.251
22 | 7.943 | -0.009 | 0.190 | 0.181 | 8.365 | -0.011 | -0.002 | 0.134 | 0.133
23 | 9.366 | -0.009 | 0.221 | 0.212 | 8.246 | -0.011 | -0.004 | 0.233 | 0.234
24 | 9.856 | -0.010 | 0.469aaModel shows anomalous carbon burning behaviour (See section 3.3) | 0.458aaModel shows anomalous carbon burning behaviour (See section 3.3) | 9.476 | -0.008 | -0.004 | 0.178 | 0.182
25 | 9.530 | -0.011 | 0.160 | 0.149 | 9.403 | -0.010 | -0.004 | 0.250 | 0.256
26 | 9.995 | -0.011 | 0.176 | 0.164 | 10.636 | -0.010 | -0.005 | 0.185 | 0.191
27 | 10.785 | -0.013 | 0.198 | 0.185 | 10.283 | -0.009 | -0.005 | 0.187 | 0.196
28 | 11.154 | -0.012 | 0.441aaModel shows anomalous carbon burning behaviour (See section 3.3) | 0.429aaModel shows anomalous carbon burning behaviour (See section 3.3) | 10.889 | -0.009 | -0.006 | 0.189 | 0.201
29 | 12.194 | -0.014 | 0.188 | 0.173 | 11.085 | -0.009 | -0.010 | 0.613 | 0.631
30 | 12.449 | -0.014 | 0.182 | 0.167 | 11.094 | -0.010 | -0.011 | 0.719 | 0.737
31 | 13.303 | -0.015 | 0.188 | 0.173 | 13.131 | -0.011 | -0.007 | 0.167aaModel shows anomalous carbon burning behaviour (See section 3.3) | 0.182aaModel shows anomalous carbon burning behaviour (See section 3.3)
32 | 14.238 | -0.017 | 0.000bbModel did not reach core-collapse. | 0.000bbModel did not reach core-collapse. | 11.632 | -0.010 | 0.012 | 0.949 | 0.997
33 | 15.031 | -0.016 | 0.000aaModel shows anomalous carbon burning behaviour (See section 3.3) | -0.016aaModel shows anomalous carbon burning behaviour (See section 3.3) | 11.969 | -0.011 | -0.013 | 0.866 | 0.890
34 | 14.870 | -0.019 | 0.488 | 0.469 | 12.818 | -0.010 | -0.014 | 0.718 | 0.747
35 | 15.031 | -0.021 | 0.733 | 0.712 | 13.252 | -0.010 | 0.025 | 1.002 | 1.073
36 | 15.469 | -0.021 | 0.786 | 0.764 | 11.650 | -0.014 | 0.575 | 0.975 | 1.595
37 | 15.789 | -0.005 | 0.945 | 0.939 | 13.159 | -0.010 | 0.879 | 0.977 | 1.912
38 | 16.196 | 0.052 | 0.000bbModel did not reach core-collapse. | 0.000bbModel did not reach core-collapse. | 13.552 | -0.010 | 0.962 | 1.008 | 2.029
39 | 16.576 | 0.120 | 0.806 | 0.926 | 13.948 | -0.011 | 1.039 | 0.971 | 2.071
40 | 16.966 | 0.186 | 0.979 | 1.165 | 14.357 | -0.011 | 1.112 | 1.063 | 2.240
41 | 17.311 | 0.277 | 1.130 | 1.407 | 14.776 | -0.011 | 1.180 | 1.062 | 2.311
42 | 17.638 | 0.379 | 1.077 | 1.456 | 15.183 | -0.012 | 1.253 | 1.021 | 2.346
43 | 17.999 | 0.458 | 1.143 | 1.600 | 15.599 | -0.013 | 1.322 | 1.134 | 2.531
44 | 18.328 | 0.554 | 1.094 | 1.648 | 15.831 | -0.013 | 1.434 | 1.069 | 2.585
45 | 18.621 | 0.667 | 1.172 | 1.839 | 16.279 | -0.013 | 1.485 | 1.138 | 2.709
Table 3:  Impact of assumptions regarding black hole formation on the core-collapse IMF weighted yield predictions normalised to the default assumption that all stars explode successfully. We assume an IMF power-law exponent $\alpha=-2.3$. In case of BH formation we assume that the carbon-rich layers fall back onto the BH. BH formation assumption | Normalised $\mathrm{{}^{12}C}$ yields
---|---
| Single stars | Binary-stripped stars
All stars explode successfully | 1.0 | 1.0
BHs form from stars with initial masses $\rm{\rm{M_{init}}>22}\,\rm{M}_{\odot}$ | 0.41 | 0.34
BHs form from stars final core masses $\rm{M_{He,final}}>7\,\rm{M}_{\odot}$ | 0.27 | 0.34
Schneider et al 2021 like formation | 0.55 | 0.84

We acknowledge helpful discussions with B. Paxton, F. Timmes, L. van Son, and
T. Wagg. The authors acknowledge funding by the European Union’s Horizon 2020
research and innovation program from the European Research Council (ERC)
(Grant agreement No. 715063), and by the Netherlands Organization for
Scientific Research (NWO) as part of the Vidi research program BinWaves with
project number 639.042.728 and a top module 2 grant with project number
614.001.501. This work was also supported by the Cost Action Program ChETEC
CA16117. This work was carried out on the Dutch national e-infrastructure with
the support of SURF Cooperative. This research has made use of NASA’s
Astrophysics Data System.

## References

  * Abt (1983) Abt, H. A. 1983, ARA&A, 21, 343
  * Akerman et al. (2004) Akerman, C. J., Carigi, L., Nissen, P. E., Pettini, M., & Asplund, M. 2004, A&A, 414, 931, doi: 10.1051/0004-6361:20034188
  * Angulo et al. (1999) Angulo, C., Arnould, M., Rayet, M., et al. 1999, Nuclear Physics A, 656, 3, doi: 10.1016/S0375-9474(99)00030-5
  * Arnett & Truran (1969) Arnett, W. D., & Truran, J. W. 1969, ApJ, 157, 339, doi: 10.1086/150072
  * Aufderheide et al. (1991) Aufderheide, M. B., Baron, E., & Thielemann, F. K. 1991, ApJ, 370, 630, doi: 10.1086/169849
  * Bensby & Feltzing (2006) Bensby, T., & Feltzing, S. 2006, MNRAS, 367, 1181, doi: 10.1111/j.1365-2966.2006.10037.x
  * Bevan et al. (2017) Bevan, A., Barlow, M. J., & Milisavljevic, D. 2017, MNRAS, 465, 4044, doi: 10.1093/mnras/stw2985
  * Braun & Langer (1995) Braun, H., & Langer, N. 1995, in Wolf-Rayet Stars: Binaries; Colliding Winds; Evolution, ed. K. A. van der Hucht & P. M. Williams, Vol. 163, 305
  * Brinkman et al. (2019) Brinkman, H. E., Doherty, C. L., Pols, O. R., et al. 2019, ApJ, 884, 38, doi: 10.3847/1538-4357/ab40ae
  * Brooker et al. (2021) Brooker, E. S., Stangl, S. M., Mauney, C. M., & Fryer, C. L. 2021, arXiv e-prints, arXiv:2103.12781.  https://arxiv.org/abs/2103.12781
  * Brott et al. (2011) Brott, I., de Mink, S. E., Cantiello, M., et al. 2011, A&A, 530, A115, doi: 10.1051/0004-6361/201016113
  * Brown et al. (2001) Brown, G. E., Heger, A., Langer, N., et al. 2001, New A, 6, 457, doi: 10.1016/S1384-1076(01)00077-X
  * Brown et al. (1999) Brown, G. E., Lee, C. H., & Bethe, H. A. 1999, New A, 4, 313, doi: 10.1016/S1384-1076(99)00023-8
  * Brown & Woosley (2013) Brown, J. M., & Woosley, S. E. 2013, ApJ, 769, 99, doi: 10.1088/0004-637X/769/2/99
  * Buchler & Yueh (1976) Buchler, J. R., & Yueh, W. R. 1976, ApJ, 210, 440, doi: 10.1086/154847
  * Burbidge et al. (1957) Burbidge, E. M., Burbidge, G. R., Fowler, W. A., & Hoyle, F. 1957, Reviews of Modern Physics, 29, 547, doi: 10.1103/RevModPhys.29.547
  * Burton & Gordon (1978) Burton, W. B., & Gordon, M. A. 1978, A&A, 63, 7
  * Cameron (1959) Cameron, A. G. W. 1959, ApJ, 130, 429, doi: 10.1086/146734
  * Carigi et al. (2005) Carigi, L., Peimbert, M., Esteban, C., & García-Rojas, J. 2005, ApJ, 623, 213, doi: 10.1086/428491
  * Carilli & Walter (2013) Carilli, C. L., & Walter, F. 2013, ARA&A, 51, 105, doi: 10.1146/annurev-astro-082812-140953
  * Cassisi et al. (2007) Cassisi, S., Potekhin, A. Y., Pietrinferni, A., Catelan, M., & Salaris, M. 2007, ApJ, 661, 1094, doi: 10.1086/516819
  * Cescutti et al. (2009) Cescutti, G., Matteucci, F., McWilliam, A., & Chiappini, C. 2009, A&A, 505, 605, doi: 10.1051/0004-6361/200912759
  * Chugunov et al. (2007) Chugunov, A. I., Dewitt, H. E., & Yakovlev, D. G. 2007, Phys. Rev. D, 76, 025028, doi: 10.1103/PhysRevD.76.025028
  * Cristini et al. (2017) Cristini, A., Meakin, C., Hirschi, R., et al. 2017, MNRAS, 471, 279, doi: 10.1093/mnras/stx1535
  * Curtis et al. (2019) Curtis, S., Ebinger, K., Fröhlich, C., et al. 2019, ApJ, 870, 2, doi: 10.3847/1538-4357/aae7d2
  * Cyburt et al. (2010) Cyburt, R. H., Amthor, A. M., Ferguson, R., et al. 2010, ApJS, 189, 240, doi: 10.1088/0067-0049/189/1/240
  * De Donder & Vanbeveren (2004) De Donder, E., & Vanbeveren, D. 2004, New A Rev., 48, 861, doi: 10.1016/j.newar.2004.07.001
  * de Jager et al. (1988) de Jager, C., Nieuwenhuijzen, H., & van der Hucht, K. A. 1988, A&AS, 72, 259
  * de Mink et al. (2009) de Mink, S. E., Pols, O. R., Langer, N., & Izzard, R. G. 2009, A&A, 507, L1, doi: 10.1051/0004-6361/200913205
  * Dray & Tout (2003) Dray, L. M., & Tout, C. A. 2003, MNRAS, 341, 299, doi: 10.1046/j.1365-8711.2003.06420.x
  * Dray et al. (2003) Dray, L. M., Tout, C. A., Karakas, A. I., & Lattanzio, J. C. 2003, MNRAS, 973, doi: 10.1046/j.1365-8711.2003.06142.x
  * Eggleton (1983) Eggleton, P. P. 1983, ApJ, 268, 368, doi: 10.1086/160960
  * Ekström et al. (2012) Ekström, S., Georgy, C., Eggenberger, P., et al. 2012, A&A, 537, A146, doi: 10.1051/0004-6361/201117751
  * Eldridge & Tout (2004) Eldridge, J. J., & Tout, C. A. 2004, MNRAS, 353, 87, doi: 10.1111/j.1365-2966.2004.08041.x
  * Eldridge & Vink (2006) Eldridge, J. J., & Vink, J. S. 2006, A&A, 452, 295, doi: 10.1051/0004-6361:20065001
  * Farmer (2018) Farmer, R. 2018, rjfarmer/mesaplot, doi: 10.5281/zenodo.1441329
  * Farmer & Bauer (2018) Farmer, R., & Bauer, E. B. 2018, rjfarmer/pyMesa: Add support for 10398, v1.0.3, Zenodo, doi: 10.5281/zenodo.1205271
  * Farmer et al. (2016) Farmer, R., Fields, C. E., Petermann, I., et al. 2016, ApJS, 227, 22, doi: 10.3847/1538-4365/227/2/22
  * Ferguson et al. (2005) Ferguson, J. W., Alexander, D. R., Allard, F., et al. 2005, ApJ, 623, 585, doi: 10.1086/428642
  * Fields et al. (2018) Fields, C. E., Timmes, F. X., Farmer, R., et al. 2018, The Astrophysical Journal Supplement Series, 234, 19, doi: 10.3847/1538-4365/aaa29b
  * Franchini et al. (2020) Franchini, M., Morossi, C., Di Marcantonio, P., et al. 2020, ApJ, 888, 55, doi: 10.3847/1538-4357/ab5dc4
  * Frerking et al. (1982) Frerking, M. A., Langer, W. D., & Wilson, R. W. 1982, ApJ, 262, 590, doi: 10.1086/160451
  * Fuller et al. (1985) Fuller, G. M., Fowler, W. A., & Newman, M. J. 1985, ApJ, 293, 1, doi: 10.1086/163208
  * Götberg et al. (2017) Götberg, Y., de Mink, S. E., & Groh, J. H. 2017, A&A, 608, A11, doi: 10.1051/0004-6361/201730472
  * Grevesse & Sauval (1998) Grevesse, N., & Sauval, A. J. 1998, Space Sci. Rev., 85, 161, doi: 10.1023/A:1005161325181
  * Griffith et al. (2021) Griffith, E. J., Sukhbold, T., Weinberg, D. H., et al. 2021, arXiv e-prints, arXiv:2103.09837.  https://arxiv.org/abs/2103.09837
  * Groh et al. (2019) Groh, J. H., Ekström, S., Georgy, C., et al. 2019, A&A, 627, A24, doi: 10.1051/0004-6361/201833720
  * Gullberg et al. (2018) Gullberg, B., Swinbank, A. M., Smail, I., et al. 2018, ApJ, 859, 12, doi: 10.3847/1538-4357/aabe8c
  * Habets (1986) Habets, G. M. H. J. 1986, A&A, 165, 95
  * Henry et al. (2000) Henry, R. B. C., Edmunds, M. G., & Köppen, J. 2000, ApJ, 541, 660, doi: 10.1086/309471
  * Herbst & van Dishoeck (2009) Herbst, E., & van Dishoeck, E. F. 2009, ARA&A, 47, 427, doi: 10.1146/annurev-astro-082708-101654
  * Hirschi et al. (2004) Hirschi, R., Meynet, G., & Maeder, A. 2004, A&A, 425, 649, doi: 10.1051/0004-6361:20041095
  * Hirschi et al. (2005) —. 2005, A&A, 433, 1013, doi: 10.1051/0004-6361:20041554
  * Hunter et al. (2007) Hunter, I., Dufton, P. L., Smartt, S. J., et al. 2007, A&A, 466, 277, doi: 10.1051/0004-6361:20066148
  * Hunter et al. (2008) Hunter, I., Brott, I., Lennon, D. J., et al. 2008, ApJ, 676, L29, doi: 10.1086/587436
  * Hunter et al. (2009) Hunter, I., Brott, I., Langer, N., et al. 2009, A&A, 496, 841, doi: 10.1051/0004-6361/200809925
  * Hunter (2007) Hunter, J. D. 2007, Computing In Science & Engineering, 9, 90
  * Iglesias & Rogers (1993) Iglesias, C. A., & Rogers, F. J. 1993, ApJ, 412, 752, doi: 10.1086/172958
  * Iglesias & Rogers (1996) —. 1996, ApJ, 464, 943, doi: 10.1086/177381
  * Itoh et al. (1996) Itoh, N., Hayashi, H., Nishikawa, A., & Kohyama, Y. 1996, ApJS, 102, 411, doi: 10.1086/192264
  * Izzard (2004) Izzard, R. G. 2004, PhD thesis, University of Cambridge
  * Izzard et al. (2006) Izzard, R. G., Dray, L. M., Karakas, A. I., Lugaro, M., & Tout, C. A. 2006, A&A, 460, 565, doi: 10.1051/0004-6361:20066129
  * Karakas & Lugaro (2016) Karakas, A. I., & Lugaro, M. 2016, ApJ, 825, 26, doi: 10.3847/0004-637X/825/1/26
  * Kippenhahn & Weigert (1967) Kippenhahn, R., & Weigert, A. 1967, ZAp, 65, 251
  * Kluyver et al. (2016) Kluyver, T., Ragan-Kelley, B., Pérez, F., et al. 2016, in Positioning and Power in Academic Publishing: Players, Agents and Agendas: Proceedings of the 20th International Conference on Electronic Publishing, ed. F. Loizides & B. Schmidt, Amsterdam: IOS Press, 87
  * Kobayashi et al. (2006) Kobayashi, C., Umeda, H., Nomoto, K., Tominaga, N., & Ohkubo, T. 2006, ApJ, 653, 1145, doi: 10.1086/508914
  * Langanke & Martínez-Pinedo (2000) Langanke, K., & Martínez-Pinedo, G. 2000, Nuclear Physics A, 673, 481, doi: 10.1016/S0375-9474(00)00131-7
  * Langer (1989) Langer, N. 1989, A&A, 210, 93
  * Langer (1991) —. 1991, A&A, 248, 531
  * Langer (2012) —. 2012, ARA&A, 50, 107, doi: 10.1146/annurev-astro-081811-125534
  * Langer & Henkel (1995) Langer, N., & Henkel, C. 1995, Space Sci. Rev., 74, 343, doi: 10.1007/BF00751420
  * Laplace et al. (2020) Laplace, E., Götberg, Y., de Mink, S. E., Justham, S., & Farmer, R. 2020, A&A, 637, A6, doi: 10.1051/0004-6361/201937300
  * Laplace et al. (2021) Laplace, E., Justham, S., Renzo, M., et al. 2021, arXiv e-prints, arXiv:2102.05036.  https://arxiv.org/abs/2102.05036
  * Lau et al. (2020) Lau, R. M., Eldridge, J. J., Hankins, M. J., et al. 2020, ApJ, 898, 74, doi: 10.3847/1538-4357/ab9cb5
  * Leung & Nomoto (2020) Leung, S.-C., & Nomoto, K. 2020, ApJ, 888, 80, doi: 10.3847/1538-4357/ab5c1f
  * Li & Draine (2001) Li, A., & Draine, B. T. 2001, ApJ, 554, 778, doi: 10.1086/323147
  * Limongi & Chieffi (2018) Limongi, M., & Chieffi, A. 2018, ApJS, 237, 13, doi: 10.3847/1538-4365/aacb24
  * Maeder (1983) Maeder, A. 1983, A&A, 120, 113
  * Maeder (1992) —. 1992, A&A, 264, 105
  * Maeder et al. (2014) Maeder, A., Przybilla, N., Nieva, M.-F., et al. 2014, A&A, 565, A39, doi: 10.1051/0004-6361/201220602
  * Mason et al. (2009) Mason, B. D., Hartkopf, W. I., Gies, D. R., Henry, T. J., & Helsel, J. W. 2009, AJ, 137, 3358
  * Meynet et al. (2006) Meynet, G., Ekström, S., & Maeder, A. 2006, A&A, 447, 623, doi: 10.1051/0004-6361:20053070
  * Meynet & Maeder (2000) Meynet, G., & Maeder, A. 2000, A&A, 361, 101.  https://arxiv.org/abs/astro-ph/0006404
  * Moe & Di Stefano (2017) Moe, M., & Di Stefano, R. 2017, ApJS, 230, 15
  * Murphy et al. (2021) Murphy, L. J., Groh, J. H., Ekström, S., et al. 2021, MNRAS, 501, 2745, doi: 10.1093/mnras/staa3803
  * Nissen et al. (2014) Nissen, P. E., Chen, Y. Q., Carigi, L., Schuster, W. J., & Zhao, G. 2014, A&A, 568, A25, doi: 10.1051/0004-6361/201424184
  * Nomoto et al. (2006) Nomoto, K., Tominaga, N., Umeda, H., Kobayashi, C., & Maeda, K. 2006, Nucl. Phys. A, 777, 424, doi: 10.1016/j.nuclphysa.2006.05.008
  * Nugis & Lamers (2000) Nugis, T., & Lamers, H. J. G. L. M. 2000, A&A, 360, 227
  * O’Connor & Ott (2011) O’Connor, E., & Ott, C. D. 2011, ApJ, 730, 70
  * Oda et al. (1994) Oda, T., Hino, M., Muto, K., Takahara, M., & Sato, K. 1994, Atomic Data and Nuclear Data Tables, 56, 231, doi: 10.1006/adnd.1994.1007
  * Paczyński (1967) Paczyński, B. 1967, Acta Astron., 17, 355
  * Paxton et al. (2011) Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3, doi: 10.1088/0067-0049/192/1/3
  * Paxton et al. (2013) Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4, doi: 10.1088/0067-0049/208/1/4
  * Paxton et al. (2015) Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15, doi: 10.1088/0067-0049/220/1/15
  * Paxton et al. (2018) Paxton, B., Schwab, J., Bauer, E. B., et al. 2018, ApJS, 234, 34, doi: 10.3847/1538-4365/aaa5a8
  * Paxton et al. (2019) Paxton, B., Smolec, R., Schwab, J., et al. 2019, ApJS, 243, 10, doi: 10.3847/1538-4365/ab2241
  * Pérez & Granger (2007) Pérez, F., & Granger, B. E. 2007, Computing in Science & Engineering, 9, 21
  * Pignatari et al. (2016) Pignatari, M., Herwig, F., Hirschi, R., et al. 2016, ApJS, 225, 24, doi: 10.3847/0067-0049/225/2/24
  * Podsiadlowski et al. (1992) Podsiadlowski, P., Joss, P. C., & Hsu, J. J. L. 1992, ApJ, 391, 246, doi: 10.1086/171341
  * Pols et al. (1995a) Pols, O. R., Tout, C. A., Eggleton, P. P., & Han, Z. 1995a, MNRAS, 274, 964, doi: 10.1093/mnras/274.3.964
  * Pols et al. (1995b) —. 1995b, MNRAS, 274, 964, doi: 10.1093/mnras/274.3.964
  * Potekhin & Chabrier (2010) Potekhin, A. Y., & Chabrier, G. 2010, Contributions to Plasma Physics, 50, 82, doi: 10.1002/ctpp.201010017
  * Prantzos et al. (2018) Prantzos, N., Abia, C., Limongi, M., Chieffi, A., & Cristallo, S. 2018, MNRAS, 476, 3432, doi: 10.1093/mnras/sty316
  * Renzo et al. (2017) Renzo, M., Ott, C. D., Shore, S. N., & de Mink, S. E. 2017, A&A, 603, A118, doi: 10.1051/0004-6361/201730698
  * Ritter (1988) Ritter, H. 1988, A&A, 202, 93
  * Rogers & Nayfonov (2002) Rogers, F. J., & Nayfonov, A. 2002, ApJ, 576, 1064, doi: 10.1086/341894
  * Romano et al. (2020) Romano, D., Franchini, M., Grisoni, V., et al. 2020, A&A, 639, A37, doi: 10.1051/0004-6361/202037972
  * Romano et al. (2017) Romano, D., Matteucci, F., Zhang, Z. Y., Papadopoulos, P. P., & Ivison, R. J. 2017, MNRAS, 470, 401, doi: 10.1093/mnras/stx1197
  * Sana et al. (2012) Sana, H., de Mink, S. E., de Koter, A., et al. 2012, Science, 337, 444, doi: 10.1126/science.1223344
  * Sander & Vink (2020) Sander, A. A. C., & Vink, J. S. 2020, MNRAS, 499, 873, doi: 10.1093/mnras/staa2712
  * Sarangi et al. (2018) Sarangi, A., Matsuura, M., & Micelotta, E. R. 2018, Space Sci. Rev., 214, 63, doi: 10.1007/s11214-018-0492-7
  * Saumon et al. (1995) Saumon, D., Chabrier, G., & van Horn, H. M. 1995, ApJS, 99, 713, doi: 10.1086/192204
  * Sawada & Maeda (2019) Sawada, R., & Maeda, K. 2019, ApJ, 886, 47, doi: 10.3847/1538-4357/ab4da3
  * Schneider et al. (2021) Schneider, F. R. N., Podsiadlowski, P., & Müller, B. 2021, A&A, 645, A5, doi: 10.1051/0004-6361/202039219
  * Schneider et al. (2018) Schneider, F. R. N., Sana, H., Evans, C. J., et al. 2018, Science, 359, 69, doi: 10.1126/science.aan0106
  * Solomon et al. (1987) Solomon, P. M., Rivolo, A. R., Barrett, J., & Yahil, A. 1987, ApJ, 319, 730, doi: 10.1086/165493
  * Sukhbold & Adams (2020) Sukhbold, T., & Adams, S. 2020, MNRAS, 492, 2578, doi: 10.1093/mnras/staa059
  * Sukhbold et al. (2016) Sukhbold, T., Ertl, T., Woosley, S. E., Brown, J. M., & Janka, H. T. 2016, ApJ, 821, 38, doi: 10.3847/0004-637X/821/1/38
  * Sukhbold & Woosley (2014) Sukhbold, T., & Woosley, S. E. 2014, ApJ, 783, 10, doi: 10.1088/0004-637X/783/1/10
  * Suwa et al. (2019) Suwa, Y., Tominaga, N., & Maeda, K. 2019, MNRAS, 483, 3607, doi: 10.1093/mnras/sty3309
  * Thielemann et al. (1996) Thielemann, F.-K., Nomoto, K., & Hashimoto, M.-A. 1996, ApJ, 460, 408, doi: 10.1086/176980
  * Timmes & Swesty (2000) Timmes, F. X., & Swesty, F. D. 2000, ApJS, 126, 501, doi: 10.1086/313304
  * Timmes et al. (1996) Timmes, F. X., Woosley, S. E., & Weaver, T. A. 1996, ApJ, 457, 834, doi: 10.1086/176778
  * Tout et al. (1996) Tout, C. A., Pols, O. R., Eggleton, P. P., & Han, Z. 1996, MNRAS, 281, 257, doi: 10.1093/mnras/281.1.257
  * Townsend (2019) Townsend, R. 2019, MESA SDK for Linux: 20190911, doi: 10.5281/zenodo.2603170
  * van der Walt et al. (2011) van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science Engineering, 13, 22, doi: 10.1109/MCSE.2011.37
  * Vartanyan et al. (2021) Vartanyan, D., Laplace, E., Renzo, M., et al. 2021, arXiv e-prints, arXiv:2104.03317.  https://arxiv.org/abs/2104.03317
  * Vink et al. (2001) Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2001, A&A, 369, 574, doi: 10.1051/0004-6361:20010127
  * Weaver & Woosley (1993) Weaver, T. A., & Woosley, S. E. 1993, Phys. Rep., 227, 65, doi: 10.1016/0370-1573(93)90058-L
  * Weingartner & Draine (2001) Weingartner, J. C., & Draine, B. T. 2001, ApJ, 548, 296, doi: 10.1086/318651
  * Wellstein & Langer (1999) Wellstein, S., & Langer, N. 1999, A&A, 350, 148.  https://arxiv.org/abs/astro-ph/9904256
  * Wolfire et al. (1995) Wolfire, M. G., Hollenbach, D., McKee, C. F., Tielens, A. G. G. M., & Bakes, E. L. O. 1995, ApJ, 443, 152, doi: 10.1086/175510
  * Wolfire et al. (2003) Wolfire, M. G., McKee, C. F., Hollenbach, D., & Tielens, A. G. G. M. 2003, ApJ, 587, 278, doi: 10.1086/368016
  * Woosley (2019) Woosley, S. E. 2019, ApJ, 878, 49, doi: 10.3847/1538-4357/ab1b41
  * Woosley et al. (1993a) Woosley, S. E., Langer, N., & Weaver, T. A. 1993a, ApJ, 411, 823, doi: 10.1086/172886
  * Woosley et al. (1993b) Woosley, S. E., Timmes, F. X., & Weaver, T. A. 1993b, Journal of Physics G Nuclear Physics, 19, S183, doi: 10.1088/0954-3899/19/S/016
  * Woosley & Weaver (1995) Woosley, S. E., & Weaver, T. A. 1995, ApJS, 101, 181, doi: 10.1086/192237
  * Yoon et al. (2017) Yoon, S.-C., Dessart, L., & Clocchiatti, A. 2017, ApJ, 840, 10, doi: 10.3847/1538-4357/aa6afe
  * Yoon et al. (2010) Yoon, S. C., Woosley, S. E., & Langer, N. 2010, ApJ, 725, 940, doi: 10.1088/0004-637X/725/1/940
  * Young & Fryer (2007) Young, P. A., & Fryer, C. L. 2007, ApJ, 664, 1033, doi: 10.1086/518081
  * Zapartas et al. (2021) Zapartas, E., Renzo, M., Fragos, T., et al. 2021, arXiv e-prints, arXiv:2106.05228.  https://arxiv.org/abs/2106.05228