The impact of pair-instability mass loss on the binary black hole mass distribution

TL;DR

This paper explores how pair-instability supernovae influence the maximum mass of black holes in binary systems, using population synthesis models to match observed gravitational wave data and predict future observational signatures.

Contribution

It introduces a population synthesis model incorporating PPISNe and PISNe effects, predicting a maximum black hole mass of 40 solar masses and estimating merger fractions involving these phenomena.

Findings

Maximum black hole mass predicted is 40 M_

Approximately 10% of mergers involve a PPISN-processed black hole

Predicted rates of PPISNe and PISNe increase with redshift up to z=2

Abstract

A population of binary black hole mergers has now been observed in gravitational waves by Advanced LIGO and Virgo. The masses of these black holes appear to show evidence for a pile-up between $30$--$45$ $M_\odot$ and a cut-off above $\sim 45$ $M_\odot$. One possible explanation for such a pile-up and subsequent cut-off are pulsational pair-instability supernovae (PPISNe) and pair-instability supernovae (PISNe) in massive stars. We investigate the plausibility of this explanation in the context of isolated massive binaries. We study a population of massive binaries using the rapid population synthesis software COMPAS, incorporating models for PPISNe and PISNe. Our models predict a maximum black hole mass of $40$ $M_\odot$. We expect $\sim 10$\% of all binary black hole mergers at redshift z = 0 will include at least one component that went through a PPISN (with mass $30$--$40$ $M_\odot$), constituting $\sim 20$--$50$\% of binary black hole mergers observed during the first two observing runs of Advanced LIGO and Virgo. Empirical models based on fitting the gravitational-wave mass measurements to a combination of a power law and a Gaussian find a fraction too large to be associated with PPISNe in our models. The rates of PPISNe and PISNe track the low metallicity star formation rate, increasing out to redshift $z = 2$. These predictions may be tested both with future gravitational-wave observations and with observations of superluminous supernovae.

Tables

Table 1. Table 1 Minimum and maximum helium core masses to undergo PPISN assumed in the various models

Model	MHe,min	MHe,max
Belczynski+ 2016	45	60
Woosley 2017	35	65
Marchant+ 2018	35	60
Woosley 2019	30	60

Table 2. Table 2 Coefficients for Equation  4

Coefficient	Value
$c_0$	7.39643451 $\times {10}^3$
$c_1$	-1.13694590 $\times {10}^3$
$c_2$	7.45060098 $\times {10}^1$
$c_3$	-2.69801221 $\times {10}^0$
$c_4$	5.83107626 $\times {10}^{-2}$
$c_5$	-7.52206933 $\times {10}^{-4}$
$c_6$	5.36316755 $\times {10}^{-6}$
$c_7$	-1.63057326 $\times {10}^{-8}$

Table 3. Table 3 Description of parameter names used in this paper

Parameter	Description
${ℛ}_{\mathrm{BBH}}$	BBH Merger Rate in Gpc-3 yr-1 at $z=0$
${ℛ}_{\mathrm{BBH}\mathrm{with}\mathrm{PPISN}}$	PPISN Formed BBH Merger Rate in Gpc-3yr-1 at $z=0$
RBBH	Predicted observed rate of BBH mergers during O1 and O2 per year
R${}_{\mathrm{BBH}\mathrm{with}\mathrm{PPISN}}$	Predicted observed rate of PPISN formed BBH mergers during O1 and O2 per year
${\Lambda }_{\mathrm{merge}}$	Proportion of BBHs merging at $z=0$ formed via PPISNe
${\Lambda }_{\mathrm{obs}}$	Proportion of observed BBHs with at least one BH formed from a progenitor which underwent PPISN
$M_{\max }$	Max mass of PPISN-formed BH (M⊙)
$\beta$	Max remnant mass formed not via PPISN (M⊙)
$Z_{\mathrm{PPISN}}$	Highest metallicity PPISN formation occurred
$Z_{\mathrm{PISN}}$	Highest metallicity PISN occurred

Table 4. Table 4 Summary of results for each of our models described in section  2.5 . The uncertainties in the quoted rates are the 1 σ 𝜎 \sigma Monte-Carlo statistical uncertainties only. Descriptions of the parameters and their units are given in Table  3 .

Parameter	Woosley 2017	Belczynski+ 2016	Marchant+ 2018	Woosley 2019	IMF-2.1	IMF-2.5
${ℛ}_{\mathrm{BBH}}$	40 $\pm 4$	40 $\pm 4$	40 $\pm 4$	40 $\pm 4$	60 $\pm 6$	20 $\pm 2$
${ℛ}_{\mathrm{BBH}\mathrm{with}\mathrm{PPISN}}$	3 $\pm 1$	2 $\pm 1$	3 $\pm 1$	8 $\pm 2$	6 $\pm 3$	2 $\pm 1$
RBBH	20	20	20	20	30	10
R${}_{\mathrm{BBH}\mathrm{with}\mathrm{PPISN}}$	7	3	5	8	10	3
${\Lambda }_{\mathrm{merge}}$	0.1	0.02	0.1	0.2	0.1	0.1
${\Lambda }_{\mathrm{obs}}$	0.4	0.2	0.4	0.5	0.4	0.3
$M_{\max }$	58.4	40.5	39.5	39.5	58.4	58.4
$\beta$	31.4	40.2	31.5	26.9	31.4	31.4
$Z_{\mathrm{PPISN}}$	0.005	0.002	0.006	0.007	0.005	0.005
$Z_{\mathrm{PISN}}$	0.001	0.001	0.002	0.002	0.001	0.001

Equations

$$\zeta =\frac{d\ln R}{d\ln m}.$$

$$M_{\mathrm{final}}=0.9\times 45{\mathrm{M}}_{\odot }.$$

$$M_{\mathrm{final}}=0.9\times M_{\mathrm{He}}.$$

$$M_{\mathrm{final}}/M_{\mathrm{He}}=\sum _{\ell =0}^7{c}_{\ell }\left(\frac{M_{\mathrm{He}}^{\ell }}{M_{\odot }}\right).$$

$${\mathcal{R}}_{\mathrm{PPISN}}(z)=\frac{{\mathrm{d}}^2{N}_{\mathrm{PPISNe}}}{\mathrm{d}{t}_{\mathrm{s}}\mathrm{d}{V}_{\mathrm{c}}}=\int \mathrm{d}Z\int \mathrm{d}\tau$$

$$\left[\frac{{\mathrm{d}}^2{N}_{\mathrm{PPISNe}}}{\mathrm{d}{M}_{\mathrm{form}}\mathrm{d}\tau }(Z)\frac{{\mathrm{d}}^3{M}_{\mathrm{form}}}{\mathrm{d}{t}_s\mathrm{d}{V}_c\mathrm{d}Z}\left(Z,t_{\mathrm{form}}\left(\tau ,z\right)\right)\right]$$

$$\psi (z)=0.01\frac{(1+z{)}^{2.77}}{1+((1+z)/2.9{)}^{4.7}}{\mathrm{M}}_{\odot }{\mathrm{Mpc}}^{-3}{\mathrm{yr}}^{-1},$$

$$⟨Z⟩=⟨Z_0⟩\times 10^{-0.23z},$$

$$\mathcal{M}=(m_1{m}_2{)}^{3/5}/(m_1+m_2{)}^{1/5}.$$

