Merger rate of stellar black hole binaries above the pair instability mass gap

TL;DR

This paper models the merger rate of stellar black hole binaries with masses above the pair instability gap, predicting detection rates for current and future gravitational wave observatories and exploring the potential for multi-band observations.

Contribution

It introduces a fiducial model for black holes beyond the pair instability mass gap, incorporating cosmic star formation and metallicity evolution, and estimates merger rates and detection prospects.

Findings

Detection rate for LIGO/Virgo: 0.4 to 7 per year.

Third-generation detectors could observe 10 to 460 per year.

LISA can detect these binaries up to thousands of years from coalescence.

Abstract

In current stellar evolutionary models, the occurrence of pair instability supernovae implies the lack of stellar black holes (BHs) with masses between about $[60, \, 120] \, \rm{M}_\odot$, resulting in the presence of an upper mass gap in the BH mass distribution. In this Letter, we develop a fiducial model describing BHs beyond the pair instability gap, by convolving the initial mass function and star formation rate with the metallicity evolution across cosmic time. Under the ansatz that the underlying physics of binary formation does not change beyond the gap, we then construct the cosmic population of merging BH binaries. The detection rate of BH binaries with both mass components above the gap, is found to range between $\simeq [0.4,\,7] \, \rm{yr}^{-1}$ for LIGO/Virgo at design sensitivity and $[10, \, 460] \, \rm{yr}^{-1}$ for third-generation ground based detectors, considering the most pessimistic and optimistic scenarios. LISA can detect individually these binaries up to thousands of years from coalescence. The number of events merging in less than four years, which enable multi-band observation in sequence, is expected in the range $[1, \, 20]$. While ET will detect all these events, LIGO/Virgo is expected to detect $\lesssim 50\%$ of them. Finally, we estimate that the gravitational-wave background from unresolved sources in the LISA band may be in principle detected with a signal-to-noise ratio between $ \simeq 2.5$ and $\simeq 80$.

Tables

Table 1. Table 1 Upper panel: Number of mergers detected per year by O1/O2, LIGO/Virgo at design sensitivity and ET for our four different models as labelled in the text. Lower panel: Number of event over four years for LISA, number of events merging within four years to unable multiband observation with LIGO/Virgo at design sensitivity and with ET in case of joint observations. Last row gives the SNR from the stochastic GW background (summing ‘across-gap’ and ‘above-gap’ sub-populations) for the four models. For each model, left (right) column refers to ‘across-gap’ (‘above-gap’) binaries.

	Models
Detection rates	mSFR-mZ		sSFR-sZ		sSFR-mZ		mSFR-sZ
	‘across’	‘above’	‘across’	‘above’	‘across’	‘above’	‘across’	‘above’
Rate O1/O2 $(S/N>8$) $\left[{\mathrm{yr}}^{-1}\right]$	$0.001$	$0.001$	$0.01$	$0.03$	$0.005$	$0.007$	$0.003$	$0.005$
Rate LIGO/Virgo design $(S/N>12$) $\left[{\mathrm{yr}}^{-1}\right]$	$0.1$	$0.4$	$0.9$	$6.9$	$0.4$	$1.9$	$0.3$	$1.6$
Rate ET $(S/N>12$) $\left[{\mathrm{yr}}^{-1}\right]$	$8.1$	$10.7$	$212.8$	$458.5$	$61.7$	$116.3$	$39.8$	$68.2$
Detected events in 4 yr
LISA events ($S/N>8$)	$0.2$	$1.3$	$2.9$	$42.5$	$0.6$	$12.9$	$0.3$	$8.6$
LISA events ($S/N>8$ & )		$0.8$	$0.5$	$19.3$		$5.8$		$4.6$
LIGO/Virgo multiband events ($S/N>12$)		$0.5$	$0.4$	$8.8$		$2.9$		$1.9$
ET multiband events ($S/N>12$)		$0.8$	$0.5$	$19.3$		$5.8$		$4.6$
S/Ngwb background in LISA	2.7		81.5		18.9		14.8
Detected events in 10 yr
LISA events ($S/N>8$)	$0.6$	$6.2$	$6.2$	$152.1$	$1.9$	$33.0$	$1.7$	$27.0$
LISA events ($S/N>8$ & )	$0.1$	$3.6$	$1.5$	$102.6$	$0.3$	$22.1$	$0.6$	$18.2$
LIGO/Virgo multiband events ($S/N>12$)	$1$	$1.6$	$1.2$	$34.7$	$0.2$	$8.6$	$0.5$	$5.6$
ET multiband events ($S/N>12$)	$0.1$	$3.6$	$1.5$	$102.6$	$0.3$	$22.1$	$0.6$	$18.2$
S/Ngwb background in LISA	3.7		117.8		26.9		21.4

