The cosmic carbon footprint of massive stars stripped in binary systems

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.


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 r.j.farmer@uva.nl(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 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. (2020Laplace et al. ( , 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 ( 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. 2011Paxton et al. , 2013Paxton et al. , 2015Paxton et al. , 2018Paxton et al. , 2019) ) to evolve massive single and binary stars from the zeroage main sequence to core collapse.Our single stars and the primary (initially most massive star) in the binary have initial masses between M init = 11 -45 M .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 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 (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. (2020Laplace et al. ( , 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 T eff > 10 4 K and surface hydrogen mass fraction X H > 0.4, Nugis & Lamers (2000) for T eff > 10 4 K and surface 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 convec-tion 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 α 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.netwhich 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 corecollapse, to compute the core structure accurately.However, we show in Appendix B.1 that models computed with approx21.netpredict similar 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, X He > 0.1 and the hydrogen mass fraction (at the same mass coordinate) is X H < 0.01.Core oxygen depletion is defined when the oxygen mass fraction at the center of the star drops below X O < 10 −4 .Finally, we define core collapse to occur when the inner regions of the star infalls at 300km s −1 .

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 M , 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 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 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 M 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 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×10 51erg/s at shock breakout while the kinetic energy at shock breakout is between 0.5-1.3×10 51erg/s.
We define the yield of an isotope as (Karakas & Lugaro 2016): Where ∆M T is the mass lost over the time interval T, X j is the surface mass fraction of isotope j, and X j,int is the initial mass fraction of isotope 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 X C12,int = 0.00244.

TOTAL CARBON
Figure 1(a) shows the total 12 C yields from our single and binary stars, from all sources of mass loss.Figure 1(b) shows the 12 C yield from winds, while Figure 1(c) shows the 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 12 C yield from the mass loss during RLOF is negative and small ≈ −0.01 M 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 CNOcycling and show depleted carbon abundances (Maeder 1983).The most massive stars show a slight decrease in 12 C ejected during RLOF as a larger fraction of the envelope is processed by CN(O)-cycling.
Figure 1(b) shows the 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 M init 35 M their winds are not 12 C enriched as compared to their initial composition, and thus not visible in Figure 1(b).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 12 C yield, taking into account reasonable assumptions for which stars eject their envelopes (Sukhbold et al. 2016;Zapartas et al. 2021).
Stars with M init 35 M have 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 M init = 37 M for single stars and M init = 35 M 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 12 C in their stellar winds than they started with.Figure 1(c) shows the 12 C yields for core-collapse ejecta, assuming that all stars eject their envelopes.The 12 C yield is relatively flat as a function of initial mass for stars with M init 27 M , at ≈ 0.2 M .Above this initial mass the yield rapidly increases up to ≈ 1.25 M .This transition between low and high carbon yields occurs between 27 M init / M 35.The increased 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).

Shock nucleosynthesis
Figure 2 shows the distribution of 12 C inside a binarystripped star with M init = 15 M 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 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 12 C near the inner boundary, and cools from a peak temperature of ≈ 10 10 K, to ≈ 3 × 10 8 K when it reaches the carbon-rich helium shell, and to ≈ 10 7 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 ≈ 1% for 12 C.The main effect on the 12 C yield is to smear the 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 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.

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 binarystripped) of the same initial mass.Figure 3 shows the time evolution of a single and a binary-stripped star with M init = 19, 37, and 45 M .Considering the M init = 19 M case, both stars start on the main se-quence 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 (τ cc − τ ) /yr ≈ 5.9.Most of the binary star's envelope is lost at this point, bringing its mass down to ≈ 8.0 M , 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 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 M init = 19 M single star the mass in the convective helium-burning core stays constant.This receding convective core leaves behind 12 C that was produced by the 3α process but had not yet been converted into oxygen in 12 C (α, γ) 16 O (Langer 1991).This left-over 12 C is outside of what will become the CO core of the star (denoted by the red lines).
Any 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 12 C + 12 C burning, or will be accreted onto the compact object during corecollapse.Thus only the 12 C that is mixed out of the core has a chance to survive until core-collapse.Some of this outwardly mixed 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 M init = 35 M .Here the wind mass loss is sufficient to remove the remaining helium layers above the core.This exposes 12 C rich material, which is then ejected in a wind, similar to the process in the binary-stripped star with M init = 19 M .Thus the single star has a net positive 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.the 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.

Carbon shell burning
In Figure 1 there is significant noise in the 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 M init = 24, 28, and 33 M for the single stars and M init = 31 M 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  Adams 2020).Figure 4 show the ignition of carbon in a 23 M single star.This star ignites 12 C at the center at log 10 (τ cc − τ ) /yr ≈ 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 12 C burning zone ignites off-center, at log 10 (τ cc − τ ) /yr ≈ 0.0 at the same time as the core ignites 20 Ne.It is this burning zone that causes the variability in the 12 C yields.
As this convection zone extends outwards, it mixes 12 C from the outer layers of the core inwards where the 12 C is then burnt.Thus the maximal extent of this zone and the amount of time that it has to mix 12 C downwards sets the final 12 C yields. Figure 4 shows that in this model a small pocket of 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 12 C and 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 ≈ 50% of the final 12 C yield, with the remaining 12 C yield coming from 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.

Sensitivity to physics choices
The sensitivity to the size and timing of the carbonburning 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 M 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 M 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.
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 ≈ 100, 000 time steps (for just carbon and oxygen burning) and ≈ 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 12 C and 16 O abundances at the center of the star.
Figure 5 shows the relative change in the total mass of 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 ≈ 2% difference in the total mass of 12 C (measured at core oxygen depletion).
Figure 5(a) shows the relative change in the total mass of 12 C compared to the relative change in the average time step taken.As the time resolution increases the total mass of 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 12 C (≈ 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 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 12 C to decrease, and the spread in 12 C values to decrease, as the amount of overshoot increases.However, it is possible to achieve the same final total mass of 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 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 12 C into the core.We note however that both predictive mixing and convective premixing assume that the mixing timescale is shorter than the computational timestep, which breaks down during carbon burning.

Sensitivity to the explosion properties
To explore the sensitivity of our results to the assumptions made we perform two sets of tests on a M init = 16 M 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 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 12 C yield.A horizontal line implies that the 12 C yield is insensitive to that parameter.Negative (M c12 − M c12,def ) /M c12,def values indicate that more 12 C was destroyed than in our default model, while positive values indicate that less 12 C was destroyed.
In Figure 6(a) and 6(b) we randomly sampled ≈ 200 times both the injected energy (between 0.5 -5 × 10 51 erg/s) and the mass cut (between 1.4 -2.0 M ). Figure 6(a) shows the effect of varying the injected energy on the total mass of 12 C ejected during the supernovae.We can see a strong correlation between injected energy and 12 C ejected, increasing the energy decreases the 12 C mass ejected.However, the change is < 1% in the total mass of 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 12 C for this range of mass cuts.These trends can be explained as follows: most of the 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 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 12 C yield due to the shock processing the small amount of 12 C that exists near the inner boundary (see Figure 2).However by the time the shock reaches the bulk of the 12 C in the helium shell the shock has cooled below 10 8 K, thus it can no longer burn the 12 C present.
In Figure 6(c) and 6(d) we randomly sampled ≈ 200 times both the time over which we inject energy during the explosion (10 −3 < T inj /s < 10 0 ) and the mass over which we inject energy (10 −2 < ∆M inj / M < 5×10 −1 ).When either the T inj is small (T inj < 0.1s) or ∆M inj is small (M inj < 0.1 M ) their effects on the 12 C yield is small.Only once M inj exceeds 0.1 M does it begin to dominate the 12 C yield.However, once T inj > 0.3s it begins to dominate the 12 C yield instead.When either T inj or ∆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 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 ≈ 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 soundcrossing 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 12 C mass of 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 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 12 C ejected is insensitive to the physical explosion parameters considered here, the amount of 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).

IMF WEIGHTED YIELDS
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 12 C yields for binary-stripped massive stars and single massive stars at solar metallicity.We define the ratio as: (2) Where we assumed a Salpeter-like IMF with different values of α (Schneider et al. 2018), M init,b and M init,s are the initial mass of the primary star in the binary and the mass of the single star in solar masses, Y b and Y s are the yields of 12 C in solar masses for the respective stars and for each type of mass loss, and f b and 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 ≈ 40% more 12 C compared to the same initial mass of single massive stars, assuming a standard Salpeter IMF α = −2.3.This is due to extra mass loss in binary systems, leading both to a higher 12 C yield in the winds and in the final supernovae for M init ≈ 30 -40 M 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 M init 22 M 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 M He 7 M .
Assuming all stars eject their envelopes with M init < 22 M , then stripped primary stars in massive binaries contribute even more to the 12 C in the Universe than single massive stars.
Next, taking the helium core mass at M He,final < 7 M , measured at core collapse, binary-stripped mas-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 α, 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 12 C yields (binary-stripped/single) sive stars eject approximately twice as much carbon than single massive stars (for α = −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 M He,final < 7 M is equivalent to a cut of M init < 19 M for single stars and M init < 22 M 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 M init ≤ 21.5 M and for stars which were stripped as case B binaries with M init ≤ 31.5 M .Then a second population of NSs are formed from more massive objects, for single stars in the range 23.5 ≤ M init / M ≤ 34.0 and for case B binaries with 34.0 ≤ M init / M ≤ 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. 1999Brown et al. , 2001)).
As the IMF α increases, thus favouring the production of more massive stars, the contribution that binarystripped makes to the 12 C yields increases.This is due the greater weight now given to the 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 binarystripped stars having 12 C positive yields at lower initial masses, and due to the higher 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 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 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 12 C.
We can combine the result shown in Figure 1(b), indicating that the 12 C yields from winds of massive stars are not positive until M init ≈ 35 M , with the expectation that only massive stars with M init 22 M eject their envelopes at core collapse.This combination suggests that stars with initial masses 22 M init / M 35 do not contribute to the net production of 12 C in the Universe.The contribution of massive stars with M init < 22 M can be well described by only their 12 C supernovae yield, while stars with M init > 35 M are well described by only their wind yields.

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 M =10 -15 M 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 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 12 C in the helium shell.This would decrease the significance of binary-stripped stars on the total 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 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 12 C above it.Thus the winds from the binary-stripped star will not be enhanced in 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 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 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(Hunter et al. , 2008(Hunter et al. , 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 12 C of a approximately a factor of three, between initial rotation rates of 0 -300km s −1 for a 20 M star.Hirschi et al. (2004Hirschi et al. ( , 2005) ) found that the 12 C yields increase by a factor 1.5 -2.5 for stars M init < 30 M , and with initial rotation rates of 0 -300km s −1 due to the more massive core that the rotating models produce.Fields et al. (2018) explore the sensitivity of a 15 M 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 ±80% compared to the model which adopts median reaction rates, driven by variations in 12 C (α, γ) 16 O and the triple-α rates.
Sukhbold & Adams (2020) showed how the type of carbon burning, radiative or convective, is altered by changing the 12 C (α, γ) 16 O rate, in 14 -26 M 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(Brown et al. , 2001)), and thus whether it will eject its envelope (Weaver & Woosley 1993;Timmes et al. 1996).Sukhbold & Adams (2020) found that variations in the ≈ 1σ uncertainty in the 12 C (α, γ) 16 O can move the location between convective and non-convection carbon ignition by ±2 M 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.

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 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 12 C yield for different mass loss processes.
Our results can be summarised as follows: 1. We find that massive stars stripped in binaries during hydrogen-shell burning are nearly twice as efficient at contributing to the cosmic 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 2and Table 1).
2. We confirm that the difference in yields between binary-stripped and single massive stars can be explained by considering the outer-most 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 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 12 C-rich layers disconnect and form a pocket where the temperature never becomes high enough for 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. Stellar winds also eject 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 > 36M and > 38M for our solarmetallicity 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. Mass loss from binary systems during conservative mass transfer does not contribute to 12 C yields, in our models.These layers are either pristine or contain CNO-processed material which is 12 C poor (c.f.de Mink et al. 2009).
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 ≈ 0.5% (See Figure 6).
6.We show that the 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 ≈ 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. We note that the variations between the predictions of 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.2021) were evolved with mesa 128.net after core oxygen depletion, and then exploded with the same method as in Section 2.

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).In Figure 7 we compare the 12 C mass ejected in the supernova calculations presented in Figure 1 to the 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.netafter 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 12 C ejected in the supernovae is similar between both sets.The largest disagreement occurs for the binary-stripped stars at M init = 21 M .This is due to the carbon shell destroying more 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 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 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 12 C that is ejected in the supernova resides in/near the helium shell.Thus the total 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 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 M init ≈ 20 M there is a reasonable agreement between the different models, and software instruments, for the final corecollapse 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 12 C yields is of the order of the spread we find due to the differences in carbon shell burning (see section 4.1).(Griffith et al. 2021).Dashed lines and empty circles have the same meaning as in Figure 1 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).

