Formation of GW190521 from stellar evolution: the impact of the hydrogen-rich envelope, dredge-up, and ¹²C(α, γ)¹⁶O rate on the pair-instability black hole mass gap Free

ABSTRACT

1 INTRODUCTION

The LIGO-Virgo collaboration (LVC) recently reported the discovery of GW190521 (Abbott et al. 2020b,e). With a primary black hole (BH) mass of |$85^{+21}_{-14}$| M_⊙ and a secondary BH mass of |$66^{+17}_{-18}$| M_⊙ (90 per cent credible intervals), this binary black hole (BBH) merger is the most massive one observed with gravitational waves to date (Abbott et al. 2016a,b,c, 2019a,b, 2017, 2020a,c,d). The mass of the primary BH is a puzzle for astrophysicists, because it lies in the middle of the pair-instability (PI) mass gap.

PI is possibly the key process to understand the maximum mass of stellar-origin BHs, because it is expected to carve a gap in the mass spectrum of BHs between ≈40–65 and ≈120 M_⊙ (e.g. Heger & Woosley 2002; Woosley, Blinnikov & Heger 2007; Belczynski et al. 2016; Spera & Mapelli 2017; Woosley 2017, 2019; Giacobbo, Mapelli & Spera 2018; Mapelli et al. 2019; Spera et al. 2019; Tanikawa et al. 2020a). The uncertainty on the boundaries of the PI mass gap depends on a number of physical processes (e.g. Takahashi 2018; Leung, Nomoto & Blinnikov 2019; Marchant et al. 2019; Stevenson et al. 2019; Farmer et al. 2019, 2020; Mapelli et al. 2020; Renzo et al. 2020; Tanikawa et al. 2020b; van Son et al. 2020).

PI starts to be efficient when the plasma inside stars reaches temperatures above 6 × 10⁸ K and densities between about 10² and 10⁶ g cm⁻³. Such conditions on internal temperature and density are reached when massive CO cores are formed in stars. Pair creation converts thermal energy into rest mass of |$e^-\, e^+$|⁠, lowering the central pressure. This causes an hydro-dynamical instability that leads to core collapse. In the collapsing core, oxygen is ignited explosively.

In stars developing a helium core mass ≈32–64 M_⊙ at the end of carbon burning, PI manifests as pulsational PI (PPI): the star starts to pulsate losing a large amount of mass, until it finds a new equilibrium. After PPI, the star will end its life with a core-collapse supernova or a direct collapse, leaving a BH. More massive stars, with helium core mass ≈64–135 M_⊙, undergo a pair-instability supernova (PISN): a powerful single pulse that totally destroys the star, leaving no compact remnant. Finally, PI triggers the direct collapse of stars with He core mass >135 M_⊙, eventually producing an intermediate-mass BH (Heger & Woosley 2002). Such theoretical picture has been confirmed by many authors in both 1D (Heger & Woosley 2002; Woosley 2017; Takahashi 2018; Marchant et al. 2019) and multi-dimensional hydrodynamical simulations (Chatzopoulos, Wheeler & Couch 2013; Chen et al. 2014). From an observational perspective, iPTF14hls (Arcavi et al. 2017) and SN 2016iet (Gomez et al. 2019) cannot be explained by current core-collapse supernova models and might be connected with PI.

If a massive star enters or not in the PI regime depends mainly on the CO core mass left after the central He burning and on its evolution during the following advanced phases. For a given zero-age main sequence (ZAMS) mass (M_ZAMS) and stellar metallicity (Z), the final CO mass can be strongly affected by many evolutionary processes, such as mass-loss, convection, overshooting, semi-convection (Kaiser et al. 2020; Clarkson & Herwig 2020; Renzo et al. 2020), rotation (Chatzopoulos & Wheeler 2012; Limongi & Chieffi 2018; Song et al. 2020; Marchant & Moriya 2020; Mapelli et al. 2020) and internal and surface magnetic fields (Haemmerlé & Meynet 2019; Keszthelyi et al. 2020). These processes may change the ZAMS mass versus CO core mass relation and, thus, the initial mass at which stars enter in the PI regime (Takahashi 2018). Another important ingredient that influences both the CO core mass and its composition is represented by nuclear reactions during the core helium burning (CHeB). The two most important ones are the triple-α reaction, ⁴He(2α, γ)¹²C, and the carbon-α one, ¹²C(α, γ)¹⁶O. The interplay of these two reactions during the CHeB phase determines the final carbon-to-oxygen ratio (¹²C/¹⁶O) in the CO core, hence the final fate of the star (Weaver & Woosley 1993). The ¹²C(α, γ)¹⁶O reaction rate is one of the most uncertain (Caughlan & Fowler 1988; Buchmann 1996; Sallaska et al. 2013; deBoer et al. 2017; Rapagnani et al. 2017), because it is difficult to measure. The study on how the ¹²C(α, γ)¹⁶O reaction affects the stellar evolution is an old problem (Iben 1967; Brunish & Becker 1990; Alongi et al. 1991). In recent years, many efforts have been done to understand its importance on the most advanced stellar phases and in the final fates (Woosley & Heger 2007; Tur, Heger & Austin 2007, 2010; deBoer et al. 2017; Fields et al. 2018; Takahashi 2018; Farmer et al. 2019; Sukhbold & Adams 2020; Farmer et al. 2020).

In particular, Farmer et al. (2019) and Farmer et al. (2020) show that the ¹²C(α, γ)¹⁶O reaction rate is one of the main sources of uncertainty on the boundaries of the PI mass gap, when pure helium stars are considered. Another major source of uncertainty on the final BH mass is the fate of the residual hydrogen envelope in case of a direct collapse: if a massive star retains a fraction of its H envelope to the very end of its life, is this envelope able to collapse together with the internal layers of the star? Or is it so loosely bound that it is expelled, even without a supernova explosion (Sukhbold et al. 2016)? If the residual H envelope collapses to the final BH, the lower edge of the PI mass gap (M_gap) might increase by ≈20 M_⊙ (Mapelli et al. 2020).

The main focus of this work is to study the effect of the uncertainty of the ¹²C(α, γ)¹⁶O reaction rate on the M_gap in metal-poor stars. With respect to Farmer et al. (2020), who consider only pure-helium stars, we also investigate stars with hydrogen envelopes. This enables us to model the PI mass gap accounting for the possible collapse of the residual H envelope. Moreover, we find that the onset of a dredge-up during the H shell burning phase plays a crucial role for the evolution of the CO core of the most massive stars. The PI mass gap might even disappear for the combined effect of envelope undershooting and low ¹²C(α, γ)¹⁶O reaction rate. Our goal is to determine under which conditions the formation of BHs like the primary component of GW190521 (Abbott et al. 2020b,e) is possible through single stellar evolution.

In Section 2, we give a description of the parsec code, we describe the method adopted to include different ¹²C(α, γ)¹⁶O reaction rates and how we test the stellar dynamical stability during the evolution. In Section 3, we present the new parsec evolutionary tracks of pure-He stars computed with varying ¹²C(α, γ)¹⁶O reaction rates and compare them with the results by other authors. In Section 4, we present the new tracks of massive stars with hydrogen envelopes and compare the maximum masses found for the M_gap with other studies. Finally, in section 6 we draw our conclusions.

2 METHODS

In this work, we use the parsec V2.0 code that has been extensively described in Costa et al. (2019) and references therein. In the following, we describe the major updates adopted for this work.

2.1 Winds of massive stars

Mass-loss by stellar winds is modelled as described in Chen et al. (2015). For massive hot stars, we take into account the dependence of mass loss on both stellar metallicity (Vink, de Koter & Lamers 2000, 2001) and Eddington ratio (Gräfener & Hamann 2008; Vink 2011). With this mass loss, the models are able to reproduce the Humphreys-Davidson limit (Humphreys & Davidson 1984) observed in the Galactic and Large Magellanic Cloud (LMC) colour-magnitude diagrams¹ (Chen et al. 2015). The metallicity dependence of mass loss is expressed as a dependence on the surface iron abundance.²