Equations

$$\frac{{\mathrm{d}}^3{n}_{\star }}{\mathrm{d}t\mathrm{d}\log M_{\star }\mathrm{d}\log Z}=\frac{f_{\star }(\alpha )}{⟨m_{\star }(\alpha )⟩}\mathrm{SFR}(\text{t})\times p(\log M_{\star })\times p(\log Z).$$

$$\begin{array}{rl}& \frac{{\mathrm{d}}^{3}n}{\mathrm{d}{t}_m\mathrm{d}\log M_1\mathrm{d}q}=\\ & \mathcal{C}{\int }_{t<t_m}\int \frac{{\mathrm{d}}^{4}n}{\mathrm{d}t\mathrm{d}\tau \mathrm{d}\log M_1\mathrm{d}q}\delta (t_m-(t+\tau ))\mathrm{d}\tau \mathrm{d}t.\end{array}$$

$$\int_{5{{\,\rm M}_{\odot}}}^{50{{\,\rm M}_{\odot}}}{\rm d}\log M_{1}\int{\rm d}q\frac{{\rm d}^{3}n}{{\rm d}t_{m}\,{\rm d}\log M_{1}\,{\rm d}q}\bigg{|}_{z=0}=50\,{\rm Gpc}^{-3}\,{\rm yr}^{-1},$$

$$\mathcal{R}(z_m)={\int }_{120{\mathrm{M}}_{\odot }}^{\infty }\mathrm{d}\log M_1\int \mathrm{d}q\frac{{\mathrm{d}}^{3}n}{\mathrm{d}{t}_m\mathrm{d}\log M_1\mathrm{d}q}.$$

$${\mathcal{R}}_{\det }=\int \mathcal{R}(z_m)\frac{1}{1+z_m}\frac{\mathrm{d}{V}_c}{\mathrm{d}{z}_m}{p}_{\det }\mathrm{d}{z}_m$$

$$\frac{\mathrm{d}N}{\mathrm{d}\log M_1\mathrm{d}q\mathrm{d}z\mathrm{d}\ln f_{\mathrm{gw}}}=\frac{{\mathrm{d}}^{3}n}{\mathrm{d}{t}_m\mathrm{d}\log M_1\mathrm{d}q}\frac{\mathrm{d}{V}_{\mathrm{c}}}{\mathrm{d}z}\frac{\mathrm{d}{t}_m}{\mathrm{d}\ln f_{\mathrm{gw}}},$$

$$\frac{\mathrm{d}{t}_m}{\mathrm{d}\ln f_{\mathrm{gw}}}=\frac{5}{96{\pi }^{8/3}}{\left(\frac{c^3}{G\mathcal{M}}\right)}^{5/3}(f_{\mathrm{gw}}(1+z){)}^{-8/3}.$$

$$(S/N{)}^2=\int \frac{|\tilde{h}(f_{\mathrm{gw}},M_1,q,z){|}^2}{S_{\mathrm{n}}(f_{\mathrm{gw}})}\mathrm{d}{f}_{\mathrm{gw}},$$