$$\frac{{\mathrm{d}}^2{N}_{\mathrm{BBH}}}{\mathrm{d}{t}_{\mathrm{obs}}\mathrm{d}z}=\frac{{\mathrm{d}}^2{N}_{\mathrm{BBH}}}{\mathrm{d}{t}_{\mathrm{s}}\mathrm{d}{V}_{\mathrm{c}}}\frac{\mathrm{d}{V}_{\mathrm{c}}}{\mathrm{d}z}\frac{\mathrm{d}{t}_{\mathrm{s}}}{\mathrm{d}{t}_{\mathrm{obs}}}{P}_{\det },$$

Full text

1. Introduction

Since beginning observing runs in 2015, Advanced LIGO (aLIGO; Aasi et al., 2015) and Advanced Virgo (AdVirgo; Acernese et al., 2015) have confirmed the detection of 11 gravitational wave (GW) events to date (Abbott et al., 2016b; The LIGO Scientific Collaboration et al., 2018b). The sources of 10 of the detections are binary black hole (BBH) mergers, and the source of the other event was the merger of a binary neutron star (BNS) system (Abbott et al., 2017; The LIGO Scientific Collaboration et al., 2018b). With the planned advancements to current GW detectors around the world, the rate of GW detections is set to vastly increase over the coming years (Abbott et al., 2018). GW observations will allow us to constrain astrophysical models (e.g. Stevenson et al., 2015; Zevin et al., 2017; Wysocki et al., 2018a).

The mass distribution of BBHs is commonly modelled as a power law (e.g. Fishbach & Holz, 2017; Kovetz et al., 2017; Talbot & Thrane, 2018; Wysocki et al., 2018b; Roulet & Zaldarriaga, 2019; The LIGO Scientific Collaboration et al., 2018a), up to some maximum BH mass. With only 10 observed BBHs, the fine details of the BBH mass distribution remain uncertain at present (The LIGO Scientific Collaboration et al., 2018a). However, even with a small number of detections, there is some evidence for certain features in the mass distribution.

One of the most prominent features to emerge so far in the BBH mass distribution is a lack of BHs more massive than $\sim 45$ M⊙ (Fishbach & Holz, 2017; Talbot & Thrane, 2018; Wysocki et al., 2018b; Roulet & Zaldarriaga, 2019; The LIGO Scientific Collaboration et al., 2018a), despite ground-based GW detectors being sensitive to such systems. Currently, the most massive BH, observed in the event GW170729, had a mass of $50.6^{+16.6}_{-10.2}$ M⊙ (The LIGO Scientific Collaboration et al., 2018b).

From the theory of stellar evolution, there is a prediction of an absence of BHs with masses $\sim 50$–$150$ M⊙. Massive stars with helium core masses in the range $\sim 50$–$150$ M⊙ are believed to become unstable to electron-positron pair production (Fowler & Hoyle, 1964; Barkat et al., 1967; Fraley, 1968). This causes the radiation pressure support in the core to drop, causing the core to contract. As it contracts, the temperature increases, triggering explosive oxygen burning. This may reverse the contraction and completely unbind the star in a pair-instability supernova (PISN) explosion, leaving no remnant behind (e.g. Kozyreva et al., 2014a, b, 2017). BH formation is expected again above a helium core mass of $\sim 150$ M⊙ (Woosley et al., 2002).

There is also tentative evidence for an excess of BHs with masses in the range $30$–$45$ M⊙ (The LIGO Scientific Collaboration et al., 2018a). Theory suggests that stars with pre-supernova helium core masses in the range $\sim 30$–$50$ M⊙ also experience the pair-instability, but the release of energy is insufficient to completely disrupt the star. These stars may instead experience pair-instability induced pulsations multiple times (e.g. Woosley, 2017; Yoshida et al., 2016; Marchant et al., 2018). Each pulsation ejects material in a supernova-like event, leading to them being typically called pulsational pair-instability supernovae (PPISNe). After every pulse, the star’s structure converges back to equilibrium. An iron core forms in equilibrium, and then collapses in a regular core-collapse supernova (CCSN; Woosley et al., 2007). These pulsations remove the envelope of the star, rather than completely unbind the star. Enhanced mass loss during the pulsations may cause all stars in this mass range to form BH remnants with masses in a narrow range.

To date there have been no confirmed observations of PPISNe or PISNe, although some fraction of superluminous supernovae (SLSN-R) are potential candidates (Woosley et al., 2007; Gal-Yam et al., 2009; Quimby et al., 2011; Cooke et al., 2012; Gal-Yam, 2012; Lunnan et al., 2016; Kozyreva et al., 2018; Gomez et al., 2019). The supernova iPTF2014hls (Arcavi et al., 2017) may potentially be a PPISN (Woosley, 2018), possibly from a stellar merger product (Vigna-Gómez et al., 2019). The supernova iPTF16eh (Lunnan et al., 2018) also shows evidence for a massive shell of material ejected $\sim 30$ years prior to explosion, consistent with a PPISN.

There are several possible alternate explanations for the apparent excess of BHs in the mass range $30$–$45$ M⊙. Due to the relatively small number of observed BBH mergers, the statistical uncertainties in the inferred BBH mass distribution are still large, and a pure power law mass distribution is only mildly disfavoured (The LIGO Scientific Collaboration et al., 2018a). It is also likely that the true mass distribution of BBHs in nature is not a perfect power law, and the preference for a model with more degrees of freedom than a power law may simply be highlighting this fact.

Alternatively, finding clusters of GW observations in different parts of parameter space (e.g. a cluster of massive BHs) may indicate that multiple formation channels are contributing to the observed population (e.g. Powell et al., 2019).

BBHs are thought to form through several evolutionary channels (Mandel & Farmer, 2018). All channels in which the BHs are of stellar origin will be subject to the effects of PISNe and PPISNe described above. These include classic isolated binary evolution (e.g. Belczynski et al., 2016a; Stevenson et al., 2017; Kruckow et al., 2018; Spera et al., 2019), dynamical formation in dense stellar environments such as open clusters (Rastello et al., 2019; Chatterjee et al., 2017), globular or nuclear clusters (Kulkarni et al., 1993; Sigurdsson & Hernquist, 1993; Rodriguez et al., 2016; Fragione & Kocsis, 2018), active galactic nuclei (Bartos et al., 2017; Stone et al., 2017), formation in triple (Antonini et al., 2017; Rodriguez & Antonini, 2018) or quadruple systems (Fragione & Kocsis, 2019), through wide binaries in the field (Michaely & Perets, 2019) or through quasi-chemically homogeneous evolution in close tidally locked binaries (Mandel & de Mink, 2016; Marchant et al., 2016; Eldridge & Stanway, 2016).

Although PISNe are still likely to operate in dense stellar environments such as globular clusters (Rodriguez et al., 2018), imposing a maximum BH mass, it may be possible to form more massive BHs in these environments. For several of the observed BBH mergers (e.g. GW150914, GW170729), the merger product would be a BH in the PISN mass gap (The LIGO Scientific Collaboration et al., 2018b). It is possible that such so called second generation BHs could again participate in mergers, filling this gap (O’Leary et al., 2016; Samsing & Ilan, 2019; Fishbach et al., 2017; Gerosa & Berti, 2017; Rodriguez et al., 2018; Kimball et al., 2019). Black holes formed from stellar merger products in star clusters may also populate the PISN mass gap (e.g. Di Carlo et al., 2019).

