Constraining the Fraction of Binary Black Holes Formed in Isolation and Young Star Clusters with Gravitational-wave Data

The initial magnitude and direction of BH spins is still a matter of debate; e.g., see Miller & Miller (2015) for a review. Overall, the dependence of the spin magnitude of a BH on the spin magnitude of the progenitor star (or stellar core) is largely unknown.

As for the direction of the spins, most binary evolution processes in isolated binaries tend to favor the alignment of stellar spins with the orbital angular momentum of the binary. In isolated binaries, supernova kicks are the leading mechanism to substantially tilt the spin axes with respect to the orbital plane (Kalogera 2000). In star clusters, dynamical exchanges tend to reset the memory of the initial binary spin. Thus, we expect the spins of exchanged BBHs to be isotropically distributed. The spin direction of original BBHs in star clusters is expected to fall somewhere in between because, on the one hand, these BBHs participate in the dynamical evolution of the star cluster (ergo dynamical encounters can affect the initial spin orientation), while on the other hand, they form from the evolution of stellar binaries (ergo binary evolution processes tend to realign the spins).

Given these considerable uncertainties on both spin magnitude and direction, we decided not to embed detailed spin models in our population synthesis and dynamical simulations. Spins are added to our simulations in post-processing, assuming simple toy models. Dimensionless spin magnitudes a (defined as $a=| J| \,c/G\,{m}_{\mathrm{BH}}^{2}$, where J is the BH spin, c is the speed of light, G is the gravitational constant, and m_BH is the mass of the BH) are randomly drawn from a Maxwellian distribution with root mean square equal to 0.1. With this choice, the median spin is a ∼ 0.15 and the distribution quickly fades off for a > 0.4. Hereafter, we refer to this model as "low-spin" (LS). We assume this distribution for spin magnitudes because the results of the first two aLIGO/aVirgo runs disfavor distributions with large spin components aligned (or nearly aligned) with the orbital angular momentum (Abbott et al. 2019a).

For comparison, we also consider a second, rather extreme case, in which spin magnitudes are drawn from a Maxwellian distribution with root mean square equal to 0.3 (we reject spin magnitudes $a\gt 0.998$). With this choice, the median spin is a ∼ 0.46. Hereafter, we refer to this model as "high-spin" (HS).

Regarding spin orientations, we assume that BH spins in both isolated BBHs and original BBHs are perfectly coaligned with the orbital angular momentum of the binary (e.g., Rodriguez et al. 2016c): there are large uncertainties on the kicks imparted on newly formed BHs, but recent work (Gerosa et al. 2018) shows that the fraction of BBHs with negative effective spins is at most ∼20%. Moreover, we neglect the effect of dynamical perturbations on original BBHs because their main properties are similar to isolated BBHs: they have nearly the same chirp mass, mass ratio, and eccentricity distribution (Di Carlo et al. 2019).

Finally, BH spins in exchanged BBHs are randomly drawn isotropically over a sphere. Our assumptions for the spin models are summarized in Table 1.

The redshift parameter was not computed self-consistently in the set of astrophysical simulations generated for this study—because we only consider a fixed value for the metallicity, and the stellar metallicity is a crucial ingredient of redshift evolution (Mapelli et al. 2017). A self-consistent redshift evolution will be included in future work (V. Baibhav et al. 2019, in preparation).

Here, we opted to exclude redshift information from our statistical analysis (see Section 4). However, we still need to prescribe a redshift probability distribution function to estimate selection effects. As this is a simple toy model, we assume that the redshift is distributed uniformly in comoving volume and source-frame time, i.e.,

$$\begin{eqnarray}&&p(z)\propto \displaystyle \frac{1}{1+z}\displaystyle \frac{{\text{}}{{dV}}_{c}}{{\text{}}{dz}}.\end{eqnarray} \tag{ 5 }$$

We consider redshifts in the range z ∈ [0, 2] for both second- and third-generation detectors, postponing more accurate modeling to future work.

We estimate selection effects using the semianalytic approach of Finn & Chernoff (1993), but see also Dominik et al. (2015), Chen et al. (2017), and Taylor & Gerosa (2018). We associate a detection probability ${p}_{\det }(\lambda )\in [0,1]$ with any given GW source with parameters λ. A source is detectable if its signal-to-noise ratio (S/N)