$$(S/{\mathrm{N}}_{\mathrm{gwb}}{)}^2=T\int \gamma (f_{\mathrm{gw}})\frac{h_{c,\mathrm{gwb}}^4}{f_{\mathrm{gw}}^2{S}_n(f_{\mathrm{gw}}{)}^2}\mathrm{d}{f}_{\mathrm{gw}}$$

$${\Omega }_{\mathrm{gw}}(f_{\mathrm{gw}})=\frac{2}{3}{\left(\frac{\pi f_{\mathrm{gw}}{h}_{c,\mathrm{gwb}}}{H_0}\right)}^2,$$

Full text

11affiliationtext: Department of Physics G. Occhialini, University of Milano - Bicocca, Piazza della Scienza 3, 20126 Milano, Italy22affiliationtext: National Institute of Nuclear Physics INFN, Milano - Bicocca, Piazza della Scienza 3, 20126 Milano, Italy

1. Introduction

During the first and second observing runs, the LIGO-Virgo scientific collaboration (Aasi et al., 2015; Acernese et al., 2014) detected the gravitational-wave (GW) signals from the coalescence of ten stellar black hole binaries (BHBs) with individual masses between $7.7^{+2.2}_{-2.6}$ and $50.6^{+16.6}_{-10.2}{{\,\rm M}_{\odot}}$ (Abbott et al., 2018). The observed events can be reproduced by stellar population synthesis codes (see e.g. Podsiadlowski et al., 2003; Postnov & Yungelson, 2014; Dominik et al., 2015; Spera et al., 2015; Belczynski et al., 2016; Schneider et al., 2017; Marassi et al., 2019) in which a key role is played by the metallicity evolution along the cosmic history. In fact, it is widely accepted that low metallicity stars experience negligible mass loss during their lifetimes, due to their weaker stellar wind, thus collapsing in heavy BH remnants, consistent with those discovered by LIGO-Virgo (Abbott et al., 2016a, 2017a).

The occurrence of pulsational pair-instability supernovae (PPISNe) and pair-instability supernovae (PISNe) in massive, low metallicity stars (with $\rm Z\lesssim 0.002$, Heger & Woosley, 2002; Woosley, 2017; Yoshida et al., 2016; Marchant et al., 2018) is expected to enhance the formation of BHs in the mass range $30\lesssim M_{\rm{rem}}/{{\,\rm M}_{\odot}}\lesssim 50$, leading to a pile up around $\sim 45{{\,\rm M}_{\odot}}$ (Stevenson et al., 2019). Current GW data indicate an excess of BHs in the interval $30-45{{\,\rm M}_{\odot}}$ (Abbott et al., 2018), that future observations can confirm or challenge (Fishbach & Holz, 2017; Talbot & Thrane, 2018). Above $\sim 50{{\,\rm M}_{\odot}}$, a cut-off or edge is expected in the BH mass function, as PISNe lead to the explosion of the star preventing the formation of a massive BH remnant. However, stars with zero-age main sequence (ZAMS) mass $M_{\rm{ZAMS}}\gtrsim 260{{\,\rm M}_{\odot}}$ and absolute metallicities $\lesssim 10^{-3}$ avoid disruption, as they develop massive CO cores that directly collapse into a BH of $M_{\rm rem}\gtrsim 100{{\,\rm M}_{\odot}}$ (Woosley et al., 2002; Uchida et al., 2019). Those systems have been invoked as viable seeds of supermassive BHs in the high redshift universe (Volonteri, 2010; Valiante et al., 2016), but have so far been ignored in population synthesis models used to interpret LIGO-Virgo detections, that customarily evolve stars only up to $100-150{{\,\rm M}_{\odot}}$ (Belczynski et al., 2016; Mapelli et al., 2019; Neijssel et al., 2019). Conversely, several alternative mechanisms to produce BHs above the PISNe cut-off have been proposed (Gerosa & Berti, 2017; Rodriguez et al., 2018; McKernan et al., 2019; Yang et al., 2019). Moreover, Spera et al. 2019 and Di Carlo et al. 2019 proposed that BHs in the pair-instability gap may originate from the direct collapse of massive stars with large envelope and small core masses, thus avoiding the pair-instability phase.