C. TABLE OF YIELDS
In Table 2 we show the total 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 M init = 33 M star shows a slightly different behaviour.Here the convection zone above the carbonburning shell extended sufficiently far that the mixing region dredged down all the 12 C that existed in/near the helium shell.Thus by the time of core-collapse the star becomes depleted in 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 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. 12 C yields in solar masses broken down by mass-loss type for both single stars and the primary star in the binary.M He,final is the final helium core mass at core-collapse .

Initial mass [ M ]
Single a Model shows anomalous carbon burning behaviour (See section 3.3) b Model did not reach core-collapse.
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 α = −2.3.In case of BH formation we assume that the carbon-rich layers fall back onto the BH.

BH formation assumption
Normalised

Figure 1 .
Figure 1.Top panel: The total 12 C yield from all ejection sources, Left Panel: The 12 C yield ejected during wind mass loss, Right Panel: The 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 12 C yield, taking into account reasonable assumptions for which stars eject their envelopes(Sukhbold et al. 2016;Zapartas et al. 2021).

Figure 2 .
Figure 2. The composition profile of a 15 M binarystripped star at the start and end of the shock propagation phase of core collapse.The colors denote the mass fraction of 12 C.The animated version shows the shock propagation and resulting nucleosynthesis with time.

Figure 3 .
Figure 3. Kippenhahn diagrams for the inner regions of single and binary-stripped stars with Minit = 19, 37, and 45 M , 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 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.

Figure 4 .
Figure 4. Kippenhahn diagram for a 23 M single star during carbon and oxygen burning.Colors show the mass fraction of 12 C at each mass co-ordinate.Hatching shows mixing regions due to convection and overshoot.

Figure 5 .
Figure5.The relative change in the final total mass of 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 M 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 f0.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 7 .
Figure 7.Total 12 C mass ejected during supernovae, comparing models from this work (dashed) with those from Laplace et al. (2021) (solid).In blue are plotted singlestar models, while red represents binary-stripped models.The models ofLaplace et al. (2021) were evolved with mesa 128.net after core oxygen depletion, and then exploded with the same method as in Section 2.
Laplace et al. 2021 . 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. 7 Note the two open symbols, which indicate models in our grid where the