$$\begin{eqnarray}&&\rho =4{\int }_{0}^{+\infty }\displaystyle \frac{| \tilde{h}(f){| }^{2}}{{S}_{n}(f)}{\text{}}{df}\end{eqnarray} \tag{ 6 }$$

exceeds a given threshold ρ_th, with $\tilde{h}(f)$ being the gravitational waveform in the frequency domain and S_n(f) the one-sided noise power spectral density of the detector. To compute $\tilde{h}(f)$, we have used the IMRPhenomD model (Khan et al. 2016), which is a phenomenological waveform model describing the inspiral, merger, and ringdown of a nonprecessing BBH merger signal. We consider the noise power spectral density curves corresponding to both the aLIGO (Abbott et al. 2018) and the Einstein Telescope (ET; Abbott et al. 2017a) at their design sensitivity. We implement a single-detector S/N threshold ρ_th = 8, which was shown to be a good approximation of more complex multidetector analysis based on large injection campaigns; see Abadie et al. (2010), Abbott et al. (2016a), and Wysocki et al. (2018), for more details. Both waveforms and detector sensitivities were generated using pycbc (Dal Canton et al. 2014; Usman et al. 2016).

For each binary in our catalogs, we estimate the optimal S/N, ${\rho }_{\mathrm{opt}}$, using Equation (6). This corresponds to a face-on source located overhead with respect to the detector. The S/N of a generic source is given by ρ = ω × ρ_opt, where ω encapsulates all the dependencies on sky-location, inclination, and polarization angle (Finn & Chernoff 1993; Finn 1996). A source located in a blind spot of the detector yields a value of ω = 0, while an optimally oriented source has ω = 1. The probability of detecting a source is then expressed as

$$\begin{eqnarray}&&{p}_{\det }(\lambda )={\mathbb{P}}(\rho \geqslant {\rho }_{\mathrm{thr}})\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&=\,{\mathbb{P}}(\omega \geqslant {\rho }_{\mathrm{thr}}/{\rho }_{\mathrm{opt}})\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&=\,1-{F}_{\omega }({\rho }_{\mathrm{thr}}/{\rho }_{\mathrm{opt}}),\end{eqnarray} \tag{ 9 }$$

where F_ω is the cumulative distribution function of ω. This function was computed via Monte Carlo methods as implemented in the Python package gwdet (Gerosa 2018). The function F_ω is set explicitly to 1 for ρ_opt < ρ_thr, which gives the expected detection probability p_det(λ) = 0 for events that are too quiet to be observed.

The noise contained in the data d of a GW detector results in errors on the measurement of the parameters λ of a GW source. From a Bayesian point of view, these errors are fully described by the posterior distribution $p(\lambda | d)$. We make use of posterior distributions for the first 10 GW events publicly released by Abbott et al. (2019b). We also generate mock observations from our catalogs to forecast future scenarios with a growing number of events. In this case, running a full injection campaign to estimate measurement errors would be too computationally expensive and out of the scope of this study. For simplicity, we approximate posterior distributions with simple Gaussians (Farr et al. 2017; Gerosa & Berti 2017):

$$\begin{eqnarray}&&p({\lambda }_{i})={ \mathcal N }(\overline{{\lambda }_{i}},{\sigma }_{i}).\end{eqnarray} \tag{ 10 }$$

The mean $\overline{{\lambda }_{i}}$ are obtained by first extracting a value λ_i^T from our astrophysical models. We then use prescriptions⁹ by Mandel et al. (2017), but see also Stevenson et al. (2017a) and Powell et al. (2019). The injected chirp mass M_c and symmetric mass ratio η are given by:

$$\begin{eqnarray}&&\overline{{{ \mathcal M }}_{c}}={{ \mathcal M }}_{c}^{T}\left[1+\alpha \displaystyle \frac{8}{\rho }{r}_{0}\right],\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&\overline{\eta }={\eta }^{T}\left[1+0.03\displaystyle \frac{8}{\rho }{r}_{0}\right],\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&\overline{{\chi }_{\mathrm{eff}}}={\chi }_{\mathrm{eff}}^{T}+\beta \displaystyle \frac{8}{\rho }{r}_{0},\end{eqnarray} \tag{ 13 }$$