If BHs above the PISNe “upper-mass” gap do indeed form, pair and coalesce in binaries, they could be potentially detectable with the Laser Interferometer Space Antenna (LISA, Amaro-Seoane et al. 2017), or from third-generation ground-based detectors such as the Einstein Telescope (ET, Punturo et al., 2010). These binaries may also contribute to the stochastic GW background (GWB) between 0.5 mHz and 20 mHz, hampering observations of individual sources close to the LISA bucket (Caprini et al., 2019).

In this Letter, we estimate the merger rate of stellar BHBs from isolated field binary evolution across the mass spectrum and beyond the upper-mass gap. Assuming the gap as a sharp cut-off at $[60,\,120]{{\,\rm M}_{\odot}}$, we distinguish three sub-populations for the binaries: the ‘above-gap’ (‘below-gap’) binaries with both components above (below) the upper (lower) edge of the mass gap and the ‘across-gap’ binaries with one component above and one below the mass gap. The population of ‘below-gap’ binaries is found to be consistent with previous studies (e.g., Sesana, 2016; Gerosa et al., 2019) and is not considered here. For the ‘across-gap’ and ‘above-gap’ sub-populations, we report detection rates with ground- and space-based detectors and estimate their contribution to the stochastic GWB in the LISA band.

2. Models

Our approach builds on knowledge of the ‘below-gap’ BHB population, extending it to BHB properties above the PISNe gap.

We evolve single stars using the stellar evolution code SEVN (Spera & Mapelli, 2017, and reference therein). The code includes up-to-date stellar winds, SN explosion models, PISNe and PPISNe prescriptions, and provides BH remnant masses as a function of the mass of the progenitor stars and of the absolute metallicity, in the range $2\times 10^{-4}<\rm{Z}<2\times 10^{-2}$ (see Fig. 2 in Spera & Mapelli 2017 for the relation between initial stellar mass and remnant mass for different metallicities).

We consider two main models for the star formation rate (SFR) and the evolution of the mean metallicity $\langle Z\rangle$ of the galaxy population across cosmic history. The first model takes both from Madau & Fragos (2017) and is labelled as “mSFR-mZ”. The second model adopts the SFR as in Strolger et al. (2004) and the metallicity for the intergalactic medium reported in Madau & Dickinson (2014), shifted to match Madau & Fragos (2017) metallicity at $z=0$. This accounts for a possible rapid decline of the metallicity between the present and redshfit $z\simeq 4$. The model is labelled as “sSFR-sZ”. The SFR and mean metallicity $\langle Z\rangle$ versus redshift are shown in Fig. 1 for these two models. We also consider two additional intermediate models, combining the SFRs and the metallicity prescriptions, thus labelled “mSFR-sZ” and “sSFR-mZ” (see also Chruslinska et al. 2018; Neijssel et al. 2019 for further discussions about uncertainties on the SFR and metallicity distribution).

We assume a stellar initial mass function (IMF) $\xi(M_{\star},\alpha)\propto M_{\star}^{-\alpha}$ between $[8,\,350]{{\,\rm M}_{\odot}}$, with $\alpha=2.7$ for the SFR from Madau & Fragos 2017 and $\alpha=2.35$ for the SFR from Strolger et al. 2004. The differential comoving volume number density formation rate of progenitor stars as a function of cosmic time, mass and redshift is

[TABLE]