These two factors are obtained by fitting models of WN and WC/WO stars in which mass-loss mainly depends on iron and, to a lesser extent, on CNO elements (Vink 2015). When metallicity decreases, the mass-loss rate lowers and flattens to values given by the CNO elements (see fig. 2 in Vink 2015). It is unlikely that the WR mass-loss depends on the helium surface content which, instead, is likely the result of a high mass-loss rate (Bestenlehner et al. 2014).

2.2 Opacity, neutrinos, and equation of state

In the high-temperature regime, 4.2 ≤ log (T/K) ≤ 8.7, we adopt the opacity tables provided by the Opacity Project At Livermore (OPAL)³ team (Iglesias & Rogers 1996, and references therein), in the low-temperature regime, 3.2 ≤ Log (T/K) ≤ 4.1, we employ the æsopus tool⁴ (Marigo & Aringer 2009). Conductive opacities are included following Itoh et al. (2008).

Energy losses by electron neutrinos are taken from Munakata, Kohyama & Itoh (1985) and Itoh & Kohyama (1983), and for plasma neutrinos we use the fitting formulae by Haft, Raffelt & Weiss (1994).

To compute the equation of state (EOS) of stellar matter, we adopt two different models, depending on central temperature. For temperature Log (T/K) < 8.5, we adopt the freeeos  ⁵ code version 2.2.1 by Alan W. Irwin. As already described in Bressan et al. (2012), this code is integrated within our evolutionary code so we can obtain the EOS in two ways: we can compute it on-the-fly, with a higher degree of accuracy, but with a higher computational cost; or we can retrieve it by means of pre-computed look-up tables with different combinations of abundances. As shown in Bressan et al. (2012), this second method is accurate enough and we prefer it for numerical speed reasons. At a given partition of elements (X_i/Z) and metallicity, the tables are divided in two sets: (i) the H-rich set in which there are 10 tables with a varying hydrogen and helium abundances, and (ii) the H-free one composed of 12 tables with a varying helium, carbon, and oxygen composition.

For temperature Log (T/K) > 8.5, we adopt the code described in Timmes & Arnett (1999), because the freeeos tables do not include the treatment of pair creation.

2.3 Nuclear reaction network

The nuclear reaction network consists in 72 different reactions and 33 isotopic elements, from hydrogen to zinc. The reaction rates and Q-values are taken from the jina reaclib database (Cyburt et al. 2010). In the network, we treat all the most important reactions from hydrogen to oxygen burning. The network includes also the reverse reactions of the α-capture. The new adopted reactions with their rates are listed in Table 1. The other reactions, already included in the previous versions of the code, can be found in Table 1 of Fu et al. (2018). In the network, we include all the reactions of the cold-CNO cycle, but not those of the hot-CNO cycle. The latter requires high temperatures (> 0.1 Gk) in H-rich regions to ignite, but these are never reached in our models.

In parsec v2.0, nuclear reaction network and element mixing are treated and solved at the same time, adopting an implicit diffusive scheme (more details in Marigo et al. 2013; Costa et al. 2019).

2.4 The ¹²C(α, γ)¹⁶O reaction

For the ¹²C(α, γ)¹⁶O reaction, we use the chw0 (v5) rate from the jina reaclib database (Fu et al. 2018), which is a re-evaluation of the Buchmann (1996) rate that includes high-lying resonances (Cyburt et al. 2012).⁶

The jina database does not include uncertainties on the chw0 (v5) rate. To estimate the uncertainties, we compared our adopted rate with the one given in the starlib database (Sallaska et al. 2013). In the starlib database, each rate is given as a mean value with a correspondent confidence error at |$1\, \sigma$|⁠. The ¹²C(α, γ)¹⁶O rate given in the starlib database is based on the rate by Kunz et al. (2002).

Fig. 1 shows the comparison of the two rates in the temperature range of the core helium burning, i.e. between about 0.1 and 0.4 GK. The mean values of the two rates are in agreement within a few percent. To reproduce the uncertainties given by the starlib  ¹²C(α, γ)¹⁶O reaction rate in the range between |$-3\, \sigma$| and |$+ 3\, \sigma$| we multiplied our adopted jina rate by a temperature independent multiplier factor, f_CO. The values of f_CO are listed in Table 2. The middle panel of Fig. 1 shows the ratios of the rates with respect to the jina standard and at different σ. In the bottom panel of Fig. 1, we also compared the adopted uncertainties with those given by deBoer et al. (2017). We find that the deBoer et al. rates are generally lower than our adopted rates in the case of |$+3\, \sigma$|⁠, |$+2\, \sigma$|⁠, and |$+1\, \sigma$|⁠. In the best value case, the two rates are very similar at low temperature but start to be different at high temperatures. Interestingly, the deBoer et al. rates with |$-1\, {}\sigma$|⁠, |$-2\, \sigma$|⁠, and |$-3\, \sigma$| show a very similar trend to the starlib ones, at temperatures > 0.2 Gk. Hence, had we used the uncertainties from deBoer et al. (2017), our main conclusions on the maximum possible mass of a BH (Section 4.5) would not have changed significantly.

Figure 1.

Upper panel: Comparison of the ¹²C(α, γ)¹⁶O rate from starlib database, deBoer et al. (2017) and jina database in red, orange, and black, respectively. The blue shaded area indicates the temperature range of helium burning in the stellar core. Middle panel: Comparison between starlib and jina reaction rates with different σ. All the rates are divided by the standard (⁠|$0\, \sigma$|⁠) jina rate. The blue, orange, and grey continuous lines are the starlib rates with |$\pm 1\, \sigma$|⁠, |$\pm 2\, \sigma$|⁠, and |$\pm 3\, \sigma$|⁠, respectively. The purple, green, and brown dotted lines indicate the jina  |$\pm 1 \, \sigma$|⁠, |$\pm 2\, \sigma$|⁠, and |$\pm 3\, \sigma$| rates, respectively. Lower panel: Comparison between deBoer et al. (2017) and jina reaction rates with different σ.

Open in new tabDownload slide

2.5 Pair creation and dynamical instability

To improve the readability of 〈Γ₁〉 listed in tables, we define the following value |$\mbox{$\langle \Gamma _{1} \rangle _{\rm -}$} \equiv {} \mbox{$\langle \Gamma _{1} \rangle _{\rm }$} - (4/3 + 0.01)$|⁠. When |$\mbox{$\langle \Gamma _{1} \rangle _{\rm -}$} \gt 0$|⁠, we consider the star globally stable, otherwise, the star is unstable. Following the evolution through hydro-dynamical phases is beyond the purpose of this paper.

3 MODELS OF PURE-HE STARS

We first analyse pure-He stars and compare our results with other studies.

3.1 Evolutionary tracks

We computed pure-helium evolutionary tracks with seven different values of the multiplier parameter, f_CO, for the ¹²C(α, γ)¹⁶O reaction rate, that correspond to varying the rate between |$-3\, \sigma$| and |$+3\, \sigma$| (Sallaska et al. 2013). Rotation is neglected in this study, since it increases the CO core and enhances the mass-loss during the evolution, with the effect to lower the final masses (Chatzopoulos & Wheeler 2012; Mapelli et al. 2020; Song et al. 2020). All tracks evolve from the helium ZAMS (He-ZAMS) to core oxygen burning, when stars remain globally stable (〈Γ₁〉_TOT- > 0).

To build the He-ZAMS, we adopted the same methodology as described by Spera et al. (2019). The new sets are composed of stars with the following masses 20, 30, 40, 50, 55, 60, 65, 70, 80, 90, and 100 M_⊙. We adopted the solar-scaled element mixture by Caffau et al. (2011), and a metallicity of Z = 0.001. We choose this particular value of metallicity to better compare our results with those found by Farmer et al. (2019), Farmer et al. (2020). For the mixing treatment, we adopt the mixing length theory (MLT, Böhm-Vitense 1958) with a fixed solar calibrated value of the mixing length, α_MLT = 1.74 (Bressan et al. 2012). We adopt the Schwarzschild criterion to define the convective regions, and a core overshooting with a mean free path of the convective element across the border of the unstable region of λ_ov = 0.4 in units of pressure scale height, H_P. Therefore, we are assuming an overshooting distance of ≈0.2 H_P above the convective core (Costa et al. 2019).