More exotic scenarios for the origin of the currently observed BBHs include primordial BHs (Bird et al., 2016), or several of the events being a single gravitationally lensed BBH merger (Broadhurst et al., 2019). The effects of PISNe are likely not relevant to these scenarios.

In this paper we focus on BBHs formed through classical isolated binary evolution. We introduce prescriptions for modelling the effects of both PISNe and PPISNe on massive stars in our rapid binary evolution code COMPAS (Stevenson et al., 2017; Vigna-Gómez et al., 2018; Barrett et al., 2018; Neijssel et al., 2019). We calculate the expected distribution of BH masses, and show that we expect to observe a maximum BH mass of around $\sim 40$ M⊙. Our models also show a mild excess of BHs in the $30$–$45$ M⊙ range due to PPISNe, with around $10$% of BBHs merging at redshift $z=0$ having at least one component which has undergone PPISN. Due to GW selection effects favouring more massive systems, this corresponds to around $20$–$50$% of observed BBH mergers having at least one BH which has undergone a PPISN. This supports the hypothesis that the more massive black holes in the BBH mergers GW150914 and GW170729 may have formed in this way. We also calculate the volumetric rate of PISNe, PPISNe, and BBH mergers , showing that they increase from redshift $z=0$ to $z=2$, tracking the low metallicity cosmic star formation rate.

In Section 2, we give a brief description of the population synthesis code and set-up used in this study. In Section 3 we calculate the volumetric rate of various astrophysical phenomena including PISNe, PPISNe and BBH mergers as a function of redshift. We also present the intrinsic mass distributions of BBHs merging at redshift $z=0$. In Section 4 we present our predicted observed BBH mass distributions, incorporating GW selection effects. We discuss our results, compare with the literature and highlight limitations of this study in Section 5. We conclude in Section 6.

2. Binary population model

We use the population synthesis software COMPAS (Stevenson et al., 2017; Vigna-Gómez et al., 2018; Barrett et al., 2018; Neijssel et al., 2019) to study a population of massive isolated binaries. COMPAS models the evolution of binary systems, incorporating approximate prescriptions for stellar evolution (Tout et al., 1996; Hurley et al., 2000), stellar mass loss (Vink et al., 2001; Belczynski et al., 2010), mass transfer and common envelope evolution (see section 2.2), supernovae (see sections 2.4 & 2.5) and GW emission (Peters, 1964). Our model neglects the effects of stellar rotation and tides.

2.1. Initial binary distributions

In COMPAS, binaries are Monte Carlo sampled from probability distributions of initial conditions based on astrophysical observations. We draw the mass of the initially more massive star from an initial mass function (IMF; Kroupa, 2001) in the mass range $5$–$150$ M⊙, whilst the secondary mass is determined from the mass ratio, drawn from a uniform distribution (Sana et al., 2012), with a lower limit of 0.1 M⊙, since stars with masses less than $\sim 0.08$ M⊙ do not burn hydrogen (Kumar, 1963; Hayashi & Nakano, 1963). We account for the low mass binaries and single stars we do not simulate in section 3.1. We discuss the uncertainties inherent in our population synthesis approach in Section 5.

We draw the initial orbital separation of the binary from a uniform in log distribution with limits of 0.1 and 1000 AU (Sana et al., 2012). We assume all binaries are initially circular.

We assume that these initial distributions are independent. Moe & Di Stefano (2017) recently presented a distribution of correlated initial parameters and binary fractions based on a compilation of observations. Uncertainties in the initial distributions of massive binaries typically translate to uncertainties in predictions of BBH rates at the factor of 2 level and do not dominate, with the exception of the initial mass function (de Mink & Belczynski, 2015; Klencki et al., 2018). We explicitly demonstrate the uncertainty in our predictions due to uncertainty in the initial mass function in section 3.

We make the simplifying assumption, in the absence of compelling evidence to the contrary, that these initial distributions do not depend on metallicity or redshift.

We use a grid of 50 metallicities distributed uniformly in log space between $Z_{\mathrm{min}}=1\times 10^{-4}$ and $Z_{\mathrm{max}}=0.03$, set by the limits of the underlying stellar models (Pols et al., 1998; Hurley et al., 2000). At each metallicity we model the evolution of $2\times 10^{5}$ binary systems, for a total of $1\times 10^{7}$ binaries.

In the following subsections we briefly summarise some of the most relevant details for the present study.

2.2. Mass Transfer Stability and Common Envelope Evolution

Stars expand during their lifetimes. If the stellar radius exceeds the Roche-radius of the star it will transfer mass to its companion star. The stability of this mass transfer depends on the response of both the stellar radius and the Roche-radius of the star. This stability is often defined in a radius-mass relationship (see for example Paczyński & Sienkiewicz (1972), Hjellming & Webbink (1987))

[TABLE]

The mass transfer is dynamically stable if $\zeta_{*}>\zeta_{RL}$, where $\zeta_{*}$ is the adiabatic response of the star, and $\zeta_{RL}$ the response of the Roche-lobe. A summary of our prescriptions for the stellar response at different evolutionary stages can be found in section 2.2.4 of Vigna-Gómez et al. (2018). If the mass transfer is stable we evolve the system in a similar manner to Hurley et al. (2002) (see also Stevenson et al. (2017)).

When the mass transfer is dynamically unstable the system will start a common-envelope event (Paczynski, 1976; Ivanova et al., 2013). In a common-envelope event, the envelope of the donor star engulfs the entire binary. In order to determine whether the binary survives the common-envelope event or results in a stellar merger, we use the “$\alpha$-$\lambda$” formalism (Webbink, 1984; de Kool, 1990; Ivanova et al., 2013). This formalism compares the orbital energy to the energy needed to unbind the envelope. We assume that all the orbital energy goes into expelling the envelope ($\alpha_{\mathrm{CE}}=1$). Our parametrization for the binding energy of the envelope (characterised by $\lambda$) is determined from the fitting formulae of Xu & Li (2010), as in StarTrack (Dominik et al., 2012). A requirement is that we know both the core and envelope mass of the donor star. Massive stars crossing the Hertzsprung gap may not have developed a clear core-envelope separation (Ivanova et al., 2013). If this is not the case, common-envelope events from these donors could always result in a stellar merger. This uncertainty can have a large impact on predictions for the rate of merging binary black holes (up to an order of magnitude c.f. Dominik et al., 2012). We assume that all common-envelope events with Hertzsprung-gap donors result in stellar mergers (Belczynski et al., 2007).

2.3. Typical evolutionary channel for BBHs

There are multiple possible channels for forming a BBH from a massive isolated binary. Here we briefly describe the typical evolutionary stages involved. At present, COMPAS does not model the formation of BBHs in close, tidally locked binaries through chemically homogeneous evolution (Mandel & de Mink, 2016; Marchant et al., 2016); work is underway to incorporate this into our model.

The dominant channels for forming BBHs in our models are detailed in Stevenson et al. (2017) and Neijssel et al. (2019). We provide a brief summary here. A typical BBH is formed from an initially wide binary with an orbital period of a few hundred days, consisting of two massive O-type stars with zero-age main-sequence (ZAMS) masses greater than $\sim 20$ M⊙. As the initially more massive star (the primary) evolves off the main sequence first, its radius increases. As it does, the primary fills its Roche lobe (Eggleton, 1983) and begins to transfer mass to its companion, until all of its envelope is removed, leaving a helium star. The helium star may undergo PPISNe before finally collapsing to a BH. The initially less massive secondary star then evolves off the main sequence and begins transferring mass onto the BH. This can either be stable mass transfer (e.g. Pavlovskii et al., 2017; van den Heuvel et al., 2017; Neijssel et al., 2019), or can result in the binary undergoing common envelope evolution (e.g. Dominik et al., 2012; Belczynski et al., 2016a; Stevenson et al., 2017; Marchant et al., 2018; Neijssel et al., 2019); both dramatically shrink the orbital period of the binary and removing the hydrogen envelope of the secondary star (see Section 2.2 for more details). Finally, the secondary star may undergo PPISNe before collapsing to a BH, forming a BBH. The orbit of the BBH then slowly shrinks due to emission of GWs (Peters, 1964) over a time period of up to the age of the universe.