Here $f_{\star}$ and $m_{\star}$ are the fraction of simulated binaries and the average IMF mass as defined by Equation (6) and (7) of Marassi et al. (2011), SFR is the cosmic star formation rate density at cosmic time $t$, and $p(\log M_{\star})$, $p(\log Z)$ are the probability densities of stellar mass and metallicity. The former is directly proportional to the IMF, while the latter is taken at each redshift to be a log-normal distribution centred along either the “mZ” or “sZ” relations (as shown in Figure 1) with dispersion $0.5$ dex. For a given interval $(\Delta\log M_{\star},\,\Delta\log Z)$, we evolve a single star with SEVN to determine its BH final mass. In this way, Equation (1) is mapped into the relic BH formation rate density, ${\rm d}^{2}n/({\rm d}t\,{\rm d}M_{1})$. The primary BH mass, $M_{1}$, of each BHB is drawn from this distribution.

To convert the formation rate of BHs into the merger rate of BHBs we make two simple assumptions: i) all BHs are in binaries with secondary BH drawn according to a flat mass ratio $q=M_{2}/M_{1}<1$ distribution in the range $[0.1,\,1]$111Binaries with $M_{2}$ falling in the mass gap are discarded and the remaining population is re-normalized to match the total mass density of BHB produced to be equal to $\int{{\rm d}\log M_{1}}\,M_{1}[{\rm d}^{2}n/({\rm d}t\,{\rm d}\log M_{1})]$., and ii) mergers occur at a time $t_{m}=t+\tau$ 222Here $\tau$ includes the evolution time of the primary star. where the delay time $\tau$ is distributed according to $p(\tau)\propto\tau^{-1}$ ) (Dominik et al., 2012; Neijssel et al., 2019) between $\tau_{\rm\,min}=50{{\,\rm Myr}}$ and $\tau_{\rm\,max}=t_{\rm Hubble}$ (Dvorkin et al., 2016), where $t_{\rm Hubble}$ is the Hubble time. We also explore the possibility for a flat mass ratio distribution in the range [0.5,1] (Fishbach & Holz, 2019). The rate density per comoving volume of merging BHBs is therefore given by:

[TABLE]

The normalization constant ${\cal C}$ is set to ensure that the intrinsic BHB merger rate in the local Universe is

[TABLE]

close to the best estimate provided by the LIGO-Virgo O2 run (The LIGO Scientific Collaboration et al., 2018). This a posteriori normalization is needed because of the very simplistic assumptions made above. We checked, however, that both the resulting BHB merger rate density as a function of redshift and the mass function of merging BHBs below the pair-instability gap are in good agreement with sophisticated population synthesis models found in the literature (e.g. Spera et al., 2019) .

We are here interested in binaries with at least one BH above the pair-instability gap. We define their merger rate as

[TABLE]

Depending on $q$, the secondary can be either below or above the mass gap, thus defining the sub-classes of ‘across-gap’ and ‘above-gap’ BHBs introduced above. The number of detections per year is then computed as

[TABLE]

where $(1+z_{m})^{-1}={\rm d}t_{m}/{\rm d}t_{m}^{\rm obs}$ accounts for the time dilation between the source and the observer frames and ${\rm d}V_{c}/{\rm d}z_{m}$ is the differential comoving volume shell. Finally, $p_{\rm det}$ represents the detection probability of a random-oriented binary with a given $M_{1}$, $q$ and $z_{m}$ for a threshold signal-to-noise ratio (S/N) (Abadie et al., 2010).

Although Equation (5) is appropriate for the detection rate of ground-based interferometers, LISA will also see persistent sources, caught several years before coalescence. The distribution of observed sources in the LISA band is simply given by

[TABLE]

where ${\rm d}t_{m}/{\rm d}\ln f_{\rm gw}$ is given by the quadrupole approximation for circular orbit (Peters, 1964) as

[TABLE]

Here $\mathcal{M}=(M_{1}M_{2})^{3/5}/(M_{1}+M_{2})^{1/5}$ is the source frame chirp mass and $f_{\rm gw}$ the observed GW frequency 333We assume circular BHBs because several processes acting during stellar evolution (e.g. tidal circularisation, common envelope evolution, etc.) and long delay times are expected to produce nearly circular BHB orbits.. For each BHB population model described at the beginning of this section, equation (6) is used to draw 10 Monte Carlo realisations of the BHB distribution across the observed frequency spectrum in the range $[10^{-4},\,10^{-1}]\,\rm Hz$. Each sample is then taken to represent the distribution of sources in the sky at the start of the LISA mission.