In the convective envelope, we adopt an undershooting distance of Λ_env = 0.7 H_P below the border of the unstable region. This value of Λ_env is calibrated on observations. In particular, it allows us to reproduce the location of the red giant branch (RGB) bumps and their luminosity function, and the extension of blue loops of intermediate-mass stars, in colour-magnitude diagrams of globular clusters (Alongi et al. 1991; Tang et al. 2014; Fu et al. 2018). Tang et al. (2014) showed that higher values of Λ_env might be needed in order to better reproduce the observed blue-to-red supergiants ratios in low metallicity dwarf irregular galaxies, while (Fu et al. 2018) suggested that an even higher value should be used to fit the RGB bumps of the most metal poor globular clusters.

Fig. 2 shows the evolution of three selected sets computed with f_CO corresponding to |$+ 3\, \sigma$|⁠, |$0\, \sigma$|⁠, and |$-3 \, \sigma$|⁠. Stars begin their evolution from the He-ZAMS. We define such phase as the moment in which the star is fully sustained by the nuclear reactions’ energy, and the star has burnt less than 1 per cent of the helium hosted in the core. During the helium main sequence (He-MS), all stars evolve to higher luminosity. Stars with M_ZAMS < 50 M_⊙ evolve directly toward hotter effective temperatures. Stars with M_ZAMS ≥ 50 M_⊙ first expand reducing their effective temperature and then turn to higher temperatures. As the evolution proceeds, the stars move to higher luminosity and hotter effective temperature and start to burn and deplete sequentially carbon, neon, and then oxygen in their cores.

Figure 2.

HR diagram of three sets of tracks of pure-He stars, computed with Z = 0.001. Solid black lines, dashed green lines and dotted brown lines indicate tracks computed with the ¹²C(α, γ)¹⁶O reaction rate |$-3\, \sigma$|⁠, with |$0 \, \sigma$| and with |$+3\, \sigma$|⁠, respectively. The dashed grey line shows the pure He-ZAMS. The symbols indicate helium (yellow), carbon (green), neon (red), and oxygen (blue) burning phases, respectively. Stars and circles symbols indicate the ignition and depletion of each element, respectively. The numbers on the right of the ZAMS line indicate the initial mass in M_⊙ units.

Open in new tabDownload slide

Fig. 2 shows that different rates do not affect the path followed by stars in the HR-diagram, but instead they may affect when and where stars ignite and deplete carbon, neon and oxygen in their cores. These small variations are due to different mass fractions of carbon and oxygen present in the cores at the end of the He-MS, that depend on the assumed rate of the ¹²C(α, γ)¹⁶O reaction. The evolution is eventually stopped if the star becomes globally unstable by checking the 〈Γ₁〉_TOT- value defined in equation (4). In the other case, we continue the evolution until the lowest possible value of central oxygen is reached. These values are shown in Table A1. At this point, we stop the evolution because the evolutionary times become too short and induce numerical instabilities in the code.

Fig. 3 shows the duration of the He-MS and the resulting mass of the CO core, M_CO, divided by the M_CO of the star with |$0\, \sigma$|⁠, as a function of the initial mass, M_HE-ZAMS, and for different reaction rates. Looking at the same initial mass, models computed with higher rates have slightly longer He-MS than models with lower rates. On the other hand, the M_CO is almost not affected by the different rates adopted, and the differences are below the 1 per cent in the worst case (as shown in the central panel). The main difference at the end of the He-MS between models computed with different rates is the carbon-to-oxygen ratio in the CO core, that is shown in the bottom panel. In the case of the 20M_⊙ stars, we find the largest variation, that goes from 1 (i.e. the CO core is composed 50 per cent of ¹²C and 50 per cent of ¹⁶O) to about 0.01 (i.e. the CO core is almost totally composed of ¹⁶O).

Figure 3.

Upper panel: He-MS lifetimes versus the initial mass of pure-He stars for different adopted rates. Middle panel: Ratio between the CO core masses and the CO core mass of the model computed with |$0 \, \sigma$| at the end of the He-MS. Lower panel: carbon-to-oxygen ratio in the CO cores at the end of the He-MS.

Open in new tabDownload slide

To show the main differences on the evolution of our pure-He stars due to different assumed rates, in Fig. 4 we plot the Kippenhahn diagrams of pure-He stars with M_ZAMS = 40 M_⊙ with different rates. The mass lost by the star through the evolution is just about 1.2 M_⊙, and is manly lost during the He-MS. During the He-MS, these stars build up massive convective cores that mix the star up to the 92 per cent of the entire mass. Different rates slightly affect both the luminosity of the model and the size of the convective core leading to relative differences of He-MS lifetimes up to 18 per cent between the model computed with |$+ 3\, \sigma$| and the one with |$- 3\, \sigma$| (appreciable also in Fig. 3).

Figure 4.

Kippenhahn diagrams of three models of pure-He stars with M_ZAMS = 40M_⊙ and with the ¹²C(α, γ)¹⁶O reaction rate |$- 3 \, \sigma$|⁠, |$0 \, \sigma$|⁠, and with |$+ 3\, \sigma$| in the left-hand, middle and right-hand panels, respectively. The black areas represent the convective unstable core, while the grey areas indicate the convective envelope or intermediate convective regions between the core and the surface of the star. Continuous yellow and blue lines show the stellar mass and the CO core. The yellow and green hatched areas indicate the helium and carbon burning regions, respectively. We plot zones that contribute at least for the 1 per cent of the L_He and L_C luminosity at a given time-step.

Open in new tabDownload slide

After the He-MS, stars are totally radiative and have massive CO cores that do not change much depending on the assumed rate of the ¹²C(α, γ)¹⁶O reaction, as shown in Fig. 3. On the other hand, the carbon-to-oxygen ratio is very different. The carbon (oxygen) mass fraction in the CO core are 0.43, 0.17, and ≈0.004 (0.52, 0.79, and 0.91) for stars with |$- 3\, \sigma$|⁠, |$0 \, \sigma$|⁠, and with |$+ 3\, \sigma$|⁠, respectively. Up to this stage, the evolution of these three models is quite similar.

Stars ignite carbon in their centre when temperatures reach 0.8–1 GK. The carbon burning phase lasts for about 21 years in the case of the |$- 3\, \sigma$| model, 10 years and 5 years for models with |$0 \, \sigma$| and |$+ 3\, \sigma$|⁠, respectively. During this phase, a helium burning shell forms above the CO core. The neon burning starts when the central temperatures are about 1.4 GK and all the neon in the central parts of the star is depleted in about 46, 6, and 2 days for models computed with |$- 3 \, \sigma$|⁠, |$0\, \sigma$|⁠, and |$+ 3\, \sigma$|⁠, respectively. During this phase, depending on the carbon mass fraction in the CO core, models may turn-on a carbon burning shell above the centre of the star, which develops a big intermediate convective region in the CO core. Such convective region may lead to a growth of the CO core, and provides fresh carbon to the C-burning shell. The energy released by the off-centre carbon burning sustains the external layers of the star, and avoids the stellar collapse due to the pair creation instability. The important role of the convective C-burning shell in preventing the collapse has been reported by other authors (such as Takahashi 2018; Farmer et al. 2020).

Finally, oxygen burning takes place at the centre of the star when the central temperature reaches ≈2 GK. At this point of the evolution, only the model computed with |$-3\, \sigma$| is globally stable and we can follow the evolution until oxygen’s depletion. The model with the |$0 \, \sigma$| becomes dynamically unstable at about half of the oxygen burning, while the one with |$+ 3\, \sigma$| becomes dynamically unstable at the beginning of oxygen burning. The large amount of energy released by oxygen burning generates a small convective core of about 1.63 M_⊙, 2.87 M_⊙, and 0.83 M_⊙ at the end of the computation for models with |$- 3\, \sigma$|⁠, |$0\, \sigma$|⁠, and |$+ 3\, \sigma$|⁠, respectively. If possible, all the oxygen is depleted leaving a core mainly composed of silicon and sulfur. The model with |$- 3\, \sigma$| burns the oxygen in about 30 days. At this point, the |$-3\, \sigma$| model is near to silicon burning and to the final core collapse.

