© ESO 2020

1. Introduction

The formation of gravitational-wave (GWs) sources is a key problem that is becoming increasingly relevant for discussion in the new era of gravitational-wave detections (Abbott et al. 2016, 2017, 2019). The compact objects, neutron stars (NSs), or black holes (BHs), whose mergers give rise to the gravitational wave chirp, represent the end products of massive stars (above about 8 M_⊙).

The way in which remnants of two stars can eventually end up in an orbit that is close enough for them to merge through the emission of gravitational waves within a Hubble time is a process that includes many aspects which are still poorly understood. Not only does it require a detailed understanding of the evolution and fate of massive stars, but also of their binary interaction. Moreover, with the reach of present-day GW detectors, we can probe nearby mergers of compact objects, but we expect their progenitors to have formed at appreciable or even large redshifts. This is because of the time delay between the formation of a double compact objects and the final merger. Thus, we expect that many of the progenitors formed out of more pristine gas, that is, gas that is less enriched by heavy elements that have been contributed by previous generations of stars. This means we must carefully understand the effect that metallicity has on the evolution of massive stars in binaries.

Several scenarios for the formation of gravitational-wave sources have been proposed. Many of these involve the interaction between two stars in a close binary system through one or multiple phases of Roche-lobe overflow (e.g., Kippenhahn & Weigert 1967; Tutukov et al. 1973; Tauris et al. 2017, and references therein) that end up stripping one star, or eventually both stars, of most of the hydrogen-rich envelope, which is about two thirds of the initial mass. The star that results from this process is what we refer to as a stripped star, hereafter.

Stripped stars are largely composed of helium and, eventually, heavier elements. As a result, one may naively expect these stars to be very compact. This is largely the case, at least during the long-lived phase of core helium fusion. However, it has been shown that stripped stars can swell and reach giant dimensions in the late stages of their evolution, depending on their mass. Some of the early numerical calculations already demonstrated this phenomenon, approximating stripped stars as pure helium stars (e.g., Divine 1965; Habets 1986a). More recent calculations confirmed this either by considering pure helium stars (Dewi et al. 2002; Dewi & Pols 2003; Ivanova et al. 2003) or by fully following the evolution of the massive star progenitor through the stripping process in a binary system (e.g., Yoon et al. 2010; Eldridge et al. 2013; Sravan et al. 2019).

Whether or not a stripped star expands and by how much is relevant for understanding their fate as the progenitors of core-collapse supernovae and possibly gravitational-wave sources. The large size of stripped stars implies that they can fill their Roche lobe anew and undergo an additional phase of mass transfer (Dewi et al. 2002; Dewi & Pols 2003; Ivanova et al. 2003). Additional phases of mass transfer can produce stars with even lower envelope masses, known as ultra-stripped stars (Tauris et al. 2015). When these stars end their lives, it is believed that the resulting low ejecta masses lead to very small supernova kicks (e.g., Podsiadlowski et al. 2004), which can prevent disruption of the binary system at the moment of explosion. This favors the formation of close NS binaries, some of which will be tight enough to merge as a result of GWs within a Hubble time (Tauris et al. 2017). The occurrence and outcome of additional phases of mass transfer thus directly impact predictions for the formation of gravitational-wave sources, particularly double neutron star mergers.

Improving our understanding of the expansion of stripped stars is also relevant in the light of upcoming electromagnetic transient surveys because they are responsible for about a third of all core-collapse supernovae (Graur et al. 2017a,b; Modjaz et al. 2019). The radius and the mass of their envelopes impact the light-curves predicted for stripped-envelope supernovae (Kleiser & Kasen 2014; Dessart et al. 2018; Kleiser et al. 2018), together with the mass-loss rate expected at late times and, thus, the circumstellar material around the progenitor at the moment of explosion (Ouchi & Maeda 2017). The case of type Ib supernova iPTF13bvn (Cao et al. 2013) is of particular interest as it appears to provide the most direct evidence we have that stripped stars can end their lives as helium giants (e.g., Fremling et al. 2014; Bersten et al. 2014; Eldridge & Maund 2016).

Recently, several studies have drawn attention to the fact that metallicity can have a large impact on the properties of stripped stars. Götberg et al. (2017) and Yoon et al. (2017) find that metallicity strongly affects the radial extent of the hydrogen layer that is left at the surface of the star after stripping. At low metallicity the reduced internal opacity makes it possible for stars to retract within their Roche lobe before the hydrogen envelope has been fully removed; see also Götberg et al. (2017). In addition, metallicity affects the stellar winds, which can strip the stars even further (see also Vink et al. 2001; Gilkis et al. 2019).

Calculations of the formation of gravitational-wave progenitors are typically made with binary population synthesis simulations, which rely on simplified assumptions for the stellar structure and interaction phases. Such simplifications are necessary because the simulations typically involve following the evolution of millions of stars in binary systems through complex phases of interaction. In the vast majority of binary population synthesis codes stars are treated with analytic fits by Hurley et al. (2000) against evolutionary models by Pols et al. (1998). A list of examples is given in Sect. 5.1. An alternative is to interpolate in pre-computed grids of models. Two examples of recent studies that use grids of pre-computed single-star models are Kruckow et al. (2018) and Giacobbo & Mapelli (2018); see also Sect. 5.2. A third alternative is to post-process extended grids of binary evolutionary models. The most prominent example of this is the BPASS code (Eldridge & Stanway 2016; Eldridge et al. 2017).

For the treatment of stripped stars, all population synthesis studies listed above (with the exception of BPASS) use interpolation in grids of models, or fits against them, of single helium stars computed at solar metallicity. They do not make use of models where the stripped star has been computed self-consistently through the Roche stripping process. This has two drawbacks: (1) they do not account for the effect that a leftover layer of hydrogen on the surface may have on the properties of stripped stars; (2) they do not fully account for the effect that metallicity may have on the properties of stripped stars.

In this paper, we present a study of the radius evolution of stars stripped in massive binaries considering solar and low metallicity. For this purpose, we compute a grid of representative progenitor models for different masses that are relevant as supernova and possible neutron star progenitors. We follow their evolution through Roche-lobe overflow with a detailed binary evolution code. We focus on the expansion phases after central helium depletion and discuss how this is linked with their interior evolution and in particular the burning phases. We then show how the radii compare to sizes usually assumed in binary population synthesis models and estimate how the differences impact the number of systems that can interact a second time through Roche-lobe overflow. We estimate and discuss the implications for core-collapse supernovae and gravitational-wave progenitors.

This paper can be considered as a companion of the paper by Götberg et al. (2017), which presents an extensive discussion of the effect of metallicity on the long-lived phase of central helium burning. In this work, we extend this study to the later evolutionary phases. The paper is structured as follows: we summarize our model assumptions in Sect. 2 and we discuss the effect of metallicity using two representative stellar models in Sect. 3. In Sect. 4, we present our full grids of evolutionary stellar models. We compare the radii obtained to those commonly used in population synthesis models in Sect. 5. We then discuss the impact of these large radii on the progenitors of supernova and GW sources in Sect. 6 together with a discussion of the uncertainties. A summary with our conclusions is provided in Sect. 7.

2. Binary evolution models