All sources are then evolved forward in time assuming GW-driven dynamics and their S/N in the detector is evaluated as

[TABLE]

where $\tilde{h}(f_{\rm gw})$ is the Fourier transform of the GW strain and $S_{\rm n}(f_{\rm gw})$ is the power spectral density of the detector 444For O2 LIGO sensitivity, we adopt the curve labelled as ’mid’ in Abbott et al. 2016b. We also consider Advanced LIGO (aLIGO, Harry & LIGO Scientific Collaboration, 2010), Advanced Virgo (AdVirgo, Acernese et al., 2014) and Einstein Telescope (ET-D, Hild et al., 2011). For LISA we adopt the curve described in Robson et al. 2019.. Note that the integral in Eq. (8) is over the frequency interval covered by the source over the observation time. For each value of $M_{1},q$ and $z_{\rm m}$, we compute the S/N randomizing over sky-position, polarization and inclination angles, and assuming non-spinning BHs. For ground-based detectors, we compute the S/N with the LALsuite (LIGO Scientific Collaboration, 2018). We model the inspiral-merger-ringdown signal with the IMRPhenomD waveform (Husa et al., 2016; Khan et al., 2016). For O1/O2, we consider an event to be detectable if $\rm S/N>8$ (Abbott et al., 2018), while for LIGO/Virgo at design sensitivity and ET we assume $\rm S/N>12$. Similarly, for LISA we compute the S/N with the IMRPhenomC waveform (Santamaría et al., 2010) with $\rm S/N>8$. To estimate the rates for multiband events, we consider only the events detected in LISA and coalescing in $T_{\rm gw}<4\,\text{yr}$, where $T_{\rm gw}$ is the merger timescale due to GW emission. We also consider a possible extended time mission of 10 years.

We compute the level of the stochastic GWB generated by the inspiralling BHBs, at each frequency by summing in quadrature the characteristic strains of all unresolved sources, i.e. binaries with $\rm S/N<8$. Then, the signal power S/Ngwb is evaluated following Thrane & Romano (2013) and Sesana (2016)

[TABLE]

where $T=4$ yr is the LISA mission required lifetime, $h^{2}_{c,\rm gwb}(f_{\rm gw})=2f_{\rm gw}S_{h}(f_{\rm gw})$ (being $S_{h}(f_{\rm gw})$ the power spectral density of the signal), and $\gamma(f_{\rm gw})=1$ (see Fig. 4 in Thrane & Romano, 2013). We estimate the strength of the GWB through its GW energy density parameter

[TABLE]

where $H_{0}$ the Hubble’s constant.

3. Populations, Rates and GW Background

In Fig. 2 we show the number density of BHBs formed per unit comoving volume versus $\mathcal{M}$, for different metallicities and sub-populations, and the corresponding merger rate for models mSFR-mZ and sSFR-sZ. Given the distribution of the mass ratio adopted, there is no evident gap in the source-frame chirp mass, and the ‘across-gap’ and ‘above-gap’ sub-populations are the result of poor-metal stars with $\langle Z\rangle<1.2\times 10^{-3}$. For the ‘below-gap’ sub-populations, the outcome of this analysis is fairly consistent with that of Spera et al. (2019) (see their Fig. D2). At $\langle Z\rangle>8\times 10^{-3}$ our maximum chirp mass is close to theirs, while at $\langle Z\rangle<8\times 10^{-3}$, we obtain larger chirp masses, in the range $[40,\,50]{{\,\rm M}_{\odot}}$. This is expected due to the difference from single to binary evolution. However we note that $\mathcal{M}\lesssim 50{{\,\rm M}_{\odot}}$ have been recovered in alternative population synthesis codes (Chruslinska et al., 2018).