Fig. 5 shows the evolution of the 〈Γ₁〉_TOT value of pure-He stars’ models with M_ZAMS = 40 M_⊙ with different rates adopted for the ¹²C(α, γ)¹⁶O reaction. This figure shows the evolution from the ignition of carbon. As already seen in Fig. 4, our models start to behave differently after carbon depletion from the core, depending on the amount of carbon-to-oxygen ratio in their CO cores. Models with f_CO corresponding to ¹²C(α, γ)¹⁶O rates between |$0\, \sigma$| and |$+3\, \sigma$| become more and more unstable as they evolve. In contrast, models with |$-1\, \sigma$|⁠, |$-2\, \sigma$|⁠, and |$-3\, \sigma$| tend to stabilize and continue their evolution until oxygen depletion, avoiding PI.

Figure 5.

Comparison of the averaged first adiabatic exponent, |$\langle \Gamma _{1} \rangle _{\rm \, }$|versus the central temperature of pure-He stars models with M_ZAMS = 40 M_⊙. The plot shows the final phases at the end of the evolution. Coloured lines indicate models computed with different ¹²C(α, γ)¹⁶O reaction rates. Symbols indicate the start and the end of core burning phases of elements with the same colour code as adopted in Fig. 2. The red horizontal line corresponds to |$\mbox{$\langle \Gamma _{1} \rangle _{\rm }$}=4/3$|⁠. The above blue shaded area indicates values between 4/3 and 4/3 + 0.01, range in which the whole star starts to be dynamically unstable.

Open in new tabDownload slide

3.2 Pre-supernova masses

Table A1 in Appendix  A shows the properties of the pure-He stars computed in this work. The final fate is assigned as follows: when |$\mbox{$\langle \Gamma _{1} \rangle _{\rm TOT-}$}\gt 0$|⁠, the star is stable, continues to evolve through the last burning phases and a final core collapse (CC) will form a black hole with mass ≤M_Pre. When |$\mbox{$\langle \Gamma _{1} \rangle _{\rm TOT-}$}\le {}0$|⁠, the star is dynamically unstable. Depending on its mass, it may start to pulsate losing mass through PPI, or it might explode after a single pulse, becoming a PISN and leaving no remnant. In Table A1, we refer to both PPI and PISN as PI. The evolution through such phases must be followed by means of hydro-dynamical simulations, as done by, e.g. Woosley (2017), Takahashi (2018) and Marchant et al. (2019). Since we did not run hydro-dynamical simulations, we conservatively assume that both stars going through PPI and stars undergoing PISN leave no compact objects. Hence, we will obtain lower limits to the maximum compact object mass.

In Appendix  B, we show a comparison between our pure-He star models and those computed by Farmer et al. (2020) with the mesa stellar evolutionary code (Paxton et al. 2011, 2013, 2015, 2018). From such a comparison we see that our stable models at the end of the computation are very similar to the models of Farmer et al. (2020) that are dynamically stable during the oxygen burning phase, which evolve directly toward the final CC.

Due to both the low metallicity adopted, that leads to mass-loss rates of ≈5 × 10⁻⁶ M_⊙yr⁻¹, and to the short lifetimes of pure-He stars, of the order of 10⁵ yr, all the pure-He models computed have pre-supernova masses, M_Pre, very close to M_He-ZAMS. As discussed before and also shown in Fig. 3, different rates for the ¹²C(α, γ)¹⁶O reaction change the carbon-to-oxygen mass fraction at the end of the He-MS, but do not influence the He and CO core mass.

Fig. 6 shows the lower edge of the PI gap as a function of the ¹²C(α, γ)¹⁶O reaction rate. For pure-He stars, the maximum He-ZAMS mass that results in a stable model and evolves until core collapse is about 68 and 19 M_⊙, in the case of ¹²C(α, γ)¹⁶O reaction rates |$-3\, \sigma$| and |$+3\, \sigma$|⁠, respectively. This trend is similar to the one shown by Farmer et al. (2020), but in our case the maximum mass is always lower by ≈10–20 M_⊙ than theirs, for each assumed rate. This happens because, on the one hand, parsec includes slightly different physical prescriptions for the computation of the evolutionary tracks (e.g. for the opacity, mass-loss, and convection treatment) and, on the other hand, we do not compute the hydro-dynamical evolution of the models after they become unstable. Hence, we conservatively assume that all stars with 〈Γ₁〉_- ≤ 0 leave no compact object, because we cannot distinguish between PPI and PISN and we cannot model mass loss during PPI. Hydro-dynamical simulations of PPI (Marchant et al. 2019; Farmer et al. 2020) show that, depending on the mass of the star, the pulses may eject just few tenths of M_⊙ before the star stabilizes. After this small mass-loss, the star continues its evolution until core collapse. Hence, including hydrodynamical evolution might stabilize some of our models and lead to a higher value of the M_gap. For the same reason, our maximum BH mass obtained adopting the standard ¹²C(α, γ)¹⁶O rate differs by ≈20 M_⊙ with respect to the maximum BH mass obtained from pure-He stars by Woosley (2017), that is 48 M_⊙, and by Leung et al. (2019), 50 M_⊙.

Figure 6.

Pre-Supernova masses of pure-He stars as a function of the ¹²C(α, γ)¹⁶O reaction rate. Circles indicate models that are dynamically stable while burning oxygen. Empty circles: models with a central oxygen mass fraction between 0.3 and 0.2; filled circles: models with a central oxygen mass fraction lower than 0.2. Black crosses indicate models that are dynamically unstable (i.e. 〈Γ₁〉_TOT- < 0). Table A1 lists all the values plotted here. Continuous blue and orange lines indicate the lower and higher mass gap edges by Farmer et al. (2020), respectively. The green star at |$0\, \sigma$| indicate the M_Gap = 48 M_⊙ found by Woosley (2017) for pure-He stars.

Open in new tabDownload slide

The analysis of pure-helium stars, the comparison with Farmer et al. (2020) pure-He models (Appendix  B), and the fact that our results are comparable to previous work confirm that our criterion for stability is a conservative one to decide whether a stellar model is stable, even if we do not follow the hydrodynamical evolution of the final stages. Since we do not follow the evolution through the pulsation PPI phase, our results are lower limits to the M_gap.

4 FULL STELLAR MODELS WITH HYDROGEN ENVELOPE

We now consider stellar models with hydrogen envelopes, computed with the same version of parsec as described in Section 3.1. Grids of evolutionary tracks are calculated with different values of the multiplier parameter, f_CO, as shown in Table 1.

4.1 Evolutionary tracks

Each grid contains stellar models with initial masses of 40, 50, 60, 70, 75, 80, 85, 90, 95, 100, 105, 110, 120, 140, and 160 M_⊙. Our models begin their evolution from the ZAMS and evolve until the end of the central oxygen burning, when possible. The initial metallicity and helium content are Z = 0.0003 and Y = 0.249, respectively (Bressan et al. 2012).

In Fig. 7, we plot the HR diagram of three selected sets of tracks computed with f_CO corresponding to |$+3\, \sigma$|⁠, |$0\, \sigma$|⁠, |$-3\, \sigma$|⁠. The ¹²C(α, γ)¹⁶O reaction does not affect the hydrogen burning phase and the models with the same mass evolve in the same way during the MS.

Figure 7.

HR diagram of the three sets of H stars computed with Z = 0.0003 and with different reaction rates. The colour code of lines and symbols is the same as in Fig. 2. The cyan stars (circles) indicate the beginning (the end) of the hydrogen burning phase. The grey dashed line is the ZAMS. The numbers on the left of the ZAMS are the values of M_ZAMS in units of M_⊙.

Open in new tabDownload slide