For our calculations of the interacting binary stars, we used the open-source 1D stellar evolution code MESA (version 10398, Paxton et al. 2011, 2013, 2015, 2018, 2019). The models are computed at solar metallicity (initial metal fraction of Z_⊙ ≡ 0.0142, based on values from Asplund et al. (2009) and sub-solar metallicity (Z = 0.001, representative of nearby low-metallicity environments, such as the Small Magellanic Cloud, Z_SMC ≊ Z_⊙/5).

Our zero age main sequence (ZAMS) models were computed by following their pre-main sequence evolution until the central helium abundance has increased by 5%. Following Tout et al. (1996), which is also consistent with Pols et al. (1998), we assumed an initial hydrogen mass fraction of X = 1 − Z − Y, where Y = 0.24 + 2Z is the helium mass fraction. Strictly speaking, given the updated abundances from Asplund et al. (2009), the helium abundances should be adapted, even though the difference is small. We computed the evolution of the stars until core carbon depletion (defined as the moment when the core carbon abundance decreases below 10⁻⁴). This is sufficient for the purpose of this work since the remaining evolutionary time is so short (less than 100 yr in these models) that the outer layers do not have time to react to quasi-hydrostatic changes in the core.

We used a nuclear network comprising 21 isotopes which follows the most prominent nuclear processes that influence the life of massive stars from hydrogen burning through the CNO cycle until silicon burning with sufficient accuracy (approx21, Paxton et al. 2011). It has the advantage of enabling fast calculations while containing the most important isotopes relevant for our study. We used the default opacity tables of MESA (Iglesias & Rogers 1993, 1996; Buchler & Yueh 1976; Cassisi et al. 2007).

Convective mixing was accounted for by using mixing-length theory (Böhm-Vitense 1958) with a mixing length parameter of α = 1.5 commonly used in stellar evolution models (e.g., Pols et al. 1998). We employed values for step overshooting above convective regions of 0.335 pressure scale heights, based on the calibration of stellar models for early B-type stars against observations valid in a comparable mass range as our models, (10 − 20 M_⊙; Brott et al. 2011). We also took into account rotational mixing (Paxton et al. 2013), semi-convection, and thermohaline mixing (Kippenhahn et al. 1980) until the end of core helium burning, even though earlier studies have shown that these have little impact on the stellar structures of stripped stars (Yoon et al. 2010; Götberg et al. 2017).

We used the theoretical wind mass-loss algorithm from Vink et al. (2001) for the main sequence evolution and the de Jager et al. (1988) prescription for stars with hydrogen mass fractions lower than X_H, s = 0.4 and effective temperatures lower than 10⁴ K. Because of the scarcity of observational constraints for the wind mass-loss from stripped stars, we employ the empirically derived wind mass-loss prescription from Nugis & Lamers (2000) based on the winds of Wolf-Rayet (WR) stars. This mass-loss rate is very close to the value derived for the observed intermediate mass stripped star (the qWR star HD 45166, Groh et al. 2008). It is possible that stripped stars have lower wind mass-loss rates than what is predicted from the extrapolated wind mass-loss scheme of Nugis & Lamers (2000) since they are not close to the Eddington limit (Bestenlehner et al. 2014). Indeed, Vink (2017) suggests that the wind mass loss rate from stripped stars is about ten times lower than what is predicted from the Nugis & Lamers (2000) prescription. However, the models from Vink (2017) assume an effective temperature for stripped stars of T_eff = 50 000 K, which is much lower than what is implied by our models for the long-lived phase of central helium burning. If stripped stars have lower wind mass-loss rates than the amount we employ in this study, they may expand more because less of the hydrogen is removed from the surface (cf. Gilkis et al. 2019).

We computed our models by employing the default spatial and temporal resolution of MESA (varcontrol_target = 10⁻⁴ and mesh_delta_coeff = 1.0). We increased the temporal resolution for phases involving the depletion of fuel in the core (e.g., core hydrogen depletion). Due to numerical issues, we lowered the temporal resolution to a maximum of 10⁻³ after core helium depletion and lower the sensitivity of the models to changes in abundances of elements in the core after the formation of an oxygen core.

We computed the interaction with a binary companion using the approach described in Paxton et al. (2015). We took the effect of tides into account using Hut (1981) and use the implicit mass-transfer scheme of Ritter (1988) for the Roche-lobe overflow. As a representative case for stable mass transfer, we treated the secondary component of the binary system as a point mass that has 80% of the primary’s mass. We assumed a conservative mass transfer. We only considered stripped-envelope stars created by the transfer of mass from stars that fill their Roche lobe due to a rapid expansion during their hydrogen shell burning phase after leaving the main sequence (case B mass-transfer, see Kippenhahn & Weigert 1967; Podsiadlowski et al. 1992). To this end, we adopted an initial orbital period of 25 d for all models. Stripped stars created through this channel at solar metallicity have very similar properties regardless of the exact choice of the initial orbital period and companion mass (Götberg et al. 2017). At low metallicity, the efficiency of stripping of the envelope is dependent on the initial orbital period (see Yoon et al. 2017; Ouchi & Maeda 2017). Therefore, at low metallicity, our results should be regarded as a representative approximation.

Stars may, alternatively, be stripped by common-envelope evolution. It is currently not known whether stars which have their envelopes ejected in this manner have different post-envelope-ejection properties than stars which are stripped by stable mass transfer (see, e.g., Ivanova et al. 2013) or even whether there is a metallicity dependence at all.

Stripped stars have relatively small radii during central helium burning (around 1 R_⊙Habets 1986a; Götberg et al. 2017; Yoon et al. 2010) but they are expected to swell up once helium is depleted in the core. In this paper, we investigate the impact of the radius expansion at the end of core helium burning. If models expand enough to fill their Roche lobe anew, they are expected to start an additional mass-transfer phase (sometimes referred to as case BB and BC mass-transfer, see e.g., Dewi et al. 2002). The size of the Roche lobe varies depending on the size of the orbit, which, in turn depends on how angular momentum and mass is transferred. This is still a major uncertainty in binary evolution (e.g., de Mink et al. 2007). To avoid these complications and derive results that are of more generic use, we followed the late evolution of stripped stars by letting them expand as much as their internal structure dictates, ignoring any limitation imposed by the finite but highly uncertain size of the Roche lobe at this stage. In other words, we effectively treat the stripped stars as single stars during their late evolutionary phases. This allows us to investigate their full expansion and simplifies the interpretation of the physical processes involved in the radial expansion. It also makes the results of our simulations more suitable for inclusion in future populations synthesis simulations.

For the analysis, we used the following software: mesaPlot (Farmer 2019), matplotlib (Hunter 2007), numpy (van der Walt et al. 2011), ipython/jupyter (Perez & Granger 2007), and ColorPy¹.

3. Evolution of two representative models

Before presenting our full grid of models, we first describe the evolution and properties of two representative stellar models with identical initial component masses and orbital periods at solar (Z  =  0.0142) and low metallicity (Z  =  0.001). The initial mass of the hydrogen-rich primary star in these models is 10.5 M_⊙, corresponding to the mean initial mass in our grid. Following the method described in the previous section, we place this star in an orbit with an initial period of 25 d, with an 8.4 M_⊙ companion. For the purposes of our investigation, the companion is removed after the mass-transfer phase, when the primary star that has become a stripped star has reached core helium depletion.

Figure 1 compares the evolution of these primary-star models (left panels: solar metallicity, right panels: low metallicity). The top pair of panels presents their evolution in Hertzsprung-Russell diagrams (HRDs). The middle panels present their Lagrangian internal structures as a function of time, known as Kippenhahn diagrams. The lower sets of plots give the time evolution of the stellar radii, along with the fractional contributions of different forms of nuclear luminosity. Key points are labelled with letters in each panel, and is discussed in the following subsections.

Fig. 1. Evolution of a high (left) and low (right) metallicity star of initially 10.5 M⊙ that is stripped due to binary interaction. Letters A to J mark evolutionary points discussed in Sect. 3. Top: evolutionary tracks on the Hertzsprung-Russell diagram. The color indicates the total mass of hydrogen present at a given moment. For comparison, we show lines of constant radii as well as the evolution of the alternate-metallicity model in gray. Middle: evolution of the stellar structure of the stars in mass coordinate (Kippenhahn diagrams), given as a function of time until core collapse (τcc). Convective and overshooting regions are marked with double- and single-hatched regions, respectively. Colors indicate zones dominated by nuclear burning (yellow) or neutrino cooling (purple) where ϵ = sign(ϵnuc−ϵν)log10(max(1.0,|ϵnuc−ϵν|)/[erg g−1 s−1]). Here ϵnuc is the nuclear energy generation rate and ϵν the neutrino energy. The black line indicates the location of the stellar surface. Bottom: luminosity produced by hydrogen-, helium- and metal-burning as a fraction of the total luminosity produced by nuclear reactions (upper) and stellar radius (lower) as a function of time until core-collapse.

3.1. Evolution until core helium burning (A–F)

The early evolution of the primary star in the binary is effectively the same as that of a single star. While we include the effects of rotation and tides, they have a negligible impact on the evolution. During their main-sequence evolution (labeled as A–B in Fig. 1) the stars burn hydrogen into helium in a convective core. As they do so, they increase in luminosity. The higher-metallicity star is cooler and larger due to a higher opacity in the outer layers (due to increased bound-free and bound-bound absorption, e.g., Maeder 1990; Schaller et al. 1992). After leaving the main sequence (B), the stars contract (see also the middle and lower panel of Fig. 1) until hydrogen is ignited in a shell (C). They then expand on their thermal time-scale until they fill their Roche lobe, which leads to transfer of matter to their companion.

During mass transfer (D–E), the stars lose the majority of their hydrogen-rich envelopes on a thermal time-scale (see, e.g., Kippenhahn & Weigert 1967, for more details). The mass-transfer phase occurs at similar effective temperatures in both models presented here since the size of the Roche lobe is the same, although the lower metallicity model has expanded slightly more relative to the higher metallicity model as of the end of the main sequence. More hydrogen is retained after the end of the Roche lobe overflow phase (B–F) in the lower metallicity model (0.27 M_⊙, compared to 0.12 M_⊙ at solar metallicity; see also the upper panels of Fig. 1). This is because the stripping process becomes less effective at low metallicity due to the lower opacity in the outermost layers (see, especially Götberg et al. 2017 and references therein; also Klencki & Nelemans 2018). With a lower opacity, more of the hydrogen-rich envelope remains within the Roche lobe once the stars detach, leaving a thicker layer of hydrogen on top of the metal-poor stripped star. This is clearest in the Kippenhahn diagrams in Fig. 1, where the lower metallicity model has more mass outside the hydrogen-burning shell. After the end of the mass-transfer phase the stripped star shrinks towards a new gravothermal equilibrium structure. Meanwhile, convective helium burning has already started in the core. From the luminosity minimum (F), the dominant driver of structural change is, once more, nuclear burning.

3.2. Core helium burning phase (F–G)

Point F in Fig. 1 marks the longest-lived phase of these stripped stars (about 10% of the stellar lifetime), with helium burning in a convective core. The metal-rich model is hotter and more compact than the metal-poor model because the metal-poor model retains more hydrogen in its envelope, which allows it to be larger (Cox & Salpeter 1961). From here onwards, the evolution varies considerably between the two models. The metal-rich star shrinks monotonically during the core helium burning phase, while the metal-poor model first shrinks, then stays approximately constant, before starting to expand again (e.g., Götberg et al. 2017).

3.3. Expansion phase after core helium depletion (G–J)

After core helium depletion, the stars begin the short last phase of their lives, which lasts for less than 1% of their total stellar lifetime. Points labelled G in Fig. 1 represent the moment at which the central helium mass fraction drops below 10⁻⁴. At this point, the metal-rich model has a total mass of 2.98 M_⊙ and has lost all its remaining hydrogen due to stellar winds, whereas the metal-poor model has a total mass of 3.90 M_⊙ and a total mass of hydrogen of 0.14 M_⊙. The structures during the post-core-helium-burning phase are presented in more detail in Fig. 2, using both mass and radius coordinate systems (the top and bottom panels, respectively).

Fig. 2. Kippenhahn diagrams showing the evolution of the stellar structure for our two example models at solar (left panels) and low (right panels) metallicity. They are shown using mass coordinate as the vertical axis (top panels) and radius coordinate (bottom panels). The horizontal axis indicates the time since helium depletion up to core carbon depletion. The black line represents the surface of the star. Regions of mixing by convection and overshooting are marked with double- and single-hatches, respectively. Colors indicate zones dominated by nuclear burning or neutrino cooling. See Fig. 1 and Sect. 3 for details and a discussion.

We see large differences between the high- and low-metallicity models in this last evolutionary phase. The radius evolution becomes very distinct, with a monotonic increase of the radius for the metal-rich model until a maximum of 15 R_⊙ is reached, and a non-monotonic increase of the radius for the metal-poor model until a maximum of 530 R_⊙ (see the upper and lower panels of Fig. 1). We further find that while the burning processes in the inner core are fairly similar, the metal-poor model develops prominent convection zones in the outer envelope and around the helium burning shell, and still experiences hydrogen shell-burning, while the metal-rich model has lost all its remaining hydrogen layer by this point due to winds (see Fig. 2 and the middle panels of Fig. 1).

Some general features are common to the evolution at both metallicities, for example, the narrowing of the helium-burning shell (in mass extent, see upper panels of Fig. 2) as the size of the carbon-oxygen (CO) core decreases (in radial extent, see lower panels of Fig. 2), and the increasing rate of neutrino cooling from the CO core as it contracts. During this phase of core contraction, both models expand overall. Since the mass of the CO core is somewhat higher in the less-stripped low-metallicity model, the evolution occurs more rapidly here.

Figure 2 shows that the low-metallicity model also retains a significant hydrogen layer which can sustain a burning shell, unlike the high-metallicity model that loses its hydrogen-rich envelope due to stellar winds. For the low-metallicity model hydrogen shell burning dominates the nuclear luminosity around the time of core helium exhaustion, and there is a local maximum in the radius evolution (G′), as shown in the lower-right panels of Figs. 1 and 2. In contrast, by point G, hydrogen burning is not relevant for the higher metallicity model and soon a shallow minimum in the radius evolution is apparent.

The radius expansion during the shell-burning phases, including the non-monotonic expansion of the low-metallicity model, might be interpreted in terms of the “mirror principle” (see, e.g., Kippenhahn et al. 2012). In examining the middle and lower panels of Fig. 1 and 2, respectively, the phases of radius expansion occur when one shell-burning source dominates the nuclear luminosity. The lower-right panels of Fig. 1 show that the first peak in the radius expansion of the low-metallicity model occurs when two shell sources are releasing roughly equivalent luminosities, analogous to the blue loop observed in some models for intermediate mass stars (see e.g., Kippenhahn et al. 2012). The layer above the helium-burning shell expands, as seen in Fig. 2, which leads to the cooling of the hydrogen burning shell. Eventually the temperature and density of the hydrogen-rich material are too low to sustain hydrogen burning and the shell is extinguished.

Carbon is later ignited in the center, leading to the development of a new convective core (point H indicates this, showing the moment when the core carbon abundance drops by 2% from the post-core-helium-burning value). This is associated with a sharp rise in the burning luminosity of elements heavier than helium, as seen in the lower panels of Fig. 1. By this time, both stars, of 2.98 M_⊙ and 3.90 M_⊙ for the metal-rich and the metal-poor model, respectively, have only one shell source, helium. The hydrogen-rich envelope of the low-metallicity model still means that the star is more than an order of magnitude larger than the solar-metallicity model.

At the end of the evolution, more than 90% of the envelope has become convective and extends over 400 R_⊙, while only containing a few tenths of a solar mass of material. In reality, we do not expect the stripped star to be able to expand this much. It will likely interact again with its companion, which is still expected to be present. This would lead to another phase of mass transfer, with potentially important consequences for the type of supernova (see, e.g. Yoon et al. 2010, 2017; Dessart et al. 2018). After core carbon depletion is reached (as marked with letter J), less than 100 yr remains before these models reach core collapse. Since the thermal timescale of the stars are much longer than 100 yrs, their radii, masses, and surface composition are not significantly affected after central carbon depletion. We can, therefore, regard the properties the stripped stars have at point J as the properties the stripped stars exhibit at the moment of explosion.

4. Comparing the metal-poor and metal-rich grids

Here we present results from two grids of evolutionary models of stars that are stripped through binary mass transfer. The grids are computed at two different metallicities (solar and Z = 0.001) and each consists of twenty-three models with different initial masses for the primary star, ranging from 8.8 to 15 M_⊙. Tables A.2 and A.3 provide an overview of the key parameters.

4.1. Evolutionary tracks in the Hertzsprung-Russell diagram

The evolutionary tracks are presented in Fig. 3, with the top panel showing the results for solar metallicity and the bottom panel the results for low metallicity. The tracks show the evolution from the onset of hydrogen burning up to the completion of helium burning in light grey; for a detailed description, see Sects. 3.1 and 3.2. Here we highlight in color the last phases of the evolution from the completion of central helium burning until central carbon depletion, where the color indicates the initial mass of the progenitor.

Fig. 3. Evolutionary tracks of stars at solar (top) and sub-solar metallicity (bottom) showing their luminosity, L as a function of effective temperature, Teff. Lines of constant radius are provided for reference. We highlight the evolution after core helium depletion in color, where the color refers to the initial mass. Yellow star symbols indicate core carbon exhaustion, and white star symbols denote the final time step calculated for models that did not reach core carbon exhaustion. The white circles in the lower panel indicate when the radius reaches a first maximum after core helium depletion. Above each panel, we indicate approximate ranges of effective temperatures typical for Wolf Rayet stars (WR), blue, yellow, and red supergiants (BSG, YSG and RSG respectively). The horizontal color scale above each panel indicates, for each effective temperature, the color as perceived by the human eye.

The most striking feature in these plots is how far the evolutionary tracks extend to the right during the late stages of their evolution, that is, how much their effective temperature decreases and their radius expands before carbon is depleted in the core (which is marked with a yellow star symbol). In our solar-metallicity models, shown in the top panel, there is a wide range of final effective temperatures and thus final radii, which is in general agreement with earlier studies (cf. Habets 1986b; Dewi & Pols 2003; Yoon et al. 2010, 2017). The lowest mass models in our grid (M_ZAMS <  9.5 M_⊙) reach the lowest final effective temperatures and largest final radii of log₁₀T_eff/K ≈ 3.8 and R ≈ 150 R_⊙, respectively. These temperatures are typical for yellow supergiants (YSG, e.g., Drout et al. 2009). For the intermediate-mass models in our grid, we find final temperatures and radii that are more typical for blue supergiants (BSG, e.g., Fitzpatrick 1988). The highest mass models in the grid reach final effective temperatures and radii of up to log₁₀T_eff/K ≈ 4.7 and as low as R ≈ 3.5 R_⊙, respectively. Their wind mass-loss rate is of the order of 10⁻⁶ M_⊙ yr⁻¹ or more. These properties are characteristic for classic Wolf-Rayet stars (WR, e.g., Crowther 2007). All the models in our low-metallicity grid end their lives as cool (log₁₀T_eff/K ≊ 3.55) and large (R >  400 − 700 R_⊙) stars that are typical for red giants or red supergiants (RSG, e.g., Groh et al. 2013).

4.2. Evolution of the radius and its connection to the nuclear burning shells

The evolution of the radii is shown in the panels in Fig. 4 as a function of time since helium depletion. The solar-metallicity models are shown on the left and low-metallicity models on the right. The color bars indicate the relevant physical quantities, which are different for each row; this is discussed in more detail further on in this paper. Each panel show a sequence of twenty-three tracks which correspond to our models for different initial masses. The tracks for the massive models can be readily identified as they evolve faster and complete their final phases of evolution in about 20 kyr. The lower mass models in our grid take 125 kyr at solar metallicity and about 80 kyr at low metallicity.

Fig. 4. Evolution of the radius as a function of time after core helium depletion for the stars in our model grids at solar (left) and lower (right) metallicity. From top to bottom, colors indicate the remaining mass of hydrogen and the normalized luminosity from nuclear burning of hydrogen, helium, and heavier elements; see text for details. More massive stars evolve faster and are thus located towards the left of the figures. White crosses indicate the moment of core carbon ignition. Yellow star symbols indicate core carbon exhaustion, and white star symbols denote the final time step calculated for models that did not reach core carbon exhaustion.

We see that the solar-metallicity models expand from about 0.5 R_⊙ up to final radii of 3 R_⊙ for the more massive models, while the lower mass models expand to radii of 180 R_⊙. In contrast, the low-metallicity models start at radii of 3 to 5 R_⊙ and reach their final radii of 400–700 R_⊙. For an overview, see also Tables A.2 and A.3.

For the low-metallicity models, the post-core-helium-burning radius expansion is clearly not monotonic with time, as we discuss in Sect. 3. The stars expand to a first maximum around 8 kyr after helium depletion for all models, followed by a contraction and re-expansion phase. The first expansion phase is most significant for the most massive models, which reach radii of up to 300 R_⊙. The total radius expansion for the low-metallicity models is very significant (two orders of magnitude). The implication of this is that these stars are expected to exceed the size of their Roche lobe and, at low metallicity, they may already do so during the first expansion phase. This would trigger at least one additional phase of mass transfer shortly after core helium depletion; those are not modeled here. Moreover, they may well fill their Roche lobe shortly before or at the moment of explosion.

To understand the origin of the difference in behavior between solar- and low-metallicity models, it is helpful to inspect subsequent physical parameters. The main difference between the solar- and low-metallicity models is the remaining hydrogen mass, which is indicated in color in the panels in the top row of Fig. 4. The solar-metallicity models contain almost no hydrogen, that is, less than about 0.05 M_⊙ after the mass transfer ceases and even less at the end their evolution (see $M_{\mathrm{H}}^{\mathrm{postMT}}$ and $M_{\mathrm{H}}^{\mathrm{f}}$ in Table A.2). In contrast, the low-metallicity models contain about 0.2 M_⊙ of hydrogen after the mass transfer ceases and the stellar winds are too weak to substantially reduce this afterwards.

We can gain further insights in considering the nuclear burning processes that contribute to the total luminosity by nuclear burning (L_nuc). Colors in the lower panels of rows 2, 3, and 4 of Fig. 4 show the relative contribution to the total nuclear luminosity resulting from the burning of hydrogen (L_H), burning of helium through the triple-alpha reaction (L_He), and the collective burning of heavier elements (L_Z), including helium burning by alpha captures onto carbon.

During the first 10 000 years after helium depletion, the radii of the solar-metallicity models change slowly, while the low-metallicity models expand rapidly. At solar metallicity, the thin hydrogen-rich layer that is left after Roche-lobe stripping, if any, is not sufficient to support hydrogen shell burning. By contrast, in the low-metallicity models enough hydrogen is retained to sustain hydrogen shell burning. Hydrogen burning dominates the nuclear luminosity during the first expansion phase, as can be seen on the second row of Fig. 4. At the peak of the first stellar expansion hydrogen burning contributes about half of the total nuclear luminosity. At that time, the contribution from the helium-burning shell is increasing, as can be seen in the panel on the third row of Fig. 4. The turning point in radial expansion occurs when the stars have roughly equivalent luminosity contributions from two shell sources. The first expansion in the low-metallicity models is thus associated with hydrogen shell-burning, and the subsequent contraction is consistent with the “double mirror effect” (see also Sect. 3.3).

Thereafter, both high- and low-metallicity models show radius expansion when helium-shell burning dominates the nuclear luminosity. The low-metallicity models reach much larger radii, and most of their expansion is during this burning phase. The plus symbols in Fig. 4 indicate the start of core carbon burning, which we define as the moment when the core carbon abundance drops by 2% below its post-core-helium-burning value. The solar-metallicity models show significant further expansion after this time, for the lowest mass models by more than an order of magnitude in radius, as a result of to the left-over hydrogen layer. The low-metallicity models have completed most of their expansion by the onset of core carbon burning.

4.3. Binding energy

In the previous sections, we discuss the expansion of stripped stars in their final evolutionary phases. As a result of this expansion, stripped stars can fill their Roche lobe again, initiating a subsequent phase of mass transfer. In the case of unstable mass transfer, the system is expected to enter a common-envelope phase. This would shrink the orbit and, thus, carry significant potential consequences for the final fate of the binary, possibly as source of gravitational waves. We briefly discuss the most important stellar property affecting the outcome of a common envelope phase that we can compute, namely, the binding energy of the envelope (e.g., Webbink 1984).

We define the binding energy E_B of the (hydrogen or helium) envelope as:

$$\begin{array}{c}\hfill E_{\mathrm{B}}=-{\displaystyle {\int }_{\mathrm{core}}^{\mathrm{surface}}(-\frac{\mathit{Gm}}{r}+ϵ(m))\mathrm{d}m,}\end{array}$$(1)

where G is the gravitational constant, m the mass coordinate, and r the radius coordinate of the star, and ϵ the specific internal energy. This internal energy includes not only the thermal energy terms, but also the potential energy stored in ionised species and dissociated molecules (see, e.g., Han et al. 1994; Ivanova et al. 2013). We define the relevant core-envelope boundary as the point at which the abundance of the dominant element in the envelope drops below 10%. The chosen definition is somewhat arbitrary, but we consider it reasonable for our current purposes (see also Appendix E).

The results for our solar-metallicity models are shown in the top panel of Fig. 5. The binding energy of the H envelope is negligible, as expected, since most hydrogen has been removed efficiently during Roche-lobe overflow and subsequently by the stellar wind. Shortly after helium depletion, the He+H envelope has a binding energy of the order of about −5 to −10 × 10⁴⁹ erg, where the more massive stars in our grid are more tightly bound (i.e., their binding energies are more negative). As time proceeds, the stars evolve and expand. Their envelopes loosen to reach values of the binding energy of about −2 to −6 × 10⁴⁹ erg at the end of their evolution, when the envelope starts to become convective.

Fig. 5. Binding energy of the envelope as a function of radius after core helium depletion at solar (top panel) and low (bottom panel) metallicity. We show values for the hydrogen-rich envelope alone, as well as the hydrogen and helium envelope combined, as labelled. Colors indicate the mass of the stripped star.

For our low-metallicity models, shown in the bottom panel, we find the same qualitative trend with mass and time. However, the binding energies are about −1 × 10⁴⁹ erg for the H envelope and about ten times more than that for the H+He envelope. All these models develop a convective envelope when their radius reaches 200–300 R_⊙. We discuss the consequences for the stability of mass-transfer and gravitational-wave progenitors in Sect. 6.2.

5. Comparison with radii for helium stars adopted in population synthesis simulations

Having accurate estimates for the minimum and maximum radii of stars is crucial in the context of binary systems. For an individual binary system, the radial expansion determines whether a star can swell to fill its Roche lobe and start a (new) phase of mass transfer. For a full population of binary systems, the radial expansion is an essential factor in determining what fraction of this population interacts.

To simulate large populations of binary systems and their exotic end products as X-ray binaries and gravitational wave progenitors, prescriptions are needed that can be evaluated quickly. Generally, two approaches are commonly adopted. Many population synthesis codes make use of the analytic formulae by Hurley et al. (2000). We compare our results with these prescriptions in Sect. 5.1. Alternatively, one can make use of interpolation in grids of pre-computed stellar evolution models. One example of this is the SEVN code, with which we compare in Sect. 5.2.

5.1. Analytic prescriptions by Hurley et al. (2000)

The Hurley fitting formulae are adopted in many population synthesis codes that are in active use. These include, but are not limited to, BSE (Hurley et al. 2002), StarTrack (Belczynski et al. 2008), Binary_c (Izzard et al. 2004, 2006), COMPAS (Stevenson et al. 2017; Vigna-Gómez et al. 2018), and the latest arrival COSMIC (Breivik et al. 2019). The discussion below thus applies to these studies and others that make use of these prescriptions.

The Hurley fitting formulae are based on evolutionary models by Pols et al. (1998) which solve for the full set of stellar structure equations using the Eggleton code (Eggleton 1971, 1972; Eggleton et al. 1973) with updates to the equation of state by Pols et al. (1995). Specifically, the fitting formulae for helium stars (provided in Sect. 6.1 of Hurley et al. 2002) are based on detailed models for single stars computed with the same input physics as the grid presented in Pols et al. (1998). These models assume a homogeneous initial composition with a helium mass fraction of Y = 0.98, a mass fraction of heavier elements of Z = 0.02 and no hydrogen, X = 0 (see Dewi et al. 2002; Dewi & Pols 2003, for a discussion).

The main disparity from our models is that we self-consistently account for the stripping process due to a binary companion, assuming a representative initial orbital period. In our models, the resulting stripped star can still contain a remaining layer of hydrogen at the surface. The effect of such a remaining layer of hydrogen is thus not accounted for in the Hurley prescriptions. A further difference is that the original stellar models behind the Hurley prescriptions were provided for only one fixed value of the metallicity. In our calculations we find large differences between stripped stars at solar and at low metallicity, which are not accounted for in the Hurley prescriptions.

The top row of Fig. 6 compares our results (in blue) with the Hurley prescriptions (in grey). For solar metallicity, shown in the top left panel, we find that minimum radii given by the Hurley prescription are comparable to the minimum radii we find in our detailed models. The maximum radii from Hurley et al. (2000) show a similar trend with mass as we find in our models, but we also observe significant and meaningful differences, in particular at the low-mass end (helium core masses of 2–2.5 M_⊙). Our models reach maximum radii of about 200 R_⊙, an order of magnitude larger than the maximum radii for the Hurley prescriptions, which reach only about 20 R_⊙.

Fig. 6. Comparison of our results (shown in blue in each panel) with (1) the widely-used prescriptions by Hurley et al. (2000; shown in grey in the panels on the top row) and (2) the PARSEC models (Bressan et al. 2012) used in the SEVN code Spera et al. (2019; shown in grey in the panels on the bottom row). Panels on the left: comparison at solar metallicity (or its closest available equivalent, i.e. Z = 0.02 for SEVN) and panels on the right: comparison for low metallicity (Z = 0.0002 for SEVN). We show the maximum radius (circles) and the minimum radius (triangles). For our models, we also mark the radius at core carbon ignition (squares).

At low metallicity, shown in the top right panel of Fig. 6, we find that the Hurley prescription under-estimates the minimum radii by about a factor of three. The differences are much larger for the maximum radii. The maximum radii in our models increase with increasing mass and reach radii of 400–600 R_⊙, which is two to three orders of magnitude larger than the maximum radii from the Hurley prescriptions. We find the largest differences in radii at the high-mass end, where our models have radii that are 500 times larger than the Hurley maximum radii. These differences are due to the presence of the remaining hydrogen layer in our stripped models, which is not accounted for in the Hurley prescriptions.

5.2. Example of grid based interpolations: SEVN code

An example of a code that uses interpolation within a pre-computed stellar model grid is SEVN, a grid-based population synthesis code (Spera et al. 2019) which has also been used to make predictions for gravitational wave progenitors. The stellar models behind this code are single stellar models computed with the PARSEC stellar evolution code (Bressan et al. 2012).

The bottom panels of Fig. 6 compare our results to those of SEVN that are closest in metallicity. In the lower left panel, we compare our solar-metallicity models (Z = 0.0142) to their models (those which assume a slightly higher value for the solar metallicity, Z = 0.02, which is the old canonical value for solar metallicity). In the lower right panel, we show the PARSEC models for Z = 0.0002, which is closest to the metallicity of Z = 0.001 that we adopted for our low-metallicity grid.

We find similar general differences as described above when making comparisons with the Hurley prescriptions. For solar metallicity, we find, again, that the minimum radii are in fairly good agreement. For the maximum radii, again we note significant differences depending on the mass. We find significantly larger maximum radii of about 200 R_⊙, compared to about 60 R_⊙, for the PARSEC models at the lower mass end (for masses below about 2.7 M_⊙). For larger masses, the models used in SEVN reach maximum radii of 10 to 50  R_⊙, which is substantially larger than the maximum radii we find (of 5 to 30 R_⊙) in our models. Understanding the origin of these differences would require further investigation. It may be due in part to the difference in metallicity, but differences in the micro-physics or treatment of convection may also play a role. At low metallicity, we again find very large differences, similar to – but even more pronounced – than the differences we find with the Hurley prescription, as shown in the lower right panel of Fig. 6. We provide new analytic fits to our models for use in population synthesis calculations in Appendix D.

6. Discussion

As shown in Sect. 5, we find systematically larger radii for stripped stars than those commonly used in population synthesis calculations. The large radii can trigger additional phases of mass and angular momentum transfer (traditionally referred to as Case BB or Case BC mass transfer, Dewi et al. 2002, and references therein).

Such additional interaction can impact the final masses and orbital separation and they are, thus, important for modeling the populations of binaries. Specifically, these later phases of interaction are thought to be very important in the formation of peculiar supernovae and gravitational wave progenitors (Ivanova et al. 2003; Dewi & Pols 2003; Tauris et al. 2013). Moreover, Zevin et al. (2019) argue that these additional mass-transfer phases are necessary to explain enrichment in globular clusters, assuming that r-process enrichment primarily originates from double-neutron star systems.

A full assessment of the implications would require extended grids of models that follow these additional phases of mass transfer self-consistently. We provide a first estimate of the additional number of systems affected compared to the widely used prescriptions in Sect. 6.1. In Sect. 6.2, we discuss the question of whether the late phase of mass transfer would be stable or would lead to a common envelope phase involving a neutron star. We examine the implications for the observability of these stars in Sect. 6.3 and for supernova progenitors in Sect. 6.4. In Sect. 6.5 we discuss the main uncertainties that affect the results presented in this work.

6.1. Expected increase in the number of binary systems that interact with a helium donor

Whether or not the stripped star will fill its Roche lobe anew depends on the size of the Roche lobe, which scales linearly with the separation a for a given mass ratio (e.g., Eggleton 1983). The separation of a particular binary system depends on its initial separation, the amount of mass that is transferred, and the amount of angular momentum that is lost from the system during the first mass transfer. Given the large uncertainties in the mass-transfer process, the distribution of separations is not well known.

In order to make a simple estimate, we make the agnostic assumption that the separations (or, more precisely, the Roche-lobe radii) are distributed according to the standard Öpik (1924) law, which is a distribution that is flat in the logarithm. We further assume that the separations span the full “range of interest”, meaning they are such that the Roche radii span from R_min to at least R_max, where,

$$R_{min}\equiv min(R_{min}^{\mathrm{this}\mathrm{study}},R_{min}^{\mathrm{Hurley}})\text{and}{R}_{max}\equiv max(R_{max}^{\mathrm{this}\mathrm{study}},R_{max}^{\mathrm{Hurley}}).$$

Here we use the superscripts to indicate the origin of the minimum and maximum radii.

Assuming Öpik’s law leads to the following expression for the relative number of systems that interact with our new estimates for the radial expansion, compared to what would have been obtained with the Hurley prescriptions:

$$\begin{array}{c}\hfill f\approx \frac{log{R}_{max}^{\mathrm{this}\mathrm{study}}-log{R}_{min}}{log{R}_{max}^{\mathrm{Hurley}}-log{R}_{min}}·\end{array}$$(2)

With this, we find that a stripped star is about twice as likely to interact relative to Hurley et al. (2000) if we consider solar metallicity progenitors in the stripped-star mass range 2 − 2.5 M_⊙. At low metallicity, this fraction grows even larger due to the larger increases in radius. With this simple estimate, we find that stripped stars between 2 and 6 M_⊙ are 2–30 times more likely to interact relative to the Hurley et al. (2000) prescriptions.

The numbers quoted here should be taken with a grain of salt. We have little reason to expect a logarithmically flat separation distribution to be realistic for systems that have already gone through a phase of interaction. Moreover, this simple estimate does not take into account the dependence of the low-metallicity models on the initial orbital parameters. For short enough orbits, stars at lower metallicity would lose their hydrogen-rich envelopes after the first binary interaction, leading to a smaller increase in radius (see, e.g., low-metallicity models with short orbital periods of Yoon et al. 2017). However, it is noteworthy that at low metallicity, 100% of all stripped stars computed with these initial orbital parameters will fill their Roche lobe anew (see Fig. 3). This is in stark contrast with what follows from the Hurley prescription that predicts no stripped stars in this mass range to fill their Roche lobe again.

6.2. Unstable mass transfer in systems with neutron star companions

The large radii of stripped stars may allow them to fill their Roche lobe anew and start to transfer mass to their companion. The case where the companion is already a neutron star is of particular interest since such a system is a possible immediate progenitor of a double neutron star system (e.g., Fragos et al. 2019).

When the stripped star fills its Roche lobe, its mass is still expected to exceed that of a typical neutron star of 1.4 M_⊙. If the mass transfer is stable, the orbit is expected to shrink because mass is transferred from a more massive star to a less massive companion (Paczynski 1971). Secondly, mass lost from the system is likely to be emitted primarily from the vicinity of the less massive neutron star. The mass lost thus likely carries a specific angular momentum that is similar to or larger than the specific orbital angular momentum of the neutron star. This is larger than the average specific orbital angular momentum and, thus, we expect the orbital separation to shrink (van den Heuvel et al. 2017).

A more dramatic shrinking of the orbit is expected when mass transfer is unstable. In this case, the neutron star becomes engulfed in the envelope of the donor (known as common envelope (CE) evolution, for a review see Ivanova et al. 2013) and the orbital separation can be shortened drastically, depending on the binding energy of the envelope and the efficiency with which it is ejected.

To know whether or not a mass transfer is unstable would require further detailed calculations. For a first estimate, we assume that unstable mass transfer occurs in these systems only if the donor star has a convective envelope. The stripped stars we consider have masses in the range of 2–6 M_⊙, hence, if the companion is a 1.4 M_⊙ neutron star, only the highest-mass of these stripped stars would canonically be expected to undergo unstable mass transfer when they have radiative envelopes.

In Fig. 7, we show various tracks for our stripped stars where we highlight systems with convective envelopes at the onset of mass-transfer with colors. Stripped stars with radii larger than about 200 R_⊙ have convective envelopes. Assuming again that these systems are distributed flat in log a implies that about a fifth of the systems would begin mass transfer with a convective envelope. However, as a caveat, we note that these stars do not develop very massive convective envelopes (see for example Fig. 2). We also note that if the remaining nuclear-burning lifetime is very short; Tauris et al. (2015) argue that there may not be sufficient time to complete the common-envelope inspiral before core collapse.

Fig. 7. Top: binding energy of the envelope plotted as a function of radius. On the top x-axis we have indicated the orbital separation assuming a neutron star companion with a canonical mass of 1.4 M⊙, for reference. Our solar models are shown in the top left panel where we show the binding energy of the helium envelope. Our low-metallicity models are shown in the top-right panel, where we show the binding energy of the hydrogen envelope (cf. Sect. 3.3 and Fig. 5). Bottom: final separation after the common envelope evolution of stripped stars with neutron star companions as a function of the initial separation at solar (left) and low (right) metallicity. Stars with a pre-CE convective envelope are marked with colors, which indicate the mass of the stripped star, assuming a common envelope efficiency parameter αCE = 1.

The top panels of Fig. 7 show the binding energy as a function of radius and, for reference, the corresponding separation if the companion were a neutron star. In the bottom panels of Fig. 7, we present the predicted post-CE separation as a function of the pre-CE separation if the mass transfer was to be unstable. For this calculation, we take the standard assumption that the orbital energy is completely converted into the binding energy of the envelope (common envelope efficiency parameter α_CE = 1, Webbink 1984) and that the envelope is ejected with exactly the necessary escape velocity. We compute the envelope structure parameter λ calculations in the Appendix E.

At low metallicity, we find that the final separation is between 2 and 3 R_⊙ for initial separations smaller than 300 R_⊙ through stable mass-transfer. This is of potential interest for studies of gravitational wave progenitors. Above these initial separations, the envelope becomes convective and the final separations reach values between 3 and 100 R_⊙ for the lowest mass models. These systems could become gravitational-wave sources if the final explosion marking the formation of the second neutron star results in a tighter orbit.

Our detailed simulations show strong metallicity effects which might affect the formation of double neutron stars. These effects are currently not included in the vast majority of binary population synthesis codes and may alter the rates and distributions of compact object mergers. Population synthesis simulations typically predict that double-neutron star merger rates are only very weakly metallicity-dependent (e.g., Neijssel et al. 2019).

6.3. Observability

Stars stripped in binaries are notoriously difficult to detect during their longest-lived phase of core helium burning with current instruments. Not only are they compact, with typical sizes of about 1 R_⊙, but most of their radiation is emitted in the extreme ultraviolet. Their companions typically outshine them in optical wavelengths (Götberg et al. 2017, 2018). However, this is no longer the case for the later evolutionary phases on which this work focuses. Their radius and luminosity increase until they reach giant sizes, a finding which is in agreement with previous studies (e.g., Habets 1986a; Yoon et al. 2010, 2012). Schootemeijer et al. (2018) discuss several systems which probably contain a helium-shell burning stripped star. Of those, the closest to the models we present is υ Sgr, with an inferred stripped-star mass of 2.5 M_⊙ (Dudley & Jeffery 1990).

All our model stars reach effective temperatures ranging from 4000 to 10 000 K, which spans typical ranges for WR stars and YSGs, while reaching the high end of the temperature range for RSGs (see Fig. 3). We can estimate their observational characteristics based on their composition, luminosity, and effective temperature. At high metallicity, typical spectra should be helium and nitrogen-rich, and hydrogen-poor. We expect surface gravity values ranging from 0.3 < log₁₀(g/[cm s⁻²]) < 4 (see Table A.2). As helium lines can be challenging to measure, these spectra could be confused with those from nitrogen-enhanced B-stars.

At low metallicity, typical spectra should be similar, but the amount of hydrogen varies depending on the mass-transfer history and the orbital separation. We find extremely low surface gravity values of log₁₀(g/[cm s⁻²]) ≊ −0.4 (see Table A.3) and expect narrow absorption features in their spectra, provided they do not experience strong mass-loss episodes. Given the lack of resolved stellar populations at such low metallicity, it may be challenging to measure spectra of systems containing such giant stripped stars with the present facilities. The Magellanic Clouds may provide examples that are sufficiently varied in composition from those of the Galaxy so that we may test the metallicity trends we describe here.

Stripped stars with the lowest effective temperatures and highest luminosities are expected to be easiest to detect at optical wavelengths. These may well dominate the total emission of their binary system. These may be observed as helium red giants (Trimble & Paczynski 1973; Yoon et al. 2012). The most promising way to detect these systems is by detecting sources that are overluminous for their (Keplerian) masses, just like the system ϕ Persei, discussed by Schootemeijer et al. (2018). These may be identified by searching for discrepancies between spectroscopic, evolutionary, and Keplerian masses, if available.

6.4. Supernova progenitors

Our findings are relevant for understanding the properties of core-collapse supernovae, in particular stripped-envelope supernovae of type Ibc and type IIb (for a pioneering work on the subject, see Podsiadlowski et al. 1992; for interpretations of the class-defining Type IIb SN 1993J, see Podsiadlowski et al. 1993; Nomoto et al. 1993; for recent studies see, e.g., Bersten et al. 2012; Eldridge et al. 2015; Dessart et al. 2018; Sravan et al. 2019).

At the moment of explosion, the supernovae from our progenitor models may be classified as type Ibc or type IIb supernovae, depending on the amount of hydrogen and helium retained (see also Yoon et al. 2017). The large radii of many of our models can be identified through shock cooling signatures in early light-curves (e.g., Schawinski et al. 2008; Yoon et al. 2010; Piro & Nakar 2013). Next we discuss two potential consequences for the evolution of these stars that, at least to our knowledge, have not been previously pointed out.

6.4.1. Circumstellar material prior to explosion

At low metallicity, we find that the stars experience two phases of expansion. The first occurs shortly after helium depletion and is associated with hydrogen shell burning. The expansion is most significant in our models for higher mass stars, which expand by more than an order of magnitude. After this phase, the stars contract until they expand again at the end of their life. The first phase of expansion is so severe for the more massive progenitors that we expect them to briefly fill their Roche lobe. If the resulting mass-transfer event is non-conservative, we would expect an ejection of mass prior to the explosion. For the most massive progenitors, this occurs about 10 000 yr before their terminal explosion. This is short enough that the ejected material may still be close enough to the star at the moment of explosion to interact with the supernova shock. If so, these systems might be progenitors of at least some type Ibn supernovae (see, e.g., Foley et al. 2007; Pastorello et al. 2007, 2008; Hosseinzadeh et al. 2019).

6.4.2. Asymmetric supernova progenitors

Due to the large radial expansion, we expect many stripped stars to fill their Roche lobe again shortly before their death or even at the moment they explode. This implies that the supernova progenitor is not spherical, but instead has a non-axisymmetric, “pear-like” shape as imposed by the shape of the Roche lobe. We expect this to have implications for the explosion, in particular, when the supernova shock reaches the outer layers where the deformation is strongest.

Afsariardchi & Matzner (2018) study explosions in non-spherical progenitors. They find that aspherical progenitors may have different shock-breakout signatures that are viewing angle-dependent. They further find that asphericity leads to collisions that would otherwise not have occurred and affect the observables. The case of a non-axisymmetric progenitor has not been modeled in their study. Although they argue the effects of asphericity in extended progenitors would be weak, they do not model Roche-lobe filling stars at explosion. Given that these configurations may be common, we encourage further detailed studies of these effects to be carried out.

6.5. Uncertainties

The evolutionary models presented in this work are affected by several uncertainties. We have considered a fixed typical initial orbital period and mass ratio. For our high-metallicity models we have verified that the exact choice of the initial orbital period and companion mass has very little effect on the maximum radius of stripped stars if varied within reasonable limits (cf. Götberg et al. 2017, 2018). At low metallicity, Yoon et al. (2010, 2017), Claeys et al. (2011), and Ouchi & Maeda (2017) show that the amount of hydrogen left at the surface of the donor star is a function of the initial orbital separation. In Appendix C, we demonstrate that for case B mass transfer, our general findings are robust against variations of the orbital period at low metallicity.

We further emphasize that we have modeled the late evolution of stripped stars by allowing them to fully expand, since our objective was to determine their maximum radii. In reality, we expect these stars to still be in orbit around a companion star. This would truncate their expansion and initiate a new phase of Roche-lobe overflow.

We consider the effects of internal mixing by convection and overshooting. Internal mixing is one of the main uncertainties in stellar evolution and also affects our results. Yoon et al. (2017) have pointed out the importance of mixing in the region above the retreating convective core of the donor star while it is on the main sequence. The choice of overshooting, efficiency of semi-convection, and other potential mixing processes affects the details of the chemical profile. This, in turn, modifies the response of the donor star to stripping process and determines the mass of the hydrogen-rich layer that is left after the stripping process. Observations of stripped stars, for example, as proposed by Götberg et al. (2018), may help to constrain these uncertainties in the future. The presence or absence of the expansion described in Sect. 6.4, which is equivalent to a blue loop, probably depends on the treatment of internal mixing (although similar additional expansion phases can be also seen in models from Yoon et al. 2017).

When estimating the number of systems that interact anew we assume a uniform distribution of orbital separations in log space. The actual distribution of separations is uncertain and depends on both the distribution of initial separations and the amount of orbital shrinking and widening during the first phase of interaction. A more advanced assumption would be to use the output of population synthesis simulations (Dewi & Pols 2003; Ivanova et al. 2003). We include estimates for the binding energy of the envelope. These depend on the chosen definition of the core-envelope boundary (e.g., Tauris & Dewi 2001) and on the energy terms that are given (e.g., Ivanova et al. 2013). This is discussed further in Appendix E.

Stellar wind mass loss is also an important uncertainty for the evolution of massive stars (Smith 2014; Renzo et al. 2017). In the cases studied here, mass is primarily removed by Roche-lobe overflow and not by the stellar winds, so we do not expect a large impact on our results. However, we stress that the mass-loss rates for stripped stars are uncertain (Götberg et al. 2017, 2018; Vink 2017; Gilkis et al. 2019). This should not have a large effect on our solar-metallicity models, except perhaps for the highest mass models in our grid. The effect of wind mass loss at low metallicity is negligible, as we demonstrate in Appendix B.

The spatial and temporal resolutions, along with the nuclear network that we have adopted, should be sufficient for our purposes. We note that these models have not been optimized as input for supernova simulations, as this requires higher resolution, a more extended nuclear network, and calculations of the final evolutionary steps (e.g., Farmer et al. 2016). However, we will present and discuss such models in a subsequent paper (Laplace et al. in prep.).

7. Summary and conclusion

In this paper, we study the radius evolution of stars that have lost most of their hydrogen-rich envelope due to an interaction with a companion. We consider stars with initial masses of 8–15 M_⊙ at solar and low metallicity. We investigate how the internal composition profile, and changes in the nuclear burning phases, are linked with the radial evolution of the star. Our results are in general agreement with previous studies. Our findings can be summarized as follows:

1. Stars stripped in binaries can swell up to giant sizes, despite having lost most of their envelope, as has also been shown in earlier studies. This implies they can fill their Roche lobe again (sometimes referred to as Case BB or BC mass transfer).

2. The maximum radius achieved strongly depends on whether or not these stars retain a hydrogen layer. At solar metallicity, mass stripping by Roche-lobe overflow is effective in removing most of the hydrogen envelope. Winds play a minor role, but they are strong enough to remove the remaining hydrogen layer. At low metallicity, all our models retain a significant hydrogen-rich envelope until the end of their evolution, which is in agreement with earlier findings.

3. At high metallicity, we find that the maximum radius (10–100 R_⊙) is inversely proportional to the stellar mass. It is only for the low-mass end (progenitors with initial masses between 8 and 9 M_⊙) that we expect stripped stars to swell enough to interact anew with their companion, which is in agreement with earlier studies.

4. At low metallicity, for case B mass transfer, stripped stars can reach sizes of up to 400–700 R_⊙, unless they fill their Roche-lobe anew. This maximum expansion is robust against variations of the wind mass loss and of the orbital period.

5. Population synthesis studies that rely on the Hurley et al. (2000) predictions or on interpolation of grids of detailed single stellar models do not properly account for the structure of stripped stars, in particular, the effect of a remaining hydrogen layer. We find large discrepancies at solar metallicity in the mass range important for neutron star progenitors. At low metallicity we find discrepancies for the full mass range and we expect progenitors of both neutron stars and black holes to be affected.

6. We estimate, with simple assumptions, that population synthesis studies under-predict the number of systems that interact by a factor of two at solar metallicity (for stripped stars of about 2–3 M_⊙). For low metallicity the discrepancy is much worse. The fraction of systems that re-interact is underestimated by a factor of 10–30 (for stripped stars of about 2–6 M_⊙).

7. We draw attention to an additional expansion phase that occurs exclusively in low-metallicity models shortly after central helium exhaustion. This phase is associated with hydrogen shell burning and only lasts a few thousand years. The star shrinks again briefly once helium shell-burning dominates, followed by the final expansion phase.

8. The first radius expansion we find at low metallicity may also have important consequences. The low binding energy of the hydrogen envelope suggests that mass may be ejected only a few tens of thousands of years prior to the final explosion, giving rise to a hydrogen-rich shell around the progenitor. This could impact the observable properties of the resulting supernova.

9. Many stripped stars are expected to be filling their Roche lobe at the moment their core collapses. This means that they would not be spherical at the moment they explode but instead they would have a “pear shape” enforced by their Roche lobe. This may carry interesting consequences for the observable characteristics of the final explosion.

10. Our results pose important concerns about the validity of rapid simulations of gravitational-wave sources. Our detailed simulations show metallicity effects that are not accounted for in rapid population synthesis simulations. Specifically, we expect the rates and channels for the formation of double neutron stars to be dependent on metallicity, in contrast to recent claims.

We anticipate further progress in this area as new observational constraints become available for core-collapse supernovae from transient surveys and for double compact-object mergers from gravitational-wave detectors. Robust model predictions will be needed to interpret these observations and learn about the physics of their binary progenitors. Our findings call for detailed investigations to better understand the outcome of additional mass-transfer phases in binary systems that have already experienced previous interactions.

Acknowledgments

We thank the referee Christopher Tout for detailed comments. This work has benefited from helpful input by M. Renzo and E. Zapartas. The authors further acknowledge N. Afsariardchi, J. J. Eldridge, A. Gilkis, J. Klencki, C. Matzner, T. Moriya, A. Piro, O. R. Pols, T. Tauris, J. Yongje, and S.-Ch. Yoon for useful discussions. 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 Organisation for Scientific Research (NWO) as part of the Vidi research program BinWaves with project number 639.042.728. YG acknowledges the funding from the Alvin E. Nashman fellowship for Theoretical Astrophysics. RF is supported by the Netherlands Organisation for Scientific Research (NWO) through a top module 2 grant with project number 614.001.501 (PI de Mink). The simulations were computed on the Dutch national e-infrastructure (Cartesius, project number 17234) with the support of the SURF Cooperative.

References

Appendix A: Parameters of the models

We provide relevant parameters of our models at solar and low metallicity in Table A.2 and A.3, respectively. Definitions for the symbols used in both tables are given in Table A.1.

Appendix B: Effects of winds on the expansion of stripped stars at low metallicity

In Sects. 3.2 and 4.2, we discuss the expansion of stripped stars. At low metallicity, the models retain a hydrogen-rich layer. The mass of that layer is linked to the maximum radius each star can achieve. However, the mass of the hydrogen-rich envelope retained by a stripped star is not only determined by the binary interaction, but can also be affected by wind mass loss. The wind mass-loss rates for stripped stars are not well known (Yoon 2015; Götberg et al. 2017, 2018; Gilkis et al. 2019). In most stellar evolution models for stripped stars, winds are typically assumed to follow the empirically-derived prescription for Wolf-Rayet stars from Nugis & Lamers (2000). At low metallicity, the effect of winds is expected to be limited due to the metallicity dependence of line-driven winds (Vink & de Koter 2005; Mokiem et al. 2007). This expectation affects predictions for the mass of the leftover hydrogen-rich layer at low metallicity and, consequently, the maximum radii of stripped-star models. We discuss the uncertainties introduced by this assumption below, showing comparisons for both our example model from Sect. 3 and the most massive model from our grid.

Figures B.1 and B.2 demonstrate the effects of increasing our assumed wind-loss rates on our exemplary model, with an initial mass of 10.5 M_⊙. Figure B.1 shows the evolution of the stellar radius after core helium burning as a function of time. We compare our default model to two models for which the mass-loss rate is increased by constant factors of three and ten. We also present a model with the mass-loss scheme behaving as if the star was at solar metallicity. For the first 20 000 years of the radius evolution, no significant impact of the wind mass-loss rate can be observed. Differences appear towards the end of the evolution. The model with a factor of ten higher mass-loss rate ends its evolution with a smaller radius, of 417 R_⊙, compared to the 515 R_⊙ of the model with our default wind mass-loss scheme. The model with the solar-metallicity wind mass-loss scheme ends its life with a radius of 516 R_⊙, a value that is very similar to that of the default model.

Fig. B.1. Radius evolution of low-metallicity stripped-star models with an initial mass of 10.5 M⊙ as a function of time after core helium depletion. The blue curve represents a model computed with our default wind scheme, while the orange and green curves are for models with mass-loss rates 3 and 10 times higher than the default, respectively. The dashed red curve represents a low-metallicity model computed with the solar metallicity mass-loss scheme.

Fig. B.2. Evolution of the wind-mass loss rate for the low-metallicity 10.5 M⊙ stripped-star model as a function of time after core helium depletion. The line styles are associated with the same models as in Fig. B.1.

Figure B.2 shows the evolution of the wind mass-loss rates of the models in Fig. B.1. For all models, the dominant mass loss is late in the evolution, when the stars have become giants. The model with a solar-metallicity wind mass-loss rate closely follows the ten-times-higher mass-loss rate for the first 20 000 yr. This is due to the metallicity dependence of the Vink et al. (2001) mass-loss rate, which is (Z)^0.85, thus, approximately ten times higher for solar metallicity. The sharp changes in the mass-loss rates can be attributed to the bistability jumps of the winds. Because the temperature at which these take place is metallicity-dependent (Vink et al. 2001), they occur earlier for the solar-metallicity model. Towards the end of the evolution, the model assuming solar-metallicity winds follows the same mass-loss rate as our default model because we assume these cool-star winds to be metallicity-independent. However, even if our model predictions had been close to the green curve our conclusions would have remained qualitatively unchanged.

Since the effect of wind-mass loss is mass-dependent, we also investigate the impact of increased wind mass-loss rates on the highest-mass model in our grid, which has an initial mass of 15 M_⊙. We present the results of this test in Fig. B.3. The impact of changing these assumptions on the radius evolution is more pronounced throughout the evolution than for the 10.5 M_⊙ model, but the differences in final radii are still not enough to affect our qualitative conclusions. The largest relative differences can be observed at the peak of the first radius expansion phase, with radii of 232, 175, 104 R_⊙ for the default wind assumptions, three-times higher, and ten-times higher mass-loss rates, respectively. However, the predicted radius at the end of the evolution is less affected. The model with a ten-times higher wind mass-loss rate has a final radius only 15% lower than the default model (and only 3% lower than the default model for the three-times higher wind mass-loss rate variation). Again, we also show the radius evolution of a model assuming solar-metallicity winds. The majority of the evolution again follows the same trend as the model with a ten-times-higher wind mass-loss rate. At the peak of the first radius expansion, the solar-metallicity wind model has a radius of 127 R_⊙. However, the final radius is only 2% smaller than that of the default mass-loss model.

Fig. B.3. Same as Fig. B.1 for the low-metallicity model with an initial mass of 15 M⊙.

From these tests, we conclude that the effect of reasonable uncertainties in wind mass loss on the maximum radius these stars are predicted to reach is very small.

Appendix C: Effects of the orbital period on the expansion of stripped stars at low metallicity

In Appendix B, we demonstrate that plausible ranges of stellar winds have only a very small impact on the expansion of stars stripped in binaries at low metallicity. The binary interaction is mainly responsible for determining the mass of the remaining hydrogen layer at low metallicity. Yoon et al. (2017) explored how the choice of orbital separation has a large impact on the final effective temperature and radius of stripped stars at low metallicity, for stars that retain a hydrogen layer with a mass that exceeds 0.15 M_⊙. In this section, we discuss how our choice of the orbital separation at low metallicity affects our results. We first show results for models with one initial stellar mass, which are stripped in a binary computed at multiple orbital separations. We then compare the maximum radii obtained for the grid of models presented in the main text with two grids at alternative orbital separations.

Figures C.1 and C.2 show the impact of the initial orbital separation on the remaining hydrogen mass and radius evolution for stellar models with the same initial mass of 11.3 M_⊙. These results are based on a set of calculations with logarithmically-spaced orbital periods ranging from 1 d to 100 d. We did not find converged solutions for the model with an initial orbital period of 1 d. We plot the full evolution of all models on the Hertzsprung–Russell diagram in Fig. C.3. The model with the shortest orbital period interacts with its companion as early as in the course of the main-sequence (case A mass transfer), which is why it has a distinct evolutionary track. We focus on the case B mass-transfer models and thus do not discuss this model below. In Fig. C.1, we display the total hydrogen mass at core helium depletion as a function of the initial orbital period. As expected, the mass of the hydrogen-rich layer increases for longer orbital periods, but it is still within a factor of two.

Fig. C.1. Total hydrogen mass at the moment of core helium depletion for low-metallicity models with an initial mass of 11.3 M⊙, shown as a function of initial orbital period, Pi.

Fig. C.2. Radius evolution of low-metallicity stripped star models with an initial mass of 11.3 M⊙ as a function of time after core helium depletion. Colors indicate the initial orbital period as specified in the legend.

Fig. C.3. Full evolution of models with an initial mass of 11.3 M⊙ with varying initial orbital periods on the Hertzsprung-Russell diagram. Colors indicate the initial orbital period as specified in the legend. All models with an initial orbital period longer than three days interact with their companions after leaving the main-sequence (case B mass-transfer).

Figure C.2 demonstrates the consequent effect of the remaining hydrogen-rich layer on the radius evolution of these low-metallicity stripped stars after core helium depletion (cf. Fig. 4). The models display a parallel radius evolution, where models with longer orbital periods have larger radii overall. However, all models reach similar final radii (between 514 and 580 R_⊙).

Because the radius evolution also depends on the total stellar mass, we compare the maximum radii of our default model grid at low metallicity with two model grids computed with initial orbital periods of 5 d and 35 d, respectively. Results from these grids are displayed in Figs. C.4 and C.5 (for comparison with Figs. C.1 and C.2, we note that the post-stripping masses from the models shown in those figures ranges from 4 to 4.35 M_⊙.) Figure C.5 presents the total hydrogen mass at the moment of core helium depletion as a function of mass. All models display a linear trend of increasing hydrogen mass as a function of mass, except for the highest mass models, which have slightly lower hydrogen mass that can be attributed to the increased effect of stellar winds. As discussed in Appendix B, the impact of stellar winds on the evolution is very small at this low metallicity. However, since the hydrogen layer has such a low mass (below 0.3 M_⊙), even a small effect can have large consequences. The total hydrogen masses for the models with initial orbital periods of 25 and 35 d are very similar, and the models with a shorter orbital periods have only slightly smaller total hydrogen masses. The total hydrogen masses are notably smaller for the models with an orbital period of 5d. We demonstrate the impact of the orbital period on the maximum radius in Fig. C.4, where we find very similar trends as in the total hydrogen mass. Overall, models with the same initial mass and for these different orbital periods reach very similar final radii.

Fig. C.4. Total hydrogen mass at the moment of core helium depletion as a function of stellar mass. Green diamonds, orange squares, and blue circles indicate grids computed with fixed orbital periods of 5, 25, and 35 d, respectively.

Fig. C.5. Maximum radius as a function of mass for the same models as in Fig. C.5.

From a population perspective, assuming the initial orbital periods are distributed uniformly in log space, only a minority of the stars in this mass range stripped by stable mass transfer in the Hertzsprung Gap would have initial orbital periods that lead to final radii which are marginally (within 100 R_⊙) different from the predictions of the grid shown in the main text. This effect does not impact our main finding, which is that the radii of such stripped stars are severely underestimated in population synthesis models, especially at low metallicity. We conclude that for case B mass-transfer, our choice of the orbital period has only a small effect on the maximum radius of stripped stars and is representative for the population. For stripped stars created by other channels (e.g., case A mass-transfer or common envelope evolution), the maximum radius at low metallicity could be significantly different because the mass of the remaining hydrogen-rich layer could be much smaller.

Appendix D: Analytic fitting functions for the radius of stripped stars

For the sake of convenience, we provide simple analytic fits to the minimum and maximum stellar radii from our stripped-star models. These fits accurately represent our models to within a few percent.

In the following formulae, M is the total mass of the stripped star at the relevant time, that is, at the minimum radius following post-stripping contraction or at the maximum radius towards the end of the nuclear burning. The difference between those masses is small for these relatively low-luminosity models and our adopted winds.

At solar metallicity, the minimum radius, R_min, for each model stripped star is well-described by:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill log\left(\frac{R_{\min }}{{\mathrm{R}}_{\odot }}\right)=max\{& -0.069\left(\frac{M}{M_{\odot }}\right)-0.23,\hfill \\ \hfill & -0.0142{\left(\frac{M}{M_{\odot }}\right)}^2+0.153\left(\frac{M}{M_{\odot }}\right)-0.744\}.\hfill \end{array}\end{array}$$(D.1)

This fit is shown in Fig. D.1, together with the ratio of the model to the fit (in all cases these ratios are for linear quantities, for example, the computed model radius in R_⊙ over fitted radius in R_⊙).

Fig. D.1. Minimum radii of stripped-envelope stars as a function of their total masses at solar metallicity. The best fit obtained (Eq. (D.1)) is shown in purple and letters indicate the best fit parameters obtained. Lower panel: residuals of the fit, which are defined as fractional (i.e., dimensionless).

The maximum radii, R_max, of the solar-metallicity stripped-star models above 2.28 M_⊙ can be expressed by:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill log\left(\frac{R_{\max }}{{\mathrm{R}}_{\odot }}\right)=min\{& 1258exp(-2.64\left(\frac{M}{M_{\odot }}\right))+0.67,\hfill \\ \hfill & -0.63\left(\frac{M}{M_{\odot }}\right)+3.76\}.\hfill \end{array}\end{array}$$(D.2)

This fit is presented in Fig. D.2. The maximum radius of the most massive progenitor was excluded from this fit due to numerical uncertainties in this model; otherwise, this fit has a maximum deviation from our models of about 2 R_⊙. For stripped stars with solar metallicity in our grid of 2.28 M_⊙ and below, we find setting log(R_max/R_⊙) = 2.8 is reasonable. The model with a radius of log(R/R_⊙) = 2.3 did not reach core carbon ignition, and so was excluded from the fit.

Fig. D.2. Maximum radii of stripped-envelope stars as a function of their total masses at solar metallicity, extended to lower mass models. The fit is shown in purple. Lower panel: fractional residuals of the fit.

For our lower-metallicity models, we fit the minimum and the maximum radii for each model with:

$$\begin{array}{c}\hfill log\left(\frac{R_{\min }}{{\mathrm{R}}_{\odot }}\right)=-9.8\times {10}^{-8}{\left(\frac{M}{M_{\odot }}\right)}^8+0.015{\left(\frac{M}{M_{\odot }}\right)}^2+0.013,\end{array}$$(D.3)

$$\begin{array}{c}\hfill log\left(\frac{R_{\max }}{{\mathrm{R}}_{\odot }}\right)=-0.016{\left(\frac{M}{M_{\odot }}\right)}^2+0.21\left(\frac{M}{M_{\odot }}\right)+2.2.\end{array}$$(D.4)

This fit to R_min is shown in Fig. D.3. The fit to R_max is in Fig. D.4, and matches our low-metallicity models to within approximately 1 R_⊙.

Fig. D.3. Minimum radii of stripped-envelope stars as a function of their total masses at low metallicity. The best fit obtained (Eq. (D.3)) is shown in red. Lower panel: residuals of the fit.

Fig. D.4. Maximum radii of stripped-envelope stars as a function of their total mass at low metallicity. The best fit obtained (Eq. (D.4)) is shown in red. Lower panel: residuals of the fit.

For these low-metallicity models, it would be mistaken to take the initial stripped-star mass as a helium-core mass since the helium core grows as a result of the helium produced by hydrogen shell burning. So for these low-metallicity models, we also provide a fit for the initial $M_{\mathrm{He}\mathrm{core}}^{\min }$ and final $M_{\mathrm{He}\mathrm{core}}^{\max }$ helium core masses, as:

$$\begin{array}{c}\hfill \left(\frac{M_{\mathrm{He}\mathrm{core}}^{\min }}{M_{\odot }}\right)=0.91\left(\frac{M}{M_{\odot }}\right)-0.047,\end{array}$$(D.5)

$$\begin{array}{c}\hfill \left(\frac{M_{\mathrm{He}\mathrm{core}}^{\max }}{M_{\odot }}\right)=0.89\left(\frac{M}{M_{\odot }}\right)-0.037.\end{array}$$(D.6)

These fits are shown in Fig. D.5 and D.6, respectively.

Fig. D.5. Helium core mass at the minimum radius as a function of their total masses at low metallicity. The best fit obtained (Eq. (D.5)) is shown in red. Lower panel: residuals of the fit defined as the difference between the fit and the models in solar masses.

Fig. D.6. Helium core mass at the maximum radius as a function of the total mass at low metallicity. The best fit obtained (Eq. (D.6)) is shown in red. Lower panel: residuals of the fit defined as the difference between the fit and the models in solar masses.

We intend our fits to be used as a first improvement to current population synthesis calculations, especially focusing on predictions for double neutron-star systems. However, we do not provide a full way to integrate with the Hurley prescriptions. Our detailed calculations show that stars at low metallicity are partially stripped and exhibit intermediate behavior, which does not fit into the powerful but relatively simple scheme provided by Hurley and collaborators. Although one could add an extra parameter to keep track of the mass of the remaining envelope, this would still not be sufficient to, for example, predict the response of such stars to mass loss. More extended grids of detailed models would be needed to encompass multiple variations of the orbital period, the metallicity, and the mass-transfer efficiency, and for an even wider mass range than we present here. The use of such dense grids of stellar models for population studies would be preferable to current population synthesis models given the approximations we point out in Sect. 5.

Appendix E: Binding energy of the envelope

In Sect. 4.3, we present the envelope binding energies computed from our models. However, the magnitude of the envelope binding energy is sensitive to the location of the boundary between core and envelope and to whether or not internal energy terms are incorporated in the calculation (see, e.g., Ivanova et al. 2013). We discuss this further below.

E.1. Binding energy with and without internal energy

When considering common-envelope evolution, the term “binding energy” is inconsistently used in the literature. Sometimes the term is used purely for the gravitational binding energy, without including either the thermal internal energy of the envelope material or the electrostatic potential energy of ionized and dissociated matter. In this context it is common to refer to the thermal energy and recombination energy terms (including molecular dissociation energy) collectively simply as “internal energy”, as we do. The difference between binding energies which ignore or include these terms can be substantial, and can qualitatively affect common-envelope outcomes (e.g., Han et al. 1994; Podsiadlowski et al. 2003; Ivanova et al. 2013, and references therein). In Fig. E.1, we compare the effect of including or excluding the internal energy terms for the calculation of the binding energy for our grids at high and low metallicity. Models for which the binding energy was computed with the internal energy terms have a lower magnitude of binding energy – that is, they are less bound – and span a smaller absolute range of binding energies at the beginning of the evolution.

Fig. E.1. Binding energy of the hydrogen-rich envelope as a function of the radius after core helium depletion at solar (top) and lower (bottom) metallicity. Here, we compare the binding energy computed with and without the internal energy terms. On the top axis, we indicate the orbital separation at which a neutron star companion would be expected if the star would fill its Roche lobe.

The dimensionless λ parameter is commonly employed to encode how the structure of an envelope affects its binding energy. This was introduced by de Kool (1990) for calculating the outcome of common-envelope evolution with the “alpha prescription”. It is defined as

$$\begin{array}{c}\hfill \lambda =\frac{G{M}_1{M}_{1,\mathrm{env}}}{E_{\mathrm{B}}{R}_1},\end{array}$$(E.1)

where M₁ is the mass of the primary star transferring mass, M_1, env the mass of its envelope, and R₁ its radius. We show the λ parameter at high and low metallicity for both definitions of the binding energy in Fig. E.2.

Fig. E.2. Envelope structure parameter λ as a function of the radius after core helium depletion at solar (top) and lower (bottom) metallicity. The dashed line indicates λ parameters computed using only the gravitational energy term for the binding energy.

E.2. Choice of core/envelope mass boundary

Subtle differences in the definition of the core and envelope boundary mass can also significantly affect the value of the binding energy (see, e.g., Han et al. 1994; Tauris & Dewi 2001; Podsiadlowski et al. 2003; Ivanova et al. 2013). We do not investigate the effect of choosing such different definitions for the core-envelope boundaries on our results.

However, for these stripped stars we do investigate the differences between choosing a hydrogen-rich or helium-rich boundary. We compare the binding energy computed using Eq. (1) for the hydrogen-rich (H) and the helium-rich (H + He) envelope in Fig. 5 as a function of time after core helium depletion at high and low metallicities. For both grids, the binding energy of the H-envelope is an order of magnitude lower than that of the the H + He envelope. This is unsurprising, given that the potential well is deeper and steeper closer to the core of the star. For both regions, the magnitude of the total absolute binding energy increases with increasing initial mass. At solar metallicity, only the lowest mass models have a hydrogen envelope shortly after core helium depletion before it disappears due to wind mass loss after about 30 kyr. In contrast, at low metallicity, all models retain a hydrogen envelope. At the moment of the first radius expansion, the binding energy of the H-envelope drops before increasing shortly, and decreasing again, while that of the H + He envelope increases, and later decreases. This is consistent with the double-mirror effect.

All Tables

All Figures

Fig. 1. Evolution of a high (left) and low (right) metallicity star of initially 10.5 M⊙ that is stripped due to binary interaction. Letters A to J mark evolutionary points discussed in Sect. 3. Top: evolutionary tracks on the Hertzsprung-Russell diagram. The color indicates the total mass of hydrogen present at a given moment. For comparison, we show lines of constant radii as well as the evolution of the alternate-metallicity model in gray. Middle: evolution of the stellar structure of the stars in mass coordinate (Kippenhahn diagrams), given as a function of time until core collapse (τcc). Convective and overshooting regions are marked with double- and single-hatched regions, respectively. Colors indicate zones dominated by nuclear burning (yellow) or neutrino cooling (purple) where ϵ = sign(ϵnuc−ϵν)log10(max(1.0,|ϵnuc−ϵν|)/[erg g−1 s−1]). Here ϵnuc is the nuclear energy generation rate and ϵν the neutrino energy. The black line indicates the location of the stellar surface. Bottom: luminosity produced by hydrogen-, helium- and metal-burning as a fraction of the total luminosity produced by nuclear reactions (upper) and stellar radius (lower) as a function of time until core-collapse.

In the text

Fig. 2. Kippenhahn diagrams showing the evolution of the stellar structure for our two example models at solar (left panels) and low (right panels) metallicity. They are shown using mass coordinate as the vertical axis (top panels) and radius coordinate (bottom panels). The horizontal axis indicates the time since helium depletion up to core carbon depletion. The black line represents the surface of the star. Regions of mixing by convection and overshooting are marked with double- and single-hatches, respectively. Colors indicate zones dominated by nuclear burning or neutrino cooling. See Fig. 1 and Sect. 3 for details and a discussion.

In the text

Fig. 3. Evolutionary tracks of stars at solar (top) and sub-solar metallicity (bottom) showing their luminosity, L as a function of effective temperature, Teff. Lines of constant radius are provided for reference. We highlight the evolution after core helium depletion in color, where the color refers to the initial mass. Yellow star symbols indicate core carbon exhaustion, and white star symbols denote the final time step calculated for models that did not reach core carbon exhaustion. The white circles in the lower panel indicate when the radius reaches a first maximum after core helium depletion. Above each panel, we indicate approximate ranges of effective temperatures typical for Wolf Rayet stars (WR), blue, yellow, and red supergiants (BSG, YSG and RSG respectively). The horizontal color scale above each panel indicates, for each effective temperature, the color as perceived by the human eye.

In the text

Fig. 4. Evolution of the radius as a function of time after core helium depletion for the stars in our model grids at solar (left) and lower (right) metallicity. From top to bottom, colors indicate the remaining mass of hydrogen and the normalized luminosity from nuclear burning of hydrogen, helium, and heavier elements; see text for details. More massive stars evolve faster and are thus located towards the left of the figures. White crosses indicate the moment of core carbon ignition. Yellow star symbols indicate core carbon exhaustion, and white star symbols denote the final time step calculated for models that did not reach core carbon exhaustion.

In the text

	Fig. 5. Binding energy of the envelope as a function of radius after core helium depletion at solar (top panel) and low (bottom panel) metallicity. We show values for the hydrogen-rich envelope alone, as well as the hydrogen and helium envelope combined, as labelled. Colors indicate the mass of the stripped star.
In the text

Fig. 6. Comparison of our results (shown in blue in each panel) with (1) the widely-used prescriptions by Hurley et al. (2000; shown in grey in the panels on the top row) and (2) the PARSEC models (Bressan et al. 2012) used in the SEVN code Spera et al. (2019; shown in grey in the panels on the bottom row). Panels on the left: comparison at solar metallicity (or its closest available equivalent, i.e. Z = 0.02 for SEVN) and panels on the right: comparison for low metallicity (Z = 0.0002 for SEVN). We show the maximum radius (circles) and the minimum radius (triangles). For our models, we also mark the radius at core carbon ignition (squares).

In the text

Fig. 7. Top: binding energy of the envelope plotted as a function of radius. On the top x-axis we have indicated the orbital separation assuming a neutron star companion with a canonical mass of 1.4 M⊙, for reference. Our solar models are shown in the top left panel where we show the binding energy of the helium envelope. Our low-metallicity models are shown in the top-right panel, where we show the binding energy of the hydrogen envelope (cf. Sect. 3.3 and Fig. 5). Bottom: final separation after the common envelope evolution of stripped stars with neutron star companions as a function of the initial separation at solar (left) and low (right) metallicity. Stars with a pre-CE convective envelope are marked with colors, which indicate the mass of the stripped star, assuming a common envelope efficiency parameter αCE = 1.

In the text

Fig. B.1. Radius evolution of low-metallicity stripped-star models with an initial mass of 10.5 M⊙ as a function of time after core helium depletion. The blue curve represents a model computed with our default wind scheme, while the orange and green curves are for models with mass-loss rates 3 and 10 times higher than the default, respectively. The dashed red curve represents a low-metallicity model computed with the solar metallicity mass-loss scheme.

In the text

	Fig. B.2. Evolution of the wind-mass loss rate for the low-metallicity 10.5 M⊙ stripped-star model as a function of time after core helium depletion. The line styles are associated with the same models as in Fig. B.1.
In the text

	Fig. B.3. Same as Fig. B.1 for the low-metallicity model with an initial mass of 15 M⊙.
In the text

	Fig. C.1. Total hydrogen mass at the moment of core helium depletion for low-metallicity models with an initial mass of 11.3 M⊙, shown as a function of initial orbital period, Pi.
In the text

	Fig. C.2. Radius evolution of low-metallicity stripped star models with an initial mass of 11.3 M⊙ as a function of time after core helium depletion. Colors indicate the initial orbital period as specified in the legend.
In the text

	Fig. C.3. Full evolution of models with an initial mass of 11.3 M⊙ with varying initial orbital periods on the Hertzsprung-Russell diagram. Colors indicate the initial orbital period as specified in the legend. All models with an initial orbital period longer than three days interact with their companions after leaving the main-sequence (case B mass-transfer).
In the text

	Fig. C.4. Total hydrogen mass at the moment of core helium depletion as a function of stellar mass. Green diamonds, orange squares, and blue circles indicate grids computed with fixed orbital periods of 5, 25, and 35 d, respectively.
In the text

	Fig. C.5. Maximum radius as a function of mass for the same models as in Fig. C.5.
In the text

	Fig. D.1. Minimum radii of stripped-envelope stars as a function of their total masses at solar metallicity. The best fit obtained (Eq. (D.1)) is shown in purple and letters indicate the best fit parameters obtained. Lower panel: residuals of the fit, which are defined as fractional (i.e., dimensionless).
In the text

	Fig. D.2. Maximum radii of stripped-envelope stars as a function of their total masses at solar metallicity, extended to lower mass models. The fit is shown in purple. Lower panel: fractional residuals of the fit.
In the text

	Fig. D.3. Minimum radii of stripped-envelope stars as a function of their total masses at low metallicity. The best fit obtained (Eq. (D.3)) is shown in red. Lower panel: residuals of the fit.
In the text

	Fig. D.4. Maximum radii of stripped-envelope stars as a function of their total mass at low metallicity. The best fit obtained (Eq. (D.4)) is shown in red. Lower panel: residuals of the fit.
In the text

	Fig. D.5. Helium core mass at the minimum radius as a function of their total masses at low metallicity. The best fit obtained (Eq. (D.5)) is shown in red. Lower panel: residuals of the fit defined as the difference between the fit and the models in solar masses.
In the text

	Fig. D.6. Helium core mass at the maximum radius as a function of the total mass at low metallicity. The best fit obtained (Eq. (D.6)) is shown in red. Lower panel: residuals of the fit defined as the difference between the fit and the models in solar masses.
In the text

	Fig. E.1. Binding energy of the hydrogen-rich envelope as a function of the radius after core helium depletion at solar (top) and lower (bottom) metallicity. Here, we compare the binding energy computed with and without the internal energy terms. On the top axis, we indicate the orbital separation at which a neutron star companion would be expected if the star would fill its Roche lobe.
In the text

	Fig. E.2. Envelope structure parameter λ as a function of the radius after core helium depletion at solar (top) and lower (bottom) metallicity. The dashed line indicates λ parameters computed using only the gravitational energy term for the binding energy.
In the text