The comparison between the two models shown in Fig. 2 highlights the impact of metallicity on the number density of heavy BHBs. Model sSFR-sZ predicts a rapid decline in the metallicity versus redshift and, as a consequence, the ‘across-gap’ and ‘above-gap’ sub-population rates are one order of magnitude higher than in model mSFR-mZ. This is also evident in the right panels showing the merger rate density of the three sub-populations for each model. Note that in both models, the ‘above-gap’ sub-population produces slightly more mergers than the ‘across-gap’ one, but the total merger rate is always heavily dominated by the ‘below-gap’ systems. Mixed models ( mSFR-sZ and sSFR-mZ, not shown) give intermediate results, as expected.

In Fig. 3 we show the overall properties of the events observed by LISA with $\rm S/N>8$. The majority of the events concentrates at $z\lesssim 0.5$ with the tail extending up to $z\approx 1.5$ in model sSFR-sZ. In addition, for both models, systems merging within the LISA lifetime (dashed lines) are detected at slightly higher redshifts. The source-frame chirp mass distributions cluster around two peaks that broadly correspond to the ‘across-gap’ (left peak) and ‘above-gap’ (right peak) sub-populations. It is also evident that the event number in the ‘above-gap’ group is nearly $\times 10$ higher than the ‘across-gap’ one. This is simply because more massive and nearly equal-mass binaries produce louder GW signal and can be seen further. Note that, in these models, we do not expect any detected ‘across-gap’ binaries merging in 4 years. Most of the detected sources have an initial $f_{\rm gw}$ around $\sim 3\times 10^{-3}$ Hz, while only for model sSFR-sZ we can observe around $\sim\mathcal{O}(1)$ sources down to $10^{-4}$ Hz. Obviously, merging BHBs (dashed lines) peak at a slightly higher frequencies with respect to persistent ones (solid lines) for both models.

Fig. 4 shows the energy density of the stochastic GWB – $\Omega_{\rm gw}$, of Equation (10)– as a function of observed frequency for all the considered models. Except for model mSFR-mZ, the signal is strong enough to be separated from the LISA detector noise. In the most optimistic model, $\Omega_{\rm gw}$ is comparable to that inferred from the ‘below-gap’ BHB and neutron star binary populations (e.g. Farmer & Phinney, 2003; Sesana, 2016), thus significantly contributing to the family of unresolved LISA astrophysical foregrounds. Our result is also consistent with the upper limits reported for O1/O2 from the LIGO/Virgo collaboration (Abbott et al., 2017b; The LIGO Scientific Collaboration et al., 2019).

Besides posing another compelling scientific case for LISA, the population of ‘above-gap’ BHBs is also important for ground-based detectors. This is quantified in Tab. 1, summarizing all the relevant figures of this study. The first three rows of the table report the merger rates, computed according to Eq. (5) for LIGO/Virgo at O2, design sensitivity and for ET. Current non detection of ‘across-gap’ and ‘above-gap’ binaries during 01/02 is consistent with our models. LIGO/Virgo at design sensitivity can detect between $\simeq 0.5$ and $\simeq 7$ ‘above-gap’ events per year for the pessimistic and optimistic models, respectively. Intermediate models predict $\approx 2$ events per year. For ET, the rate increases by more than an order of magnitude for all sub-populations, independently of the model. The number of detected events range from $\simeq[10,\,460]\,{{\,\rm yr}^{-1}}$. Due to the improved sensitivity, ET would also be able to detect several ‘across-gap’ binaries per year in all models.