After the MS, the stellar cores rapidly contract until their central temperature is high enough for helium ignition. Stars cross the so called Hertzsprung gap while the effective temperatures decrease. Whether a massive star ignites helium as blue supergiant (BSG) or as a yellow/red supergiant (YSG/RSG) star depends on many stellar and physical parameters (such as the stellar mass, the metallicity, the mass-loss, the opacity, and previous mixing efficiency). In our sample, all models with M_ZAMS ≤  90 M_⊙ ignite helium as BSG stars, while more massive ones ignite helium as YSG/RSG stars.

During the CHeB, stars develop a convective core that stops growing when the central He mass fraction falls below ≈0.5. In this phase, the carbon fraction initially increases up to a maximum and then decreases because of the effects of the ¹²C(α, γ)¹⁶O reaction. As the ¹²C(α, γ)¹⁶O rate increases (from |$- 3\, \sigma$| to |$+ 3\, \sigma$|⁠), more carbon is converted into oxygen, eventually leaving an almost carbon-free core at the end of the CHeB phase in the case of the model computed with |$+ 3 \, \sigma$|⁠. However, these differences do not affect much the CHeB lifetimes and the mass of the CO core, as shown in Fig. 8. During the CHeB, the star develops several intermediate convective regions in the envelope, that remain detached from the He core for the remaining evolution. As shown in Fig. 8, the CHeB lifetimes are only slightly affected by the assumption of different f_CO values, with a maximum difference of ≈5 × 10⁴ yr in the case of the M_ZAMS = 40 M_⊙ star. Moreover, adopting different f_CO values has no impact on the CO core masses left after the CHeB phase. Only in stars with M_ZAMS ≥ 120 M_⊙ the ¹²C(α, γ)¹⁶O rate does affect the CO mass, because of the presence of an efficient envelope undershooting, as it will be further discussed in Section 4.3. Fig. 8 shows that the major effect of changing f_CO is on the central carbon-to-oxygen ratio left at the end of the CHeB, as in pure-He stars.

Figure 8.

Same as Fig. 3, but for stars with hydrogen envelopes. The upper panel shows the MS lifetimes of hydrogen stars versus the ZAMS mass. The second panel from the top shows the CHeB lifetimes. The third panel shows the CO core mass. Different colours indicate different ¹²C(α, γ)¹⁶O rates (several of these lines overlap). Table A2 lists all the relevant quantities. The lower panel shows the carbon-to-oxygen ratio in the CO cores at the end of CHeB.

Open in new tabDownload slide

At such low metallicity, all our models evolve through advanced burning phases as RSG stars. This is due to the low stellar winds that allow stars to keep their massive H-rich envelope.

After the CHeB phase, the core starts to contract again increasing its density and temperature. Meanwhile, neutrino losses tend to cool down the core, but, when the central temperature reaches ≈1 GK, carbon is ignited. Stars computed with lower ¹²C(α, γ)¹⁶O rates have longer lifetimes than models computed with higher ¹²C(α, γ)¹⁶O rates.

After carbon depletion, neon photo-disintegration begins at the stellar centre when the temperature rises up to about 1.2–1.5 GK. During carbon and neon burning, our models have radiative cores and convective envelopes. As already seen for pure-He stars, during the central neon burning phase, carbon burns off-centre in a shell that, in case of low ¹²C(α, γ)¹⁶O rates, is able to turn on an intermediate convective region.

This convective zone inside the CO core continuously provides fresh carbon to the burning shell that sustains the overlying layers, thus delaying the ignition of oxygen in the core. As the evolution proceeds further, 〈Γ₁〉_TOT may become lower than 4/3 + 0.01, indicating that the star becomes globally unstable. In this case, the computation is stopped. Otherwise, we continue the evolution until the end of central oxygen burning, when eventually the code encounters some convergence problems. Apart from a few isolated cases, the values of the final central oxygen mass fractions in stable models are between 0.01 and 0.3 (Table A2).

After the central oxygen burning phase, the evolution of the core and of the envelope become decoupled. Therefore, in analogy with the pure-He models, we can assume that if the star remains dynamically stable during the core oxygen burning, it will evolve towards the advanced silicon and iron burning phases and finally, to the core collapse.

Fig. 9 shows the Kippenhahn diagrams of stars with M_ZAMS = 100 M_⊙ computed with f_CO corresponding to ¹²C(α, γ)¹⁶O rate |$+ 3\, \sigma$|⁠, |$0 \, \sigma$|⁠, and |$- 3\, \sigma$|⁠. The three models evolve in a similar way during the MS and the CHeB phases. They build up a convective core that decreases over time in the MS. At the end of the MS, stars have an He core of ≈50 M_⊙. During the CHeB, the models develop a convective core of ≈47 M_⊙, in which all the helium is converted to carbon, oxygen, and neon.

Figure 9.

Upper panels: Same as Fig. 4, but for hydrogen stars with M_ZAMS = 100 M_⊙. The continuous grey lines indicate the total stellar mass. The cyan-coloured areas indicate the hydrogen burning regions, which contribute at least for the 1 per cent of the L_H. Lower panels: Luminosity as a function of time from H-burning, He-burning, C-burning, and neutrinos in cyan, orange, green, and black, respectively.

Open in new tabDownload slide

When temperature and density are high enough, carbon is ignited in the core. It burns for 12 yr, 2.8 yr, and only 4 months in the case of models computed with f_CO corresponding to |$- 3\, \sigma$|⁠, 0 and |$+ 3 \, \sigma$|⁠, respectively. After the carbon burning, neon photo-disintegration takes place in the star centre, and it lasts for about 20, 2, and 1.4 days for models computed with |$- 3\, \sigma$|⁠, 0, and |$+ 3\, \sigma$|⁠, respectively. At this point, only the model computed with |$- 3\, \sigma$| remains globally stable and burns central oxygen in a small convective core of about 2 M_⊙. At the end of the computation, this 2 M_⊙ core is composed of 15 per cent ¹⁶O, 45 per cent ²⁸Si, 30 per cent ³²S, and 10 per cent heavier elements. This star will continue its evolution through the most advanced burning phases and then will collapse forming a BH. The other models computed with |$0\, \sigma$| and |$+ 3\, \sigma$| do not have sufficient ¹²C in the CO core to prevent PI, hence they become unstable before the ignition of oxygen.

Fig. 10 shows the evolution of the averaged first adiabatic exponent, 〈Γ₁〉_TOT, of the models with 100 M_⊙ and computed with different ¹²C(α, γ)¹⁶O rates. Depending on the assumed rate, the star enters or not in the unstable regime. In this case, only models with |$-3\, \sigma$| and |$-2 \, \sigma$| remain stable during the oxygen burning phase. The other models become unstable shortly after the end of the neon burning phase.

Figure 10.

Comparison of the 〈Γ₁〉_TOT values versus the central temperature of stars with hydrogen envelopes and M_ZAMS = 100 M_⊙. Here, we use the same colour code for lines and symbols as in Fig. 5. The lines corresponding to +2 and |$+3\, \sigma$| are barely visible, because they almost perfectly overlap with the |$+1\, \sigma$| line.

Open in new tabDownload slide

4.2 Pre-supernova masses

Table A2 in Appendix  A shows the results for the stars with H envelopes. The columns list the same quantities as shown in Table A1, with the exception of the first adiabatic exponent (defined in Section 2.5), that here is computed in two ways: (i) for whole star, 〈Γ₁〉_TOT- = 〈Γ₁〉_TOT - (4/3 + 0.01), and (ii) for the He core, 〈Γ₁〉_Core- = 〈Γ₁〉_Core - (4/3 + 0.01). We find that 〈Γ₁〉_Core is always slightly lower than 〈Γ₁〉_TOT. The envelope of massive stars could be near the dynamical instability too, depending on the ionization state of the plasma (Stothers 1999). 〈Γ₁〉_Core and 〈Γ₁〉_TOT are not significantly different, because 〈Γ₁〉 is a weighted integral and the external parts of the envelope contribute less than the internal parts of the star. We use the same methodology adopted for the pure-He stars to assign the final fate of stars. In the case of 〈Γ₁〉_TOT- > 0, the star is considered globally stable and evolves until the final CC. On the other hand, if 〈Γ₁〉_TOT- ≤ 0, the star becomes dynamically unstable and the PI-induced collapse leads to a PPI or a PISN. This fate is indicated with the label PI in the Table. We find that the final fate does not change if we assume 〈Γ₁〉_TOT or 〈Γ₁〉_Core.