2.4. Supernovae

Massive stars end their lives in supernovae, resulting in either the formation of a neutron star, a BH, or the complete destruction of the star. From radio observations of isolated Galactic pulsars, neutron stars are known to receive velocity kicks at birth of up to several hundred km s*-1* (Hobbs et al., 2005). Our assumptions for neutron star natal kicks are described in Section 2.2.3 of Vigna-Gómez et al. (2018). It is unclear both observationally and theoretically whether BHs also receive such kicks (Repetto et al., 2012; Janka, 2013; Mandel, 2016; Repetto et al., 2017). For canonical CCSNe we use the ‘delayed’ supernova engine from Fryer et al. (2012) to determine the remnant masses, with kicks drawn from a Maxwellian with a 1D dispersion of 265 km s*-1*, as for neutron stars (Hobbs et al., 2005; Vigna-Gómez et al., 2018). These kicks are then suppressed by the fraction of mass falling back on to the BH, as described in Fryer et al. (2012). In practice, this means that most BHs formed in our models receive no natal kick.

2.5. PPISN & PISN Models

In order to investigate the occurrence rates and properties of PPISNe and PISNe, we need to be able to relate the final properties of the remnant to those of its progenitor. This involves several assumptions and is very model dependent. In this section, we therefore describe several different models for determining the remnant mass $M_{\mathrm{final}}$ for BHs formed from progenitors earlier undergoing PPISN as a function of their pre-collapse helium core mass $M_{\mathrm{He}}$. These models are based on Belczynski et al. (2016b) (section 2.5.1), Woosley (2017) (section 2.5.2), Marchant et al. (2018) (section 2.5.3) and Woosley (2019) (section 2.5.4). Other studies of PPISNe in the literature include Leung et al. (2019), which finds similar results to Marchant et al. (2018) and Woosley (2017), and Moriya & Langer (2015) and Yoshida et al. (2016), which do not present a dense enough grid of models to be useful for the present study.

For each model, we present a simple analytic fit to the results of the more detailed models. We apply this fit between the limiting values for the helium core mass given in Table 1. We visualise both the original data and our fits in Figure 1. In all models, we assume that stars with helium cores more massive than $M_{\mathrm{He,max}}$ and less massive than 135 M⊙ (Woosley, 2017) undergo a PISN and leave no remnant. Cores more massive than 135 M⊙ are assumed to directly collapse to BHs (Woosley et al., 2002; Woosley, 2017; Fryer et al., 2012). In this paper, we only study stars with ZAMS masses below 150 $M_{\odot}$, and so we never form such massive cores. We leave an analysis of BBHs with component masses above the PISN mass gap to future work. See section 5 for a discussion of issues with stellar models at such high masses. We do not explicitly model the multiple mass loss kicks (Blaauw, 1961; Boersma, 1961) from the multiple PPISNe mass loss episodes (Woosley, 2017; Marchant et al., 2018); since the mass lost in most episodes is small, we do not expect this to significantly impact our results.

2.5.1 Belczynski et al. 2016

Belczynski et al. (2016b) assume that all stars with a helium core more massive than 45 M⊙ lose all of their exterior mass through PPISN. They give the final BH mass as

[TABLE]

This prescription is not a particularly accurate approximation to the more detailed models discussed below. We include it here to demonstrate the effect of using this model on our predictions. This model has been adopted by other authors investigating the BH mass distribution (e.g. Rodriguez et al., 2018).

2.5.2 Woosley 2017

Woosley (2017) simulates a grid of pure helium stars without mass loss with masses between 30 and 64 M⊙. They also evolve a grid of hydrogen-rich stars with masses between 70 and 150 M⊙ and metallicity $Z=0.1\,Z_{\odot}$, including mass loss. We perform a linear regression to fit the remnant BH mass as a function of the pre-supernova helium core mass as

[TABLE]

This fit does not attempt to capture any possible turnover in the relation at around 60 M⊙.

2.5.3 Marchant et al. 2018

Marchant et al. (2018) simulate a grid of single, non-rotating, pure helium stars with metallicity $Z=Z_{\odot}/10$ with masses between 40 and 100 M⊙ using MESA (Paxton et al., 2011, 2018) including mass loss. These models demonstrate a turnover in the relation between pre-supernova helium core mass and final mass (see Figure 1). We find that a 7th order polynomial is a good fit to the data of Marchant et al. (2018), giving

[TABLE]

The $c_{\ell}$ coefficients are given in Table 2. We bound the ratio between 0 and 1. We use this model as our reference model for Figures 2, 3, 4 & 7.

2.5.4 Woosley 2019

Woosley (2019) also studies PPISNe in massive helium stars, including mass loss. The remnant masses calculated by Woosley (2019) agree reasonably well with the models of Marchant et al. (2018), although the former find a lower minium helium core mass which undergoes PPISN of $\sim 30$ M⊙. Figure 1 shows that the models of Woosley (2019) are in broad agreement with those of Marchant et al. (2018). We therefore choose to use the same polynomial fitting formula given by Equation 4, with different values for $M_{\mathrm{He,min}}$ and $M_{\mathrm{He,max}}$ given in Table 1.

2.6. Caveats

For our models based on both Marchant et al. (2018) and Woosley (2019) we apply an additional factor of 0.9 to the masses to account for the difference between the baryonic and gravitational masses of the final BH, to be consistent with supernova engine adopted from Fryer et al. (2012). This avoids introducing an artificial mass gap at a mass $M_{\mathrm{He,min}}$ that we do not believe to be physical. However, the 10% mass loss for BHs in Fryer et al. (2012) is ad hoc and may be an overestimate111If the peak neutrino luminosity is $\sim 10^{53}$ erg s*-1*, and the BH formation occurs on a timescale of $\sim 1$ s (Mirizzi et al., 2016; Chan et al., 2018) (substantially shorter than the Kelvin-Helmholtz timescale of the proto-neutron star), then the maximum energy loss can be a few $\times 10^{53}$ erg = 0.1 M⊙ c2. (Bernhard Müller, private communication, 2019). This will impact our estimates of the location of the maximum BH mass by up to $\sim 5$ M⊙.

Since the detailed PPISN models we use in Section 2.5 are only calculated for a single metallicity, we assume that the relation between pre-supernova helium core mass and BH mass is independent of metallicity. This can be updated in the future if significant metallicity effects are discovered in the detailed models.

In our models, we assume that the natal kick received by a BH formed from a progenitor which has undergone PPISNe is the same as one which did not. In our Fryer et al. (2012) natal kick prescription, this means that PPISN progenitors do not receive a natal kick in practice. We do not model orbital changes from individual pulsations.

3. Results