where ${r}_{0}\sim { \mathcal N }(0,1)$, β = 0.1 and α takes values of 0.01, 0.03, and 0.1 for η^T > 0.1, 0.1 > η^T > 0.05, and η^T < 0.05, respectively, and we convert $(\overline{{{ \mathcal M }}_{c}},\overline{\eta })\to (\overline{{M}_{t}},\overline{q})$. Finally, the values for the standard deviations σⁱ are set using the same prescriptions as Mandel et al. (2017) given in the previous set of equations.

Measurement errors for ET are obtained by rescaling the aLIGO results using the S/N

$$\begin{eqnarray}&&{\sigma }_{\mathrm{ET}}={\sigma }_{\mathrm{aLIGO}}\displaystyle \frac{{\rho }_{\mathrm{aLIGO}}}{{\rho }_{\mathrm{ET}}},\end{eqnarray} \tag{ 14 }$$

as expected in the large-S/N limit (Poisson & Will 1995).

The general expression for the rate of a given model with parameters θ can be written as

$$\begin{eqnarray}&&\displaystyle \frac{{\text{}}{dN}}{{\text{}}d\lambda }(\theta )=N(\theta )p(\lambda | \theta ),\end{eqnarray} \tag{ 15 }$$

where N(θ) is the total number of sources predicted by the model and $p(\lambda | \theta )$ is the normalized model distribution or rate.

From the catalog of sources presented in Section 2.4, we approximate the normalized rates p using kernel density estimation (KDE) methods. Gaussian kernels with a bandwidth parameter of 0.05 on the data set $\{{M}_{t},q,{\chi }_{\mathrm{eff}}\}$ are capable of accurately reproducing the distributions in Figure 1 for both formation channels and spin models.

Statistical inference is implemented with a standard Bayesian hierarchical model. Our analysis is based on the formalism already presented by Loredo (2004), Mandel et al. (2019), and Taylor & Gerosa (2018). In a nutshell, the posterior distribution on the model parameters θ marginalized over N(θ) is

$$\begin{eqnarray}&&p(\theta | d)\propto p(\theta )\displaystyle \prod _{k=1}^{{N}_{\det }}\displaystyle \frac{\int p(\lambda | \theta )p(\lambda | d)/p(\lambda )\,{\text{}}d\lambda }{\int {p}_{\det }(\lambda )\,p(\lambda | \theta )\,{\text{}}d\lambda },\end{eqnarray} \tag{ 16 }$$

where N_det is the number of entries in the detection catalog, $p(\lambda | \theta )$ describes the astrophysical model, p(θ) is the prior on each astrophysical model, $p(\lambda | d)$ is the posterior of an individual GW event, p(λ) is the prior used in the single-event analysis, and p_det(λ) describes selection effects.

If the posterior $p(\lambda | d)$ is provided in terms of Monte-Carlo samples λ_i, as in Abbott et al. (2019b), we can rewrite Equation (16) as

$$\begin{eqnarray}&&p(\theta | d)\propto p(\theta )\displaystyle \prod _{k=1}^{{N}_{\det }}\displaystyle \frac{{\displaystyle \sum }_{i}p({\lambda }_{i}| \theta )/p({\lambda }_{i})}{\int {p}_{\det }(\lambda )\,p(\lambda | \theta )\,{\text{}}d\lambda }.\end{eqnarray} \tag{ 17 }$$

For the case of our mock Gaussian posteriors, Equation (16) becomes

$$\begin{eqnarray}&&p(\theta | d)\propto p(\theta )\displaystyle \prod _{k=1}^{{N}_{\det }}\displaystyle \frac{\int p(\lambda | \theta ){ \mathcal N }({\lambda }_{k},{\sigma }_{k})/p(\lambda )\,{\text{}}d\lambda }{\int {p}_{\det }(\lambda )\,p(\lambda | \theta )\,{\text{}}d\lambda }.\end{eqnarray} \tag{ 18 }$$

The denominator $\beta (\theta )\equiv \int {p}_{{\det }}(\lambda )\,p(\lambda | \theta )\,{\rm{d}}\lambda$ depends only on the model θ, and not on the event parameters λ. In practice, we estimate β(θ) by generating values for (M_t, q) using rejection sampling from our KDE approximation of the astrophysical rates, and extracting values for the aligned components of the spins and the redshift as described in Sections 2.3.1 and 2.3.2. As the dynamical catalog is formed of both exchanged and originals BBHs, we consider the two subpopulations separately and combine them in the proportion predicted by our dynamical simulations.