Fig. 11 shows the pre-supernova masses of stars with H envelope as a function of the ¹²C(α, γ)¹⁶O rate. The lower edge of the mass gap (M_gap) ranges from ≈60 to ≈92 M_⊙ in the case of models computed with |$+ 3\, \sigma$| and |$-2\, \sigma$|⁠, respectively. Interestingly, we find no mass gap in the case of models computed with |$- 3\, \sigma$|⁠. The maximum M_Pre is ≈150 M_⊙. Moreover, we find that stars computed with |$-2\, \sigma$| and with initial masses of 140 and 160 M_⊙ do not become unstable and burn oxygen non-explosively, with final M_Pre masses of ≈131 and ≈150 M_⊙ respectively.

Figure 11.

Pre-supernova masses as a function of the ¹²C(α, γ)¹⁶O reaction rate for stars with H envelopes. Circles indicate models that are dynamically stable while burning oxygen. Red empty circles: Models with a central oxygen mass fraction higher than 0.3. Black empty circles: Models with a central oxygen mass fraction between 0.3 and 0.2. Filled circles: Models with a central oxygen mass fraction lower than 0.2. Black crosses indicate models that are dynamically unstable (i.e. 〈Γ₁〉_TOT- < 0). The PI mass gap is completely filled in the |$-3\, \sigma$| case, because of the effect of the dredge-up (see Section 4.2 for details). Continuous blue and orange lines indicate the lower and higher mass gap edges by Farmer et al. (2020), respectively. The green star at |$0\, \sigma$| indicates the M_gap = 65 M_⊙ found by Woosley (2017) for stars with hydrogen envelope.

Open in new tabDownload slide

The disappearance of the mass gap for models with |$-3\, \sigma$| is due to a dredge-up during the CHeB, that reduces the mass of the CO core with the effect of stabilizing the core. This process is described in detail below.

In this case, the lower edge of the PI window is higher than that found by Farmer et al. (2020), because we consider stars that still have their hydrogen envelopes. The range of maximum masses we find is in agreement with the maximum BH mass ≈65 M_⊙, found by Woosley (2017) for stars with hydrogen envelope, and overlaps with the range obtained by Takahashi (2018), that is between 100 and 270 M_⊙.

4.3 Envelope undershooting and dredge-up during the advanced phases

The dredge-up is an opacity-driven process that begins when the star approaches the location of RSG. The envelope expands and cools, the opacity rises, and the radiative temperature gradient becomes superadiabatic in the inner envelope layers. As a consequence, the base of the convective envelope extends inwards, crossing the H–He discontinuity, and deepening into the He core (envelope undershooting).

As described in Section 3.1, in parsec we adopt the undershooting at the bottom of the outermost convective region. As a consequence, only models that develop fully convective envelopes, extending from the surface to the He core, undergo the dredge-up.

Models with M_ZAMS < 120 M_⊙, ignite and deplete helium as BSG/YSG stars (Fig. 7), regardless of the ¹²C(α, γ)¹⁶O reaction rate assumed. Such models show a complex structure with many intermediate convective regions: they do not develop fully convective envelopes and avoid the dredge-up (see Fig. 9).

In contrast, all our models with M_ZAMS ≥ 120 M_⊙ become RSG stars after the MS, regardless of their ¹²C(α, γ)¹⁶O rate. They develop a large convective envelope and undergo a dredge-up after the ignition of helium, which extracts a few M_⊙ of material from the core (Fig. 12).

Figure 12.

Left-hand panels: Same as Fig. 9, but the evolution is zoomed on the CHeB phase of the M_ZAMS = 140 M_⊙ star computed with the ¹²C(α, γ)¹⁶O reaction rate |$+3\, \sigma$|⁠. Time until the final collapse is in linear scale. In the lower panel, the dotted lines show the contribution to the He luminosity from the CHeB phase main reactions, the triple-α, the ¹²C(α, γ)¹⁶O, and the ¹⁶O(α, γ)²⁰Ne in purple, red, and green, respectively. Right-hand panels: Same plots as those shown in the left-hand panels, but for a model with M_ZAMS = 140 M_⊙ and |$- 3\, \sigma$|⁠.

Open in new tabDownload slide

In the following stages, models with high ¹²C(α, γ)¹⁶O rates have a very different evolution from models with low ¹²C(α, γ)¹⁶O rates (Fig. 12). Models with high ¹²C(α, γ)¹⁶O rates tend to perform blue loops, as in the case of less massive stars (Iben 1974; Bertelli, Bressan & Chiosi 1985; Brunish & Becker 1990). For example, the 140 M_⊙ model with |$+3\, \sigma$|  ¹²C(α, γ)¹⁶O rate begins the blue loop when the ¹²C(α, γ)¹⁶O luminosity exceeds the triple-α one. At higher effective temperatures, the H-burning shell re-ignites and the convective envelope becomes thinner, while some detached intermediate convective zones appear (left-hand panels of Fig. 12). Envelope undershooting is not applied to such intermediate convective regions. Thus, the penetration into the He-core is not boosted even if, sometimes, convection extends down to the bottom of the H-rich envelope. We will investigate these cases in a follow-up study.

In contrast, models computed with low ¹²C(α, γ)¹⁶O reaction rates deplete He in the RSG branch, maintaining their large convective envelopes. Hence, the dredge-up process continues (right-hand panels of Fig. 12), leading to a significant decrease of the He core mass during CHeB. This keeps these models out of the PI regime.

Fig. 13 shows the complete Kippenhahn diagram and the luminosity evolution until the end of the O-burning of a star with M_ZAMS = 140 M_⊙ and ¹²C(α, γ)¹⁶O rate |$-3\, \sigma$|⁠. During the MS, the star has a big convective core, which decreases over time, and a radiative envelope. As the core recedes, small intermediate convective zones form in the envelope. At the end of the MS, the star moves towards the red part of the HR-diagram expanding and cooling. At the beginning of the CHeB, the star develops a large convective region that, after about 5 × 10⁴ yr, extends from the stellar surface to the bottom of the H-rich envelope (upper-right panel of Fig. 12). During this phase, the H-burning shell quenches out and the undershooting at the base of the envelope penetrates in the He core bringing helium up to the surface (dredge-up).

Figure 13.

Same as Fig. 9, but for the M_ZAMS = 140 M_⊙ star with |$- 3\, \sigma$|⁠.

Open in new tabDownload slide

During CHeB phase, more than ≈20 M_⊙ of helium are extracted from the core (Table A2). The temperature at the bottom of the H-burning shell (T_b) reaches 40 Mk. At the end of this stage, the envelope is He-rich (with an helium mass fraction of 0.52) and has super-solar metallicity (Z_CNO ≈ 0.13).

After CHeB, the core contracts and the H-shell ignites again, being the only nuclear energy source for a while. The excess luminosity causes a further quick penetration of the external convection corresponding to a second dredge-up. This brings the H discontinuity very close to the CO core leaving only a very thin He shell. In the meantime, the core contracts and both neutrino and hydrogen luminosities rise. The latter causes a further expansion of the external envelope accompanied by a continuous convective envelope undershooting into the CO core. T_b reaches about 70 MK. The undershooting mixing and subsequent hot bottom burning continue until carbon ignites in the centre, increasing appreciably the surface metallicity. At this stage, the convective envelope has a mass of ≈105 M_⊙ and a composition of 3.5 per cent ¹H, 23 per cent ⁴He, 48 per cent ¹⁴N, 23 per cent ¹⁶O, and 2.5 per cent of ²⁰Ne and heavier elements. At the end of the central carbon burning phase, lasting for 5.7 × 10² yr, the CO core mass is ≈24.8 M_⊙. Subsequent central neon burning lasts for about 2 months, while carbon burns in a shell above the core.