Massive stars lose a substantial fraction of their mass through line driven winds (Castor et al., 1975) and eruptive mass loss (Humphreys & Davidson, 1994). Line driven winds are thought to be quenched in low metallicity environments (Vink et al., 2001). It is therefore expected that it is possible to form more massive helium cores, and thus more massive BHs, in lower metallicity environments (e.g. Belczynski et al., 2010; Spera et al., 2015). We show in Figure 2 that BBHs form at all metallicities $Z\lesssim 0.02$ in our models. Both PPISN and PISN occur down to our lowest metallicity of $Z=1\times 10^{-4}$, and we expect that they would continue to occur at lower metallicities (e.g. Fryer et al., 2001; Yoon et al., 2012; Chatzopoulos & Wheeler, 2012). The highest metallicities at which PISNe and PPISNe occur in our models are around $Z\sim 0.002$ and $Z\sim 0.006$ respectively222We find a very small number of PISN and PPISN events at high metallicity $Z\gtrsim 0.01$ due to very particular mass transfer histories. Although we include these systems in computing rates, the maximum metallicities quoted in Table 4 exclude these rare systems (see Figure 2 and Table 4 for details). These are in reasonable agreement with values found in previous studies for PISNe (Belczynski et al., 2016b; Spera & Mapelli, 2017). However, Spera & Mapelli (2017) find that PPISN occur in their models up to metallicities of $Z\sim 0.018$ (see their Fig. 3), much higher than the value we find. This discrepancy is likely due to different assumed stellar wind mass loss prescriptions, which are highly uncertain (Renzo et al., 2017).

3.1. Merger and Supernovae Rates - Methods

We calculate the volumetric rates of 7 different astrophysical phenomena:

the rate of stars experiencing PPISNe 2. 2.

the rate of PISNe 3. 3.

the total rate of CCSNe including those stars which went through PPISNe before collapsing 4. 4.

formation rate of BBHs 5. 5.

formation rate of BBHs in which at least one of the component BHs formed from a progenitor that underwent PPISN 6. 6.

BBH merger rate 7. 7.

rate of merging BBHs in which at least one of the component BHs formed from a progenitor that underwent PPISN.

Each rate is calculated in a similar way. For example, the volumetric rate of PPISNe $\mathcal{R}_{\mathrm{PPISN}}$ at redshift $z$ is given by

[TABLE]

where $Z$ is the metallicity and $t_{\mathrm{form}}$ is the age of the universe when a given source would need to have formed if there is delay $\tau$ between the time of formation and the event of interest.

The first fraction in the square brackets in Equation 5 is the number of events (PPISNe) formed per unit star-forming mass per unit delay time at a metallicity Z. We show the number of events per unit star-forming mass as a function of metallicity in Figure 2. In agreement with other authors (e.g. Giacobbo et al., 2018; Spera et al., 2019; Neijssel et al., 2019) we find that the formation efficiency of BBHs is a strong function of metallicity, decreasing with increasing metallicity. At $Z=1\times 10^{-3}$, the formation efficiency of BBHs is $\sim 5\times 10^{-5}$ M⊙-1. We also find that the formation efficiency of PISNe is a strong function of metallicity, decreasing from $\sim 10^{-5}$ M⊙-1 at $Z=1\times 10^{-4}$ to $\sim 10^{-6}$ M⊙-1 at $Z=1\times 10^{-3}$.

Since we did not include single stars333Our widest binaries have separations of up to $1000$ AU, and so do not interact or low mass binaries (see Section 2.1), we must account for the mass we did not simulate. Using a Monte Carlo simulation with COMPAS over the full mass range 0.1–150 M⊙ and assuming a binary fraction of 70% (Sana et al., 2012) independent of primary mass, we find that we must multiply the total mass of our simulation by a factor of $\sim 5$ in order to account for the stars we did not simulate. Since we do not include single stars, our calculated rates are only for those events occurring in binaries. This means that the rates of supernovae in our model should be interpreted as lower limits. In addition, our adopted binary fraction is appropriate for the massive stars we are interested in, but we note that the binary fraction decreases with decreasing mass (e.g. Raghavan et al., 2010; Moe & Di Stefano, 2017). However, we do not expect this simplification to have a large impact on our results, as we focus on the outcome of the most massive stars’ lives (c.f. Klencki et al., 2018).

The second fraction in the square brackets in Equation 5 is the metallicity-specific star formation rate. Since we are interested in transient events with a strong metallicity dependence (as shown in Figure 2), we need to know not only the total star formation rate, but the amount of star formation occurring at a given metallicity. Both the total star formation rate, and the fraction of star formation occurring at a given metallicity evolve with time (or equivalently redshift) throughout the history of the universe (Madau & Dickinson, 2014).

We assume an analytic, phenomenological model for the metallicity specific star formation rate introduced by Neijssel et al. (2019), calibrated to match GW observations. The star formation rate at redshift $z$ is

[TABLE]

where we take the star formation rate at redshift $z=0$ from Madau & Fragos (2017). We assume that metallicities are log-normally distributed with a spread of 0.39 dex about a mean metallicity $\langle Z\rangle$ which scales with redshift $z$ as

[TABLE]

where $\langle Z_{0}\rangle=0.035$ is the mean metallicity of the universe at redshift $z=0$.

We convert between lookback times and redshifts using the standard expressions (e.g. Hogg, 1999)444We use a flat $\Lambda$-CDM cosmological model, using the WMAP9 cosmological parameters (Hinshaw et al., 2013) as implemented in astropy’s cosmology module. We convert between redshift and lookback time using astropy’s lookback_time.

Since the typical lifetime of massive stars is only a few Myr, we neglect any time between binary formation and supernovae ($\tau=0$), and assume that these occur at the same redshift (e.g. Zapartas et al., 2017). The PISN rate was calculated using the same method as above, with the substitution of $\mathrm{d}N_{\mathrm{PISN}}$ for $\mathrm{d}N_{\mathrm{PPISNe}}$ in Equation 5. We do the same when calculating the BBH formation rate, substituting instead $\mathrm{d}N_{\mathrm{BBH}}$.

When calculating the BBH merger rate, one must take more care, since the delay times $\tau$ between binary formation and GW driven merger can be long (Gyrs; e.g. Dominik et al., 2012). We therefore use $t_{\mathrm{form}}=t_{\mathrm{merge}}(z)-\tau$, where $t_{\mathrm{merge}}(z)$ is the age of the universe when the BHs merge (Barrett et al., 2018). Binaries with a merger time longer than the age of the universe at a given redshift are excluded from the analysis. We also calculate both the formation and merger rates of BBHs where at least one of the components was formed from a progenitor which underwent a PPISN.

3.2. Merger and Supernovae Rates - Results

We show the rates of various astrophysical phenomena in Figure 3. At redshift $z=0$, the total CCSN rate is $\sim 2\times 10^{4}$ Gpc*-3* yr*-1*.

Both PISNe and PPISNe are extremely rare events. At redshift $z=0$, we find the volumetric rate of PISNe to be less than $10^{-2}$ Gpc*-3* yr*-1* and the rate of PPISNe to be $\sim 10^{-1}$ Gpc*-3* yr*-1*. These events are rare at low redshift in our model because we assume that the mean metallicity of the local universe is super-solar (c.f. Equation 7).

At higher redshifts, where the mean metallicity of the universe is lower, the rates of these events increase. We find the rate of PPISNe at redshift $z=4$ is $\sim 2\times 10^{3}$ Gpc*-3* yr*-1* and the rate of PISNe at redshift $z=4$ is $\sim 10^{-1}$ Gpc*-3* yr*-1*. As we discuss in section 5, our rate of PPISNe quoted here should be considered a lower limit.