We want to assess the number of observations N_obs needed to discriminate between the two models, assuming that one of them is a faithful representation of reality. To understand how each parameter impacts the analysis, we run our statistical pipeline assuming we only measure either $\{{M}_{t}\}$, $\{{M}_{t},q\}$, or $\{{M}_{t},q,{\chi }_{\mathrm{eff}}\}$. In practice, this implies that the integral in Equation (18) is evaluated on the selected variables while marginalizing over the others. Regarding the production of mock data: we generated 10³ sets of observations for each value of N_obs, in order to produce a statistical estimation on our results. Each of the mock observations was sampled from the model and included in the observation set with probability p_det.

In Figure 2, we show values of the odds ratio O_AB where the set of observed events is generated from model A, so that O_AB is expected to increase as N_obs grows. We present results for observations generated from the dynamical (orange) and isolated (blue) models, both for aLIGO/aVirgo (left panels) and ET (right panels).

Zoom In Zoom Out Reset image size

Let us focus first on the top row, where only the total mass is considered in the analysis. In the case in which the dynamical model is assumed to be true (orange) with aLIGO/aVirgo, the upper limit of the 90% credible interval presents high values ${ \mathcal O }\gt {10}^{10}$ even for small N_obs. This is because the isolated model only predicts merging BBHs with M_t < 70 M_⊙. As a result, any observation where the entire support of the total mass posterior distribution is above 70 M_⊙ can only be described by the dynamical model.

Another interesting feature is that the lower bound of the credible interval has values of ${ \mathcal O }\gt {10}^{10}$ starting at ${N}_{\mathrm{obs}}\sim 40$, which is consistent with our models. In fact, as the number of dynamical BBHs with M_t > 70 M_⊙ in our dynamical catalog is equal to 16 (out of 229), the probability of not observing any massive BBHs in 50 observations generated from the dynamical model is equal to p = 1.9 × 10⁻³. It is also interesting to see that the results obtained are similar when observing with ET, suggesting that the differences in total mass range of our models is the dominant effect in the analysis.

In contrast, the evolution of the odds ratio for the isolated case presents a steady increase of ${ \mathcal O }$ with N_obs, as expected. This indicates that, unlike the dynamical case, there is no strong feature in the total mass spectrum of the isolated model that can drive the odds ratio toward very high values with a few observations. In this case, measurements with ET improve the analysis: the median number of observations yielding ${ \mathcal O }={10}^{6}$ (5σ level) decreases from 65 with aLIGO down to ∼30–35.

The second row of Figure 2 presents results obtained when performing the analysis considering both the total mass and the mass ratio. We observe that including the mass ratio helps to discriminate faster between the two models for all cases. This is particularly true for ET, because in the case where the dynamical model is true, the lower bound of the 90% credible interval reaches the 5σ level for N_obs = 20, while we had a value of N_obs = 40 for the analysis considering only the total mass. Similarly, in 90% of the cases, we find that only 30 observations generated from the isolated model and observed by ET give values of odds ratio corresponding to a 5σ level (60 observations were needed when considering only M_t) .

Finally, in the last two rows, we present results obtained by simulating measurements with the total mass, the mass ratio, and the effective spin parameter, adopting either the low-spin (LS, third row) or high-spin (HS, fourth row) models. When the dynamical model is true, only a few observations are needed to push the value of the lower bound of the credible interval toward ${ \mathcal O }\gt {10}^{10}$ (similar results were found in Stevenson et al. (2017a)). This results from our assumption for the spin models, as only the dynamical model has support for χ_eff < 0, meaning that observations of negative values give full support to the dynamical model. As the errors on the spins are already relatively low, with ${\sigma }_{{\chi }_{\mathrm{eff}}}=0.1$ for aLIGO/aVirgo (see Section 3.1.2), we do not see much difference when repeating the analysis with even smaller error for ET. Finally, when the isolated model is true, the 5σ level is attained in 90% of the cases only after ∼15–20 observations for the LS and HS models. Measurements with ET are slightly better and bring the number of necessary observations down to ∼10–15 for both models. Once more, we highlight here that this study was restricted to sources with redshift z ≤ 2, while ET is also expected to see a significant number of sources at high redshift.