At this point, the star, which is well outside the PI region, evolves through the advanced burning phases until core collapse. This is shown in Fig. 14 where the upper panel shows the 〈Γ₁〉_TOT value as a function of the central temperature, while the lower one shows the evolution of the central density and temperature. These plots show models with M_ZAMS = 140 M_⊙ computed with different ¹²C(α, γ)¹⁶O reaction rates. The branching moment between models happens just after the beginning of the CHeB, when stars with rate |$-3\, \sigma$| and |$-2\, \sigma$| undergo the dredge-up, which erodes significantly the He core. The model with |$-3\, \sigma$| has a second branching, when the dredge-up starts to extract matter from the CO core. The |$-2\, \sigma$| model does not undergo this second branching. In all the other cases, the star evolves becoming unstable during or shortly after the neon burning phase in the core.

Figure 14.

Upper panel: 〈Γ₁〉 values versus the central temperature of stars with M_ZAMS = 140 M_⊙. Lower panel: Central density and temperature evolution of models with an initial mass of 140 M_⊙. In both the upper and the lower right-hand panels, we use the same colour code for lines and symbols as in Figs 5 and 10.

Open in new tabDownload slide

Table A2 shows if a star undergoes the dredge-up during its evolution. It is worth to note that not all the stars that experience the dredge-up avoid the PI.

4.4 Comparison with previous works

Takahashi (2018) found that very massive stars avoid pulsational PI if they develop a carbon burning shell in the core. This happens if they assume low reaction rates.

Clarkson & Herwig (2020) studied the interaction between convective H shell and He core in massive stars. They found that H shell and He core interact in models adopting the Schwarzschild criterion, while they do not if the Ledoux criterion is assumed. Models with shell interaction can produce up to 1000× more ¹⁴N in the shell with respect to models that do not have a shell interaction.

Kaiser et al. (2020) studied the role of overshooting and semi-convection in stars with masses between 15 and 25 M_⊙, adopting different criteria for convection. They found relative uncertainties up to about 70 per cent on the core masses and lifetimes depending on the mixing assumption. Moreover, the relative importance of semi-convection decreases with an increasing amount of overshooting. The amount of overshooting is the main source of uncertainty in all phases.

Recently, Farrell et al. (2020) discuss the role of the H-shell–He-core interaction in metal-poor massive stars. They find that uncertainties related to convective mixing, mass-loss, H–He interactions, and PPI may increase the M_gap up to ≈85 M_⊙.

In our case, the dredge-up is driven by envelope undershooting. If we turn-off the undershooting, by assuming a Λ_env = 0 H_P, the dredge-up disappears and all models with M_ZAMS = 140 M_⊙ become globally unstable and undergo PI (as shown in Fig. 15). The fact that other authors (e.g. Woosley 2017; Takahashi 2018; Marchant et al. 2019; Farmer et al. 2020) do not observe such process in their models is due to different assumptions in the treatment of convection. For instance, stars evolved adopting the Schwarzschild criterion and without semi-convection tend to develop bigger convective regions with respect to stars computed with the Ledoux criterion, due to molecular weight gradient.

Figure 15.

Left-hand panel: Estimated mass of the BH (M_BH) as a function of the ZAMS mass of the progenitor star. Different colours indicate stars computed with varying ¹²C(α, γ)¹⁶O reaction rates. Solid lines indicate the total final mass of the star. Dotted lines indicate the final He core mass. As discussed in the text, the actual M_BH mass will be between the He core mass (pessimistic case, dotted lines) and the total final mass of the star (optimistic case, solid lines). The green area indicates models that experience dredge-up during the CHeB. Right-hand panel: Same as the left-hand panel but for a set of tracks computed without envelope undershooting.

Open in new tabDownload slide

4.5 BH mass spectrum

The left-hand panel of Fig. 15 shows the BH mass as a function of M_ZAMS and for different ¹²C(α, γ)¹⁶O reaction rates. Since we do not know what fraction of the hydrogen envelope is lost during the final collapse, we give a pessimistic and an optimistic estimate for the BH mass in Fig. 15, which bracket this uncertainty. In the pessimistic case (dotted line), the mass of the BH is the same as the final He core mass, which corresponds to assuming that the H envelope is completely lost during the collapse. In the optimistic case (solid line), the mass of the BH is the same as the entire final mass of the star, including all the residual H envelope.

For stars computed with f_CO corresponding to |$+ 3\, \sigma$|⁠, |$+ 2\, \sigma$|⁠, and |$+ 1\, \sigma$|⁠, we find a maximum BH mass M_BH = 60 M_⊙. For stars computed with |$0 \, \sigma$|⁠, we find a maximum BH mass of 68 M_⊙.

The PI mass gap for models computed with |$- 1\, \sigma$| is between ≈80 and ≈150 M_⊙. In fact, the dredge-up allows stars with ZAMS mass M_ZAMS = 160 M_⊙ to leave BHs with mass ≈150 M_⊙.

The effect of dredge-up becomes stronger if we consider lower rates for the ¹²C(α, γ)¹⁶O reaction. In the |$-2\, \sigma$| case, the PI mass gap is between M_BH = 92 M_⊙ and 110 M_⊙. Finally, models computed with |$- 3\, \sigma$| do not show any mass gap. The mass gap is completely removed by the effect of dredge-up.

The right-hand panel of Fig. 15 shows the same plot of the estimated BH mass but for stars evolved without envelope undershooting (Λ_env = 0 H_P). In this case, stars more massive than 110 M_⊙ do not undergo the dredge-up during the CHeB phase. Without the dredge-up, the mass gap re-appears, even for low ¹²C(α, γ)¹⁶O rates. The maximum stellar mass that ignites oxygen non-explosively is about 112 M_⊙, in the case with rate |$- 3\, \sigma$|⁠. Hence, if undershooting is suppressed, the lower edge of the mass gap M_gap is at 60 M_⊙ for rates +3, +2, and +1, while M_gap ≈ 68, 80, 92, and 112 M_⊙ for rates 0, −1, −2, and |$-3\, \sigma$|⁠, respectively.

Table 3 lists the results obtained for the lower edge of the mass gap for stars with H-rich envelopes, computed with and without envelope undershooting and adopting the optimistic case.

5 DISCUSSION

The recent detection of GW190521 (Abbott et al. 2020b,e), with a primary BH mass of |$85^{+21}_{-14}$| M_⊙ and a secondary BH mass of |$66^{+17}_{-18}$| M_⊙ (90 per cent credible intervals), challenges current models of BH formation.

Here, we have shown that the lower edge of the PI mass gap, M_gap, can be as high as ≈70 M_⊙, if we assume that the H envelope collapses to a BH directly (Fig. 15). Moreover, if we not only assume that the H envelope collapses but we also take into account the uncertainty on the ¹²C(α, γ)¹⁶O nuclear reaction rate, M_gap can be even higher. For the standard ¹²C(α, γ)¹⁶O rate |$-1\, {}\sigma$|⁠, the lower edge and the upper edge of the mass gap are ≈80 and ≈150 M_⊙, respectively. For the standard ¹²C(α, γ)¹⁶O rate |$-2\, {}\sigma {}$|⁠, the lower edge and the upper edge of the mass gap are ≈92 and ≈110 M_⊙, respectively. Finally for the standard ¹²C(α, γ)¹⁶O rate |$-3\, {}\sigma$|⁠, the PI mass gap completely disappears, because of the envelope undershooting effect we just discussed.

The main reasons why previous work has neglected the effect of the H envelope, calculating only models of pure-He stars, is that the H envelope is usually lost during the evolution of an interacting binary star, because of Roche lobe mass transfer and common-envelope episodes. Only metal-poor single stars, or metal-poor stars in loose binary systems (with initial semi-major axis ≳10³ R_⊙, see e.g. Fig. 9 of Spera et al. 2019) can preserve most of their H envelope to the very end of their life.