For a 4 year LISA mission, the number of detected ‘above-gap’ events ranges between $\simeq 1.5$ and $43$ for the two limiting models. The number of events detected with $\rm S/N>8$ and merging in the mission lifetime is reduced but it’s still of order unity even in the pessimistic model. In particular, the case sSFR-sZ predicts $\simeq 20$ events, while the two intermediate models predict around $\simeq 5$ events. These sources are the best candidate for multiband detection. If LISA will join LIGO/Virgo, we expect to detect from $\simeq 0.5$ to $\simeq 10$ multiband binaries. If ET will be operative instead, all the sources detected in LISA will also be detected at later times by ET. If LISA reaches the 10 years mission goal, the number of detected events increase by a factor of $\gtrsim 3$. This can be explained considering that the increase in SNR from a longer time mission goes as $\sim\sqrt{10/4}$, which translates in almost a factor of four in the accessible volume and, therefore, in the number of sources, assuming a constant merger rate.

We also test the case for a flat mass ratio distribution for $q\in[0.5,\,1]$, which seems to be favoured by current LIGO-Virgo detections (Abbott et al., 2018). Under this assumption, we lose the ‘across-gap’ sub-populations, since ‘across-gap’ systems would necessarily have $q<0.5$, but we find a significant increase of $\approx 2$ in the rate for the ‘above-gap’ sub-population.

Finally, we run a Fisher Matrix code to estimate the uncertainties in the parameter estimation of these populations in the LISA band. The parameter estimation is performed with an 8 post-Newtonian (8PN) frequency domain waveform for circular precessing BHs (Klein et al., 2014). Due to the larger number of cycles in band and to LISA orbital motion, we expect to be able to localize multiband sources with a median precision of $\lesssim 10\,\rm deg^{2}$ and determine single component redshifted masses to better than $1\%$ percent. The luminosity distance is determined with a median uncertainty of $\lesssim 30\%$, while the spins would be essentially unresolved.

4. Conclusion

The existence of stellar BH above the pair-instability gap is an uncharted territory, as there are many uncertainties on the formation and evolution of very massive stars, on the gravitational collapse and nuclear energy production, on binary formation in galactic fields and dynamical systems. Any detection of BHBs belonging to this population will be of capital importance in the understanding of the extreme physics governing the evolution of massive stars. In this Letter we performed a first attempt to quantify the population of ‘above-gap’ BHBs observable with current and future ground- and space-based detectors. Since it works under the ansatz that ‘above-gap’ BHBs form and evolve abiding the same physics governing the evolution of ‘below-gap’ systems, our approach can be considered agnostic. As such, although employed for the field formation scenario, it can be in principle extended to alternative formation channels (e.g. dynamical capture), where massive BHs form from multiple GW-driven mergers (Rodriguez et al., 2019).

We find that prospects for discovering ‘above-gap’ BHBs are interesting, with several systems detectable either from the ground or from space (see numbers in Tab. 1). Moreover, ‘above-gap’ binaries are primary candidates for multi-band detection, with up to ${\cal O}(100)$ observable systems for a detector network featuring a 10 yr LISA mission plus 3G interferometers on the ground (Kalogera et al., 2019). A critical ansatz underlying those numbers is the extension of the Kroupa IMF up to $350{{\,\rm M}_{\odot}}$. The formation of such massive stars is still puzzling, although there is some observational evidence that the stellar IMF might in fact extend to $M\gtrsim 300{{\,\rm M}_{\odot}}$ (Crowther et al., 2010). We note that if massive stars do not form, the merger rate for ‘above-gap’ binaries from this channel has to be zero. On the other hand, future detections of ‘above-gap’ binaries can empirically prove that massive stars form.

Finally, we stress that the numbers presented here are subject to large uncertainties, stemming from several pieces of poorly known underlying physics, including the detailed evolution of massive star, the cosmic evolution of stellar metallicity and its dependence on the galactic environment (Pezzulli et al., 2016; Valiante et al., 2017), the mass and mass-ratio distribution of binaries and so on. We considered in this Letter a minimal set of models, bracketing some of the critical uncertainties. An extended study, exploring the whole parameter space is ongoing and will be the subject of a future publication.

The authors acknowledge P. Madau, R. Schneider, C. Berry, A. Lupi and the anonymous referee for insightful discussions and comments. A. M., M. B. and M. C. acknowledge partial financial support from the INFN TEONGRAV specific initiative.