We now apply the formalism described in the previous section to aLIGO/aVirgo BBH data from the catalog in Abbott et al. (2019b). We have used both the posterior and prior samples for each event released by the LVC, to compute the expression of the posterior in Equation (17). In Table 2, we present the values obtained for ${ \mathcal O }$ and σ for the favored model under the same parameter configurations presented in the last section.

Table 2.  Results Obtained on the Set of 10 BBHs from Abbott et al. (2019b).

Parameters	Favored Model	${ \mathcal O }$	σ-equivalent
Mt	Dynamical	42	2.3
Mt,q	Dynamical	1.7 × 103	3.4
Mt,q,${\chi }_{{\mathrm{eff}}_{{\rm{L}}}}$	Dynamical	39	2.2
Mt,q,${\chi }_{{\mathrm{eff}}_{{\rm{H}}}}$	Dynamical	4.1 × 1016	x
Mt,q [1]	Dynamical	16	1.9
Mt,q,${\chi }_{{\mathrm{eff}}_{{\rm{H}}}}$ [2]	Dynamical	$2.7\times {10}^{6}$	5.1

Note. We run our statistical analysis varying the set of parameters considered. The configurations [1] and [2] correspond to cases where we respectively exclude GW170729 and the set of events (GW1701014, GW170818, GW170823) from the analysis. Among these models, we prefer a dynamical origin, mainly because of the hard cuts in the mass and spin distributions (see Section 5).

Download table as:  ASCIITypeset image

For all simulations, the values of the odds ratio are ${ \mathcal O }\lesssim {10}^{4}$, corresponding to ≲4σ and favor a dynamical origin. If we only take into account the mass parameters, the analysis done with $\{{M}_{t},q\}$ gives a value of ${ \mathcal O }=1.7\times {10}^{3}$ (corresponding to 3.4σ) in favor of the dynamical model. At the fixed metallicity considered here, our isolated model has a hard cutoff on the total mass at 70 M_⊙, which cannot accommodate the presence of GW170729 (but we stress here that this is not the case for all metallicities). As a result, if we discard this event from the observed catalog, the odds ratio drops to ${ \mathcal O }=1.6\times {10}^{1}$ when considering {M_t, q}.

Our results depend strongly on the underlying spin model. In the LS case, we get ${ \mathcal O }=39$, with a slight preference for the dynamical model. This conclusion changes drastically in the HS case, as the value of the odds ratio is now larger than 10¹⁶. This is because, in the HS case, the dynamical model has a much wider range of possible negative values for χ_eff, extending down to χ_eff − 0.5 (compared to χ_eff − 0.2 for LS case). As some of the events (namely GW1701014, GW170818, GW170823) have a posterior distribution with some support for values of χ_eff < −0.25, the HS dynamical model is significantly favored over the HS isolated case. When removing these events from the analysis, we found that the odds ratio is reduced by 10 orders of magnitude, down to 2.7 × 10⁶.

We generate mock observations assuming a "true" mixing fraction ${{f}}_{T}=0.7$. To generate N_obs events from a mixed model, we first generated N_obs events both for isolated and dynamical models. For each entry of the mixed model, we draw a random number $\epsilon \in { \mathcal U }[0,1]$ and associate an event from the preprocessed isolated (dynamical) set if $\epsilon \lt {{f}}_{T}$ ($\epsilon \gt {{f}}_{T}$).

In Figure 3, we report the values of the medians (square), 90% (straight line), and 99% (dashed line) credible intervals for a set of mock observations with ${N}_{\mathrm{obs}}=\{10,100,500\}$. The values correspond to the current number of events (N_obs = 10), an optimistic prediction for O3 (N_obs = 100) and a high-statistic case (N_obs = 500). For simplicity, we restrict this study to the aLIGO/aVirgo detector case. For each set of observations, we perform the analysis with the combination of parameters $\{{M}_{t}\}$ (orange), $\{{M}_{t},q\}$ (purple), $\{{M}_{t},q,{\chi }_{{\mathrm{eff}}_{{\rm{L}}}}\}$ (green), and $\{{M}_{t},q,{\chi }_{{\mathrm{eff}}_{{\rm{H}}}}\}$ (red).

Zoom In Zoom Out Reset image size