Loose binary systems cannot lead to BBH mergers via isolated binary evolution, because the initial orbital separation of the BBH is too large to permit coalescence by GW emission. Hence, even in our most optimistic case, it is hard to form systems like GW190521 through isolated binary evolution.

However, in a dense stellar system, the evolution of a single BH (formed from the collapse of a single massive star) and that of a loose BBH (formed from the evolution of a loose binary star) can be very different from that of isolated binary systems (e.g. Sigurdsson & Phinney 1995; Portegies Zwart & McMillan 2000; Banerjee, Baumgardt & Kroupa 2010; Downing et al. 2010; Tanikawa 2013; Samsing, MacLeod & Ramirez-Ruiz 2014; Ziosi et al. 2014; Morscher et al. 2015; Mapelli 2016; Rodriguez, Chatterjee & Rasio 2016; Banerjee 2017, 2020; Samsing 2018; Fragione & Kocsis 2018; Arca Sedda, Askar & Giersz 2018; Kumamoto, Fujii & Tanikawa 2019, 2020; Di Carlo et al. 2019, 2020a, b; Fragione, Loeb & Rasio 2020; Rizzuto et al. 2020). The single BH can pair up with another BH via three-body encounters and dynamical exchanges. A loose BBH can harden (i.e. its semi-major axis shrinks) via three body encounters, speeding up its merger by gravitational waves. The more massive a BH is, the more it is effective in acquiring companions via exchanges and in undergoing dynamical hardening (Heggie 1975; Hills 1976; Hills & Fullerton 1980). For this reason, dynamics favours the formation and merger of the most massive BBHs. Thus, the formation of a BH with mass  85 M_⊙ from the collapse of a single star can easily lead to a system like GW190521 in a dense stellar system. In particular, Di Carlo et al. (2020b) have already shown that the dynamical formation of a merging BBH with masses similar to GW190521 (88 and 48 M_⊙) is possible in a young dense star cluster.

Moreover, a ≈85 M_⊙ BH born from the direct collapse of a star will preserve most (if not all) the spin of its progenitor star, leading to a fast rotating BH (Mapelli et al. 2020). If this BH acquires a companion via dynamical interactions, the spins of the primary and of the secondary BH will be isotropically oriented, possibly leading to a large precession spin (χ_p). This feature is consistent with the properties of GW190521 (Abbott et al. 2020b,e).

Recently, Belczynski (2020) even suggested that isolated binary star evolution at low metallicity can produce BBHs with masses consistent with GW190521, when the uncertainties on the ¹²C(α, γ)¹⁶O rate are taken into account (assuming −2.5 σ rates from Farmer et al. 2020) and when extreme stellar masses (M_ZAMS ≈ 180 M_⊙) are considered. However, the isolated binary evolution channel predicts a lower χ_p than the one estimated by Abbott et al. (2020e) within the 90 per cent credible interval.

6 CONCLUSIONS

We investigated the lower edge of the pair-instability (PI) BH mass gap (M_gap), by means of the parsec evolutionary code.

We evolved (i) a set of pure-He stars with initial masses between 20 and 100 M_⊙ and with metallicity Z = 0.001, (ii) a set of stars with H envelopes, with masses ranging from 40 to 160 M_⊙and with Z = 0.0003. For both sets, we have varied the ¹²C(α, γ)¹⁶O rates from |$-3\, {}\sigma {}$| to |$+3\, \sigma {}$|⁠, where |$1\, \sigma$| is the given confidence error by the starlib library. We have followed the evolution of the models until oxygen burning, if the star remains globally stable. We find that, changing the ¹²C(α, γ)¹⁶O rates directly affects the carbon-to-oxygen ratio at the end of the CHeB.

In the case of pure-He stars, the lower edge of the mass gap (M_gap) varies from 68 M_⊙ (⁠|$-3\, \sigma {}$|  ¹²C(α, γ)¹⁶O rate) to 20 M_⊙ (⁠|$+3\, \sigma$| rate). For the standard ¹²C(α, γ)¹⁶O rate, we find M_gap ≈ 30 M_⊙.

Since we use a hydro-static code, we do not follow the evolution through the PPI phase. We conservatively assume that both stars going through PPI and stars going through PISN leave no compact objects. Since stars going through PPI might be able to stabilize and form BHs by direct collapse (Marchant et al. 2019), our findings are robust lower limits for the lower edge of the mass gap. For this reason, we find lower values of M_gap with respect to previous works (Woosley 2017; Leung et al. 2019; Farmer et al. 2020) in case of pure-He stars.

In the case of stellar models with an hydrogen envelope, the main uncertainty on the final BH mass is represented by the fate of the envelope (Mapelli et al. 2020). Here, we consider a pessimistic estimate in which only the He core collapses to the final BH and an optimistic approach in which the H envelope entirely collapses to the final BH.

Adopting the standard ¹²C(α, γ)¹⁶O rate and assuming an optimistic approach, we find M_gap = 68 M_⊙. Assuming the |$+3\, \sigma$|  ¹²C(α, γ)¹⁶O rate, this value lowers to M_gap = 59 M_⊙.

If we consider lower values than the standard ¹²C(α, γ)¹⁶O rate, we find that the interaction between the hydrogen envelope and the He core may result in a very different evolution of the star. Models computed with low ¹²C(α, γ)¹⁶O rates (−1, −2, |$-3\, \sigma$|⁠) and with an initial mass higher that 110 M_⊙ may experience an important dredge-up during the CHeB that extracts matter from the core enriching the envelope.

Stars experiencing such dredge-up are more stable and ignite oxygen in a non-explosive way, despite their high mass. The most remarkable result is that the PI mass gap narrows in the −1 and |$-2\, \sigma$| cases, and it completely disappears in the |$-3\, \sigma$| case.

In the |$-3\, \sigma$| case, all our models evolve until the core collapse, without undergoing PI. In case of |$-2\, \sigma$| and |$-1\, \sigma$| we find a mass gap between 92–110 M_⊙ and 80–150 M_⊙, respectively.

For higher rates (0, +1, +2, |$+3\, \sigma {}$|⁠), the dredge-up does not happen and the lower edge of the mass gap is M_gap ≈ 60–70 M_⊙ in the optimistic case (M_gap ≈ 30–40 M_⊙ in the pessimistic case).

We confirm that the ¹²C(α, γ)¹⁶O reaction has a strong impact on the final fate of a massive star, as found by other authors (Takahashi 2018; Farmer et al. 2019, 2020). In addition, we show that another challenging phenomenon could be the occurrence of an efficient dredge-up at the bottom of the H-rich envelope, in the advanced phases. Stars that experience such a dredge-up may develop envelopes heavily enriched by N and O. In the case such enriched envelopes are ejected, their contribution to stellar yields, especially N, may be significant. This could be relevant for the chemical evolution of nitrogen at low metallicity (Pettini et al. 2008; Vincenzo et al. 2016).

Taking into account the uncertainties on the collapse of the H envelope and on the ¹²C(α, γ)¹⁶O reaction rate, we can explain the formation of massive BHs such as the primary component of GW190521 (Abbott et al. 2020b,e) with the collapse of a single massive metal-poor star. If such massive BH forms in a dense stellar cluster, it might have a chance to dynamically pair up with another massive BH, leading to the formation of a GW190521-like system, with a large precession spin.

ACKNOWLEDGEMENTS

We are grateful to the anonymous referee for their useful comments. We thank Rob Farmer for the fruitful discussions and for providing us mesa pure-He models. MM, GC, and GI acknowledge financial support from the European Research Council for the ERC Consolidator grant DEMOBLACK, under contract no. 770017. AB acknowledges support from PRIN-MIUR 2017 prot. 20173ML3WW. This research made use of NumPy (Harris et al. 2020), SciPy (Virtanen et al. 2020), IPython (Perez & Granger 2007). For the plots, we used Matplotlib, a Python library for publication of quality graphics (Hunter 2007).

DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author.

Footnotes

REFERENCES

APPENDIX A: TABLES

Tables A1 and A2 show the results obtained for pure-He and full stars. More details in Sections 3.2 and 4.2.