We find that including PPISNe does not change the overall BBH merger rate appreciably compared to models not including PPISNe, in agreement with Belczynski et al. (2016b). The BBH merger rate at redshift $z=0$ in our models is $\mathcal{R}_{\mathrm{BBH}}\sim 40$ Gpc*-3* yr*-1* (see Table 4), in agreement with empirical estimates from GW observations of $10$–$100$ Gpc*-3* yr*-1*(The LIGO Scientific Collaboration et al., 2018a). The total BBH merger rate is sensitive to the metallicity-specific star formation history of the Universe (see Neijssel et al. 2019 for a detailed discussion) and additional assumptions about binary evolution (see Section 5 for more details). Therefore, we henceforth focus on the relative contribution of PPISN remnants as a fraction of all BBH mergers.

The rates of BBH mergers, PISNe and PPISNe all follow the low-metallicity star formation rate, initially increasing with redshift and peaking around redshift $z=2$ before decreasing to higher redshifts. This is due to a combination of both the overall star formation rate increasing, as well as the average metallicity decreasing between $z=0$ and $z=2$ (Madau & Dickinson, 2014).

The BBH merger rate at a given redshift is typically an order of magnitude smaller than the BBH formation rate at the same redshift in our models, as a large fraction of BBHs have delay times longer than the age of the universe and never merge (see also Spera et al., 2019). The BBH merger rate peaks at later times (lower redshifts) than the BBH formation rate due to the significant time delays between formation and merger.

We define $\Lambda_{\mathrm{merge}}$ as the fraction of BBHs merging at a given redshift in which at least one BH was formed from a progenitor that went through a PPISN. We similarly define $\Lambda_{\mathrm{form}}$ to be the fraction of BBHs forming at a given redshift in which at least one BH was formed from a progenitor that went through a PPISN. We show the value for both $\Lambda_{\mathrm{merge}}$ and $\Lambda_{\mathrm{form}}$ as a function of redshift in Figure 4. In our standard model based upon Marchant et al. (2018), around $\sim 10$% of BBHs merging at redshift $z=0$ have a component BH formed from a progenitor which underwent PPISN (see Table 4 for the value for other models).

Our value for $\Lambda_{\mathrm{merge}}$ is towards the low end of the broad, empirically determined value $\Lambda_{\mathrm{LVC}}=0.4^{+0.3}_{-0.3}$ (The LIGO Scientific Collaboration et al., 2018a, also shown in Figure 4). Our predicted rate of merging BBHs in which one of the BHs formed from a progenitor that went through PPISN is $\sim 3$ Gpc*-3* yr*-1*, which is similarly lower than the rate of the high mass component inferred from GW observations of $10$–$40$ Gpc*-3* yr*-1* (The LIGO Scientific Collaboration et al., 2018a)55590% credible interval obtained from publicly available posterior samples https://dcc.ligo.org/LIGO-P1800324/public.

This result suggests that standard binary evolution models may struggle to produce such a large fraction of BBHs through PPISNe, and another interpretation for the excess of massive black holes may be warranted.

Our rates are sensitive to the assumed IMF (e.g. de Mink & Belczynski, 2015), which is observationally uncertain (e.g. Schneider et al., 2018; Farr & Mandel, 2018). To demonstrate this, we also include two models where we change the power law slope of the IMF for masses greater than 0.5 M⊙ from our default value of $-2.3$ (Kroupa, 2001) to $-2.1$ or $-2.5$, using our PPISN model based on Marchant et al. (2018) (Equation 4). We find that varying the IMF in this range leads to a factor of 2 uncertainty in the total BBH merger rate at redshift $z=0$ (see also de Mink & Belczynski, 2015). The fraction of BBH mergers at redshift $z=0$ in which at least one component underwent a PPISN is relatively robust to uncertainties in the IMF (see Table 4).

3.3. Maximum BH mass

The maximum BH mass $m_{\mathrm{max}}$ for each of our models is given in Table 4. Our models based on Belczynski et al. (2016b), Marchant et al. (2018) and Woosley (2019) have an upper mass limit of $\sim 40$ M⊙ which is consistent with the maximum BH mass inferred from GW observations (Fishbach & Holz, 2017; Talbot & Thrane, 2018; Wysocki et al., 2018b; Roulet & Zaldarriaga, 2019; The LIGO Scientific Collaboration et al., 2018a).

Our fit to the models of Woosley (2017) results in a maximum BH mass of $\sim 53$ M⊙. This higher limit can be explained by the fact that our fit assumes that the remnant mass increases linearly with increasing helium core mass, resulting in larger remnant masses $>50$ M⊙ close to maximum helium core mass $M_{\mathrm{He,max}}$. This can be seen in Figure 1.

Excluding the Belczynski et al. (2016b) model, the most massive BHs to be formed from progenitors without undergoing a PPISN ($\beta$ in Table 4) are in the range of $\sim 27-32$ M⊙ for all our models, with exact numbers shown in Table 4. In these models, BHs with masses in the range $35$–$45$ M⊙ are exclusively formed from progenitors that previously underwent PPISN events.

3.4. BBH mass distribution

We show the intrinsic BBH mass distribution in Figures 5 and 6. Figure 5 shows the intrinsic distribution of the mass of the more massive BH in a merging BBH at redshift $z=0$ for 4 of our models. The sharp spike at $\sim 41$ M⊙ in the Belczynski et al. (2016b) model is consistent with the features of the Belczynski et al. (2016b) fit shown in Figure 1, where a flattening of the predicted remnant mass occurs at $\sim 41$ M⊙.

The distribution is not consistent with a power law in mass across the entire mass range. Our model populations based on Woosley (2017), Marchant et al. (2018) and Woosley (2019) produce BH mass distributions (for the more massive BH) which rise from 5 to $\sim 15$ M⊙, and do not resemble a power law distribution above $\sim 15$ M⊙. Above $\sim 34$ M⊙ there is an increase in the amount of BBHs in the $35-40$ M⊙ region peaking at a final remnant mass of $\sim 39$ M⊙.

A well measured combination of the component masses in a BBH is the chirp mass

[TABLE]

We show the intrinsic chirp mass distributions predicted by our models of BBHs merging at redshift $z=0$ in Figure 6. Our models show that even BBHs with chirp masses as low as $\sim 10$ M⊙ may have one component BH formed from a progenitor which underwent a PPISN (in a binary with a lower mass BH). We discuss the predicted observed mass distributions in Section 4. The models presented in this section can be used to construct phenomenological models of the BBH mass distribution and provide physically motivated priors for their parameters.

4. Observed binary black hole mass distribution

GW observations have strong selection effects favouring more massive systems; the maximum distance a BBH merger is observable to is approximately proportional to $\mathcal{M}^{5/6}$ (for low masses), meaning the observable volume is proportional to $\mathcal{M}^{5/2}$ (Stevenson et al., 2015). This means that the observed distribution is expected to over-emphasize the high-mass tail, allowing it to be explored more readily.

Our method for incorporating GW selection effects is described in Barrett et al. (2018). We assume a representative strain sensitivity for the aLIGO detectors during O1 and O2666We assume each detector has a strain sensitivity equivalent to that of the Hanford LIGO detector during the observation of GW150914, publicly available at https://dcc.ligo.org/LIGO-P1500238/public (Abbott et al., 2016b, a). We have verified that assuming a different representative sensitivity from O2, publicly available at https://dcc.ligo.org/LIGO-P1900011/public, does not significantly change our results.. We assume a single interferometer detection threshold signal-to-noise ratio of 8. We calculate the predicted observed rate as