For any set of parameters, we observe a reduction of the width of the credible intervals for higher values of N_obs, as expected. In fact, for N_obs = 10, the 99% credible interval spans almost the entire range of values for ${f}$. For N_obs = 500, this interval is reduced down to 0.2 for M_t, and it narrows to only ∼0.15 when including q and χ_eff in the analysis. Another feature shown in Figure 3 is that the width of the credible interval gets smaller when including more parameters, which is expected as more parameters provide more constraints on the model selection analysis.

In conclusion, Figure 3 suggests that we can already, with 100 detections (an optimistic scenario for the end of O3), constrain the value of the fractional errors on the mixing fraction to an interval smaller than 0.5 using $\{{M}_{t},q\}$. This result is in agreement with previous studies (Vitale et al. 2017; Stevenson et al. 2017a; Zevin et al. 2017; Talbot & Thrane 2017). Furthermore, with even higher statistics of a few hundred detections, the fractional errors on the mixing fraction will go down to 20%, if we consider only the mass parameters. The inclusion of the effective spin parameter reduces this value even further down to 10%. Finally, we have repeated this study for a value of ${{f}}_{T}=0.3$, and found similar predictions.

In this section, we describe the results obtained when applying our mixed model analysis to aLIGO/aVirgo data (Abbott et al. 2019b). Figure 4 shows the posterior distributions obtained when doing the analysis using $\{{M}_{t}\}$ (orange), $\{{M}_{t},q\}$ (purple), $\{{M}_{t},q,{\chi }_{{\mathrm{eff}}_{{\rm{L}}}}\}$ (green), and $\{{M}_{t},q,{\chi }_{{\mathrm{eff}}_{{\rm{H}}}}\}$ (red).

Zoom In Zoom Out Reset image size

First, the two posterior distributions obtained when using the mass parameters are slightly shifted toward values of ${f}\lt 0.5$ (pure dynamical scenario), with median values of 0.40 and 0.36 for {M_t} and {M_t, q} respectively. In addition, the upper limit of the 99% credible interval is equal to 0.94 ({M_t} case) and 0.90 ({M_t, q} case), thus excluding our pure isolated scenario. As before, this is due to the fact that some of the detected BBHs have support for M_t > 70 M_⊙, while our isolated model strictly prevents the occurrence of these high masses. In particular, for GW170729, the lower 90% credible interval of M_t is equal to 74.2 M_⊙. For testing purposes, we ran an analysis where the total mass of GW170729 was not included (blue curve); in this case, the median was 0.55 and the upper limit of the 99% credible interval is much higher, with values of 0.99. As 0.99 is also the 99% upper limit on our prior, this shows that our posterior does not disfavor the isolated scenario any less than the prior.

Including the effective spin parameter in the analysis has a significant impact on the posterior distributions. In the LS model, the distribution of the mixing fraction is centered toward higher values of ${f}$ (with a median of 0.58), favoring the isolated scenario. Once again, this result is dominated by a single event, GW151226, which presents an effective spin parameter with a clearly positive median value χ_eff = 0.18 that is not well-supported by the dynamical scenario. We reran the analysis for the LS model excluding this event (cyan curve) and found that the value of the median is reduced down to a value of 0.48. It is interesting to highlight that, while GW170729 has even higher values of effective spin (median of χ_eff = 0.36), this event gives only limited support to the isolated scenario, as the masses are very high (unlike GW151226). In the HS case, the dynamical scenario is favored with a median value for the mixing fraction equal to 0.27. This can be understood by taking in account the fact that the support for the effective spin parameter in the dynamical case extends between −0.6 and 0.7 (see Figure 1), which is the range of all the events currently observed by aLIGO and aVirgo, while the isolated scenario now struggles to capture events with negative values of χ_eff.

In conclusion, our results suggest that O1+O2 LVC data (Abbott et al. 2019b) exclude a pure isolated scenario as described by our population-synthesis simulations, which include hard cutoffs on both mass and spin distributions. It is important to point out that the width of the posterior distribution is still quite large, and it is in agreement with the results obtained in the last section. Moreover, we stress that the models considered in this paper refer to a single metallicity: at lower progenitor metallicity, even the isolated scenario includes higher-mass BHs (with M_t up to 80–90 M_⊙, Giacobbo & Mapelli 2018). Thus, in a follow-up study (Y. Bouffanais et al. 2019, in preparation), we will explore the importance of metallicity and other population-synthesis parameters.