[TABLE]

where the first term is given by Equation 5, the second term is the differential comoving volume at redshift $z$ (e.g. Hogg, 1999)777We use astropy’s differential_comoving_volume function., $t_{\mathrm{s}}$ is time in the source frame and $t_{\mathrm{obs}}=(1+z)t_{\mathrm{s}}$ is time in the observer’s frame. $P_{\mathrm{det}}$ is the detection probability for a binary with a given set of masses and a given distance, averaged over sky locations and source orientations.

We show the predicted observed redshift distribution in Figure 7. We see that in our model most BBH mergers during O1 and O2 would have been expected to be observed at redshift $z\sim 0.1$.

Our models predict a total observed rate of $\sim 20$ BBHs per year of observing time, integrated over redshift. Since O1 and O2 combined constitute around 0.5 yr of observing time (The LIGO Scientific Collaboration et al., 2018b), our model predicts $\sim 10$ BBH detections during O1 and O2, in agreement with observations.

In Figure 8 we show the predicted observed distribution of the more massive black hole. We define $\Lambda_{\mathrm{obs}}$ as the fraction of observed BBHs in which at least one of the BHs was formed from a progenitor which experienced a PPISN, similarly to Section 3. We find that $\Lambda_{\mathrm{obs}}=0.2$–$0.5$ in our models (see Table 4). This suggests that we would expect $2$–$5$ of the observed BBH mergers to be such systems, strengthening our hypothesis that the more massive black holes in both GW150914 and GW170729 may have formed this way.

We also show the predicted observed chirp mass distributions for our models in Figure 9. Marchant et al. (2018) use their models, along with a toy model for the BH mass distribution, to investigate the expected observed BH chirp mass distribution. Accounting for GW selection effects, their model predicts that the distribution of observed chirp masses will show a double peak structure. We do not see such a feature in our chirp mass distributions shown in Figure 9.

5. Discussion and Caveats

5.1. Results Interpretation

BBHs merging at redshift $z=0$ where at least one BH formed from a progenitor that underwent a PPISN are intrinsically rare, with a volumetric rate of $\sim 3$ Gpc*-3* yr*-1*(see Figure 3). Meanwhile, the merger rate of massive BBHs has been empirically determined through GW observations to be $\sim 10$–$40$ Gpc*-3* yr*-1* (The LIGO Scientific Collaboration et al., 2018a). In addition to the absolute rate, we also present the fraction $\Lambda_{\mathrm{merge}}$ of BBHs in which at least one BH formed from a progenitor that underwent a PPISN. In our models, $\Lambda_{\mathrm{merge}}\sim 0.1$ at redshift $z=0$, meaning that around 10% of merging BBHs form in this way (see Figure 4). This value is towards the low end of the empirically estimated value $\Lambda_{\mathrm{LVC}}=0.4^{+0.3}_{-0.3}$ (The LIGO Scientific Collaboration et al., 2018a).

Our models suggest that it is unlikely that BBHs formed from progenitors which underwent PPISNe can contribute more significantly than this to the population merging at redshift $z=0$. We suggest that an alternative explanation may be required, such as a model for the distribution of the mass of the more massive black hole that is not a pure power law across the entire mass range.

Although BBHs merging at redshift $z=0$ with at least one component formed from a progenitor which underwent PPISNe are intrinsically rare, GW observations strongly favor massive systems. Our models predict that $20$–$50$% of observed BBHs will have at least one component formed from a progenitor that underwent PPISNe (see e.g. Figure 8). This lends support to the hypothesis that some of the observed BBH systems may have formed this way.

We have also calculated the volumetric rates of CCSNe, PISNe and PPISNe predicted by our models as a function of redshift (see Figure 3). The rate of both PISNe and PPISNe track the low metallicity star formation rate. The volumetric rate of CCSNe (including the stars which undergo PPISNe) at redshift $z=0$ in our model is $\sim 2\times 10^{4}$ Gpc*-3* yr*-1*, whilst the rate of PISNe is $\lesssim 10^{-2}$ Gpc*-3* yr*-1*. We find the rate of stars undergoing PPISNe to be $\sim 0.1$ Gpc*-3* yr*-1*. Our rates for CCSNe, PPISNe and PISNe should be treated as lower limits from our model populations, since there are evolutionary channels we have not included, such as supernovae from stellar merger products (Vigna-Gómez et al., 2019) and supernovae in unbound systems. Our CCSNe rate includes those events where a BH forms through complete fallback; such an event may lead to a ‘failed’ (and therefore not observed) supernova (Gerke et al., 2015; Adams et al., 2017).

In addition, our rate of systems undergoing PPISNe is likely lower than the actual rate of observed PPISNe, since each star may undergo multiple widely-spaced PPISN events before collapsing (Woosley et al., 2007; Woosley, 2017; Marchant et al., 2018).

PISNe, PPISNe and BBH mergers are inherently rare events in our populations. This leads to large statistical uncertainties on some of our predictions (see Table 4). Broekgaarden et al. (2019) have recently shown how binary population synthesis can much more effectively explore rare populations of astrophysical phenomena. These techniques can be applied to study BBH mergers, PISNe and PPISNe.

5.2. Impact of metallicity specific star formation rate

The BBH merger rate tracks the low metallicity star formation rate, increasing from redshift $z=0$ to $z=2$. The fraction of star formation occurring at low metallicities is poorly constrained observationally (e.g. Madau & Dickinson, 2014), particularly at the high redshifts ($z>2$) where BBHs may form in large numbers.

Neijssel et al. (2019) show that uncertainties in the metallicity specific star formation rate can change predictions for the BBH merger rate by an order of magnitude (see also Chruslinska et al., 2019), comparable to uncertainties in binary evolution (Dominik et al., 2012). Since the occurrence of PISNe and PPISNe depends strongly on metallicity (see Figure 2), these uncertainties also have a large impact on our predictions of these events. For example, in combination with our default binary evolution assumptions, using the star formation rate from Madau & Dickinson (2014) and the metallicity distribution from Langer & Norman (2006) leads to $\Lambda_{\mathrm{merge}}$ being a factor of $10$ lower than in our default model. This primarily arises due to an increase in the merger rate of low mass ($\mathcal{M}<10M_{\odot}$) BBH mergers. This metallicity-specific star formation history model over-predicts the total BBH merger rate (Neijssel et al., 2019), compared to GW observations The LIGO Scientific Collaboration et al. (2018a).

5.3. Binary evolution

There are also significant uncertainties in massive binary evolution. Previous studies have shown that uncertainties in black hole kicks, the initial mass function, stellar wind mass loss, common envelope evolution and mass transfer, may lead to up to an order of magnitude uncertainty in BBH merger rates (e.g. Dominik et al., 2012). In the models we have presented here, we have only investigated uncertainties in our modelling of PISNe and PPISNe. We have not attempted to quantify the uncertainies in the model predictions due to uncertain physics. We will investigate what constraints current GW observations place on binary evolution in the future.

5.4. Stellar models

Another source of uncertainty inherent in our population synthesis method is the reliance on pre-computed stellar evolutionary models. The stellar models implemented in COMPAS (Pols et al., 1998; Hurley et al., 2000) were originally computed for the mass range of $5$–$50$ M⊙. This means that we needed to extrapolate up to $150$ M⊙ for the present study. Such an extrapolation, while standard practice in population synthesis (e.g. Dominik et al., 2012; Neijssel et al., 2019), may need to be revisited and possibly refined in future work. Additionally, there have been significant updates in massive stellar evolution since the models we rely on were published. For example, the role of rotation has been extensively investigated and shown to be important in driving the evolution of a star, especially at low metallicities (cf. Brott et al., 2011; Szécsi et al., 2015). Another recently studied effect that may influence our results is the possibility of ‘envelope inflation’ in stars with $\gtrsim$ 70 M*⊙* (Sanyal et al., 2015, 2017). This may lead to differences in when and how a binary interacts and thus to differences in the final outcome of our population synthesis code. The update of our code in order to be able to interpolate on-the-fly from pre-computed stellar tracks (cf. Kruckow et al., 2018; Spera et al., 2019) in order to reflect these recent developments in massive stellar evolution is currently under way.

5.5. Literature comparison

We now turn to comparing our predicted event rates to those in the literature. Eldridge et al. (2019) assume the star formation rate from Madau & Dickinson (2014) and the metallicity distribution form Langer & Norman (2006) and find the volumetric rate of CCSNe at redshift $z=0$ in their models to be $\sim 10^{5}$ Gpc*-3* yr*-1*, whilst they predict the PISN rate to be $\sim 10$ Gpc*-3* yr*-1*, and their BBH merger rate is $\sim 10^{2}$ Gpc*-3* yr*-1* at the same redshift. They do not provide a rate for PPISNe. Eldridge et al. (2018) also provide chirp mass distributions using a simplified method of accounting for GW selection effects.

Nicholl et al. (2013) argue that based on a lack of unambiguous observations of local PISNe, the local PISN rate must be less than $6\times 10^{-6}$ of the local CCSN rate. This is compatible with our predicted rates at redshift $z=0$ which give $\mathcal{R}_{\mathrm{PISN}}/\mathcal{R}_{\mathrm{CCSN}}\lesssim 10^{-6}$. Eldridge et al. (2019) find $\mathcal{R}_{\mathrm{PISN}}/\mathcal{R}_{\mathrm{CCSN}}\sim 10^{-4}$ in their model.

PISNe and PPISNe may account for a small fraction of the observed SLSNe (Kozyreva & Blinnikov, 2015; Aguilera-Dena et al., 2018). SLSNe are observable to cosmological distances of redshifts $z\sim 3$–$4$ (e.g. Cooke et al., 2012; Moriya et al., 2018). The rate of SLSNe inferred from observations at redshift $z\sim 0$ is $\mathcal{R}_{\mathrm{SLSNe}}\sim 100$ Gpc*-3* yr*-1*, rising to $\mathcal{R}_{\mathrm{SLSNe}}\sim 1000$ Gpc*-3* yr*-1at $z=2$ before falling again to $\mathcal{R}_{\mathrm{SLSNe}}\sim 500$ Gpc-3* yr*-1*at $z=3$–$4$ (Quimby et al., 2013; Cooke et al., 2012; Prajs et al., 2017; Moriya et al., 2018). Since it is currently unclear what fraction of SLSNe are due to PISNe and PPISNe, a direct comparison is not possible.

Marchant et al. (2018) suggest that in close binaries which go on to form BBHs, the expansion of the star during a PPISN may lead to an additional period of mass transfer, or even common envelope evolution. We have neglected this effect in this work.

6. Conclusion

In this paper we have investigated the effect of PPISNe and PISNe on the BBH mass distribution, which is being unveiled through GW observations (The LIGO Scientific Collaboration et al., 2018b, a). These observations are providing evidence that the maximum mass of BHs in merging BBHs is $\sim 45$ M⊙, and are hinting at an excess of BHs in the mass range $30$–$45$ M⊙ (The LIGO Scientific Collaboration et al., 2018b, a). In this paper, we have investigated whether it is plausible to attribute these features to PISNe and PPISNe in massive binaries.

We use simple analytic fits to more detailed models to describe the relation between pre-supernova helium core mass and final remnant mass (Belczynski et al., 2016b; Woosley, 2017; Marchant et al., 2018; Woosley, 2019) coupled with the population synthesis code COMPAS (Stevenson et al., 2017; Vigna-Gómez et al., 2018; Barrett et al., 2018; Broekgaarden et al., 2019; Neijssel et al., 2019; Bavera et al., 2019).

Our models predict a maximum BH mass of $40$–$50$ M⊙, in agreement with GW observations. More massive pre-supernova helium cores undergo PISNe and are completely unbound, leaving no remnant. BHs in the mass range $35$–$45$ M⊙ form exclusively from progenitors which underwent PPISNe in our models (except in the model of Belczynski et al. (2016b)), and can be readily identified from GW observations. However, GW measurement uncertainties for such massive BBH mergers are typically $>10$ M⊙ (see e.g. The LIGO Scientific Collaboration et al., 2018b), hindering definite association in many cases. We suggest that the more massive black hole in the BBH merger GW150914 (Abbott et al., 2016b) may have formed from a progenitor which underwent PPISNe, along with potentially both components of the BBH merger GW170729 if the ratio of the mass of the two black holes is close to 1 (The LIGO Scientific Collaboration et al., 2018b).

We find that PPISNe lead to a more gradual transition into the PISN mass gap, which will complicate efforts to determine the exact mass coordinate of the gap. Sharp cutoffs can be measured with an accuracy that scales as $1/N$, not $1/\sqrt{N}$ (e.g. Mandel et al., 2014).

As ground-based GW detectors continue to increase in sensitivity, the number of observations of BBHs will continue to grow (Abbott et al., 2018). Their increased reach will enable measurements of the BBH merger rate as a function of redshift (Fishbach et al., 2018), which can then be compared to the models presented in this paper. A subpopulation of BBHs formed from progenitors which underwent PPISNe may be uncovered. This will allow us to study these types of supernovae indirectly using GWs, increasing our understanding of massive stellar evolution.

In this study, we have used the 10 BBH observations from LIGO and Virgo, published in the GWTC-1 catalog (The LIGO Scientific Collaboration et al., 2018b). Recently, Zackay et al. (2019) and Venumadhav et al. (2019) have almost doubled the number of BBH candidates observed during O1 and O2 using different data analysis techniques. Future comparisons to GW observations should include all observations in order to avoid biased astrophysical conclusions.

Additionally, in this paper we have focused on the mass distribution of BBHs; we present our predictions for the spin distribution of binary black holes in Bavera et al. (2019).

Acknowledgements

We thank Bernhard Müller, Eric Thrane and Alexandra Kozyreva for useful discussions and comments on this paper. We also thank the anonymous referee for insightful comments and suggestions. The authors are supported by the Australian Research Council Centre of Excellence for Gravitational Wave Discovery (OzGrav), through project number CE170100004. This work made use of the OzSTAR high performance computer at Swinburne University of Technology. OzSTAR is funded by Swinburne University of Technology and the National Collaborative Research Infrastructure Strategy (NCRIS). D.Sz. accepts funding from the Alexander von Humboldt Foundation. AVG acknowledges funding support from CONACYT. This research made use of Astropy,888http://www.astropy.org a community-developed core Python package for Astronomy (Astropy Collaboration et al., 2013; Price-Whelan et al., 2018).