The first gravitational-wave source from the isolated evolution of two 40-100 Msun stars

References and Notes

The Methods

Our Monte Carlo evolutionary modeling is performed with the StarTrack binary population synthesis code . In particular, we incorporate a calibrated treatment of tidal interactions in close binaries , a physical measure of the common envelope (CE) binding energy , and a rapid explosion supernova model that reproduces the observed mass gap between neutron stars and black holes . Our updated mass spectrum of black holes shows a strong dependence on the metallicity of the progenitor stars (Extended Data Fig. 5). In galaxies with metallicities similar to the Milky Way ($Z={\rm~Z}_{\odot}=0.02$) black holes formed out of single massive stars (initial mass $M_{\rm ZAMS}=150{\rm~M}_{\odot}$) reach a maximum mass of $M_{\rm BH}=15{\rm~M}_{\odot}$, while for very low metallicity ($Z=0.0001=0.5\%{\rm~Z}_{\odot}$) the maximum mass becomes $M_{\rm BH}=94{\rm~M}_{\odot}$. The above input physics represents our standard model (M1), which is representative of our classical formation scheme for double compact objects (BH-BH, BH-NS, and NS-NS).

We have adopted specific values for a number of evolutionary parameters. Single stars are evolved with calibrated formulae based on detailed evolutionary calculations . Massive star winds are adopted from detailed studies of radiation driven mass loss . For the Luminous Blue Variable phase the high mass loss rate is adopted ($1.5\times 10^{-4}{\rm~M}_{\odot}$ yr-1). Binary interactions, and in particular the stability of RLOF, is judged based on binary parameters: mass ratio, evolutionary stage of donor, response to mass loss, and behavior of the orbital separation in response to mass transfer. The orbital separation is additionally affected by gravitational radiation, magnetic braking, and angular momentum loss associated with systemic mass loss. During stable RLOF we assume that half of the mass is accreted onto the companion, while the other half ($1-f_{\rm a}=0.5$) is lost with specific angular momentum ($dJ/dt=j_{\rm loss}[J_{\rm orb}/(M_{\rm don}+M_{\rm acc})](1-f_{\rm a})dM_{\rm RLOF}/dt$ with $j_{\rm loss}=1.0$ ). CE is treated with energy balance with fully effective conversion of orbital energy into envelope ejection ($\alpha=1.0$), while the envelope binding energy for massive stars is calibrated by a parameter $\lambda$ that depends on star radius, mass, and metallicity. For massive stars $\lambda\approx 0.1$ is adopted . During CE compact objects accrete at $10\%$ Bondi-Hoyle rate as estimated by recent hydrodynamical simulations . Our CE is done “instantaneously”, so the time at the beginning and end of CE is exactly the same (see Fig. 1); the time duration of CE has no impact on our results.

We consider two extra variations of the binary evolution input physics. In one model (M2) we test highly uncertain CE physics and we allow for Hertzsprung gap (HG) stars to initiate and survive CE evolution. This is an optimistic assumption, since these stars may not allow for CE evolution , nor survive as a binary if CE does happen . For comparison, in our standard model we allow only evolved stars with a deep convective envelope (core Helium burning stars) to survive CE.

In the opposite extreme, we employ a model (M3) where black holes receive high natal kicks. In particular, each BH gets a natal kick with its components drawn from a $1$-D Maxwellian distribution with $\sigma=265$ km s-1, independent of BH mass. Such high natal kicks are measured for Galactic pulsars . This is a pessimistic assumption, as high natal kicks tend to disrupt BH-BH progenitor binaries. This assumption is not yet excluded based on electromagnetic observations . In contrast, in our standard model BH natal kicks decrease with BH mass. In particular, for massive BHs that form through direct collapse of an entire star to a BH with no supernova explosion ($M_{\rm BH}\gtrsim 10{\rm~M}_{\odot}$ for solar metallicity; $M_{\rm BH}\gtrsim 15{\rm~M}_{\odot}$ for $Z=10\%{\rm~Z}_{\odot}$; and $M_{\rm BH}\gtrsim 15$–$30{\rm~M}_{\odot}$ for $Z=1\%{\rm~Z}_{\odot}$) we assume no natal kicks . We have also calculated a series of models with intermediate BH kicks (see Extended Data Fig. 4): $\sigma=200$ km s-1 (model M4), $\sigma=130$ km s-1 (model M5), $\sigma=70$ km s-1 (model M6).

For each evolutionary model we compute $2\times 10^{7}$ massive binaries for each point on a grid of $32$ sub-models covering a wide range of metallicities: $Z=0.0001$, $0.0002$, $0.0003$, $0.0004$, $0.0005$, $0.0006$, $0.0007$, $0.0008$, $0.0009$, $0.001$, $0.0015$, $0.002$, $0.0025$, $0.003$, $0.0035$, $0.004$, $0.0045$, $0.005$, $0.0055$, $0.006$, $0.0065$, $0.007$, $0.0075$, $0.008$, $0.0085$, $0.009$, $0.0095$, $0.01$, $0.015$, $0.02$, $0.025$, $0.03$. We assume that stellar evolution at even lower metallicities proceeds in the same way as the evolution at $Z=0.5\%{\rm~Z}_{\odot}$. However, we note that stars with very low metal content (e.g., Population III) may evolve differently than metal-rich stars .

Each sub-model is computed with initial distributions of orbital periods ($\propto(\log P)^{-0.5}$), eccentricities ($\propto e^{-0.42}$), and mass ratios ($\propto q^{0}$) appropriate for massive stars . We adopt an initial mass function that is close to flat for low mass stars ($\propto M^{-1.3}$ for $0.08\leq M<0.5{\rm~M}_{\odot}$ and $\propto M^{-2.2}$ for $0.5\leq M<1.0{\rm~M}_{\odot}$) and top heavy for massive stars ($\propto M^{-2.3}$ for $1.0\leq M\leq 150{\rm~M}_{\odot}$), as guided by recent observations . The adopted IMF generates higher BH-BH merger rate densities as compared with the steeper IMF ($\propto M^{-2.7}$ for $1.0\leq M\leq 150{\rm~M}_{\odot}$) adopted in our earlier studies as there are more BH-BH merger progenitors in our simulations .

A moderate binary fraction ($f_{\rm bi}=0.5$) is adopted for stars with mass $M_{\rm ZAMS}<10{\rm~M}_{\odot}$, while we assume that all more massive stars are formed in binaries ($f_{\rm bi}=1.0$) as indicated by recent empirical estimates ).

We adopt an extinction corrected cosmic star formation rate based on numerous multi-wavelength observations :

This SFR declines rapidly at high redshifts ($z>2$). This may be contrasted with some SFR models that we have used in the past which generated a greater number of stars at high redshifts. This revision will thus reduce the BH-BH merger rate densities at all redshifts. Even though the formation of BH-BH binaries takes a very short time ($\sim 5$ Myr), the time to coalescence of two black holes may take a very long time (Fig. 1 and Extended Data Fig. 2).

In our new treatment of chemical enrichment of the Universe we follow the mean metallicity increase with cosmic time (since Big Bang till present). The mean metallicity as a function of redshift is given by

with a return fraction $R=0.27$ (mass fraction of each generation of stars that is put back into the interstellar medium), a net metal yield $y=0.019$ (mass of new metal created and ejected into the interstellar medium by each generation of stars per unit mass locked in stars), a baryon density $\rho_{\rm b}=2.77\times 10^{11}\,\Omega_{\rm b}\,h_{0}^{2}\,{\rm~M}_{\odot}\,{\rm Mpc}^{-3}$ with $\Omega_{\rm b}=0.045$ and $h_{0}=0.7$, a star formation rate given by eq. 1, and $E(z)=\sqrt{\Omega_{\rm M}(1+z)^{3}+\Omega_{\rm k}(1+z)^{2}+\Omega_{\Lambda})}$ with $\Omega_{\Lambda}=0.7$, $\Omega_{\rm M}=0.3$, $\Omega_{\rm k}=0$, and $H_{0}=70.0\,{\rm km}\,{\rm s}^{-1}\,{\rm Mpc}^{-1}$. The shape of the mean metallicity dependence on redshift follows recent estimates , although the level was increased by $0.5$ dex to better fit observational data . At each redshift we assume a log–normal distribution of metallicity around the mean, with $\sigma=0.5$ dex . Our prescription (Extended Data Fig. 6) produces more low-metallicity stars than previously . Since BH-BH formation is enhanced at low-metallicity , our new approach increases the predicted rate densities of BH-BH mergers.

Here we discuss caveats of evolutionary calculations. First, we only consider isolated binary evolution, and thus our approach is applicable to field stars in low density environments. It is possible that dynamical interactions enhance BH-BH merger formation in dense globular clusters , offering a completely independent channel.

Second, our predictions are based on a “classical” theory of stellar and binary evolution for the modeling of massive stars that we have compiled, developed, and calibrated over the last $15$ years. We do not consider exotic channels for the formation of BH-BH mergers, such as the one from rapidly rotating stars in contact binaries .

Third, our modeling includes only three evolutionary models: a “standard” model consisting of our best estimates for reasonable parameters, as well as “optimistic” and “pessimistic” alternate models. The optimistic model consists of only one change from the standard model: we allow all stars beyond the main sequence to survive the common envelope phase. Alternatively, the pessimistic model also consists of only one change: larger BH natal kicks. We have not investigated other possible deviations from the standard model (e.g., different assumptions of mass and angular momentum loss during stable mass transfer evolution) nor have we checked inter-parameter degeneracies (e.g., models with high BH kicks and an optimistic common envelope phase). Albeit with low statistics and limited scope, precursor versions of these computationally demanding studies have already been performed ; these calculations indicate that our three models are likely to cover the range of interesting effects.

Fourth, our observations are severely statistically limited. We are attempting to draw inferences about our models based on a single detection (GW150914).

In was argued that the formation of GW150914 in isolated binary evolution requires a metallicity lower than $50\%{\rm~Z}_{\odot}$. This was based on single stellar models ; stars in close binaries are subject to significant mass loss during RLOF/CE, and they form BHs with lower mass than single stars. Thus in binaries the metallicity threshold for massive BH formation is lower than in single stellar evolution. For example, formation of a single $30{\rm~M}_{\odot}$ BH requires $Z<25\%{\rm~Z}_{\odot}$ (Extended Data Fig. 5), while formation of two such BHs in a binary requires $Z<10\%{\rm~Z}_{\odot}$ (Extended Data Fig. 1). The value of this threshold depends on assumptions for the model of stellar evolution, winds, and BH formation processes. The physical models we have adopted yield a threshold of $Z<10\%{\rm~Z}_{\odot}$, the same as the one obtained with MESA for homogeneous stellar evolution . Our model was calibrated on known masses of BHs, and in particular we do not exceed $15{\rm~M}_{\odot}$ for ${\rm~Z}_{\odot}$ (the highest mass stellar BH known in our Galaxy). In contrast, single stellar models used to derive the high metallicity threshold produce $25{\rm~M}_{\odot}$ for ${\rm~Z}_{\odot}$ . The highest threshold obtained with binary evolution was reported at the level of $50\%{\rm~Z}_{\odot}$ . Such high value of metallicity threshold for the progenitor of GW150914 implies that stars at approximately solar metallicity ($Z=0.014$) produce BHs as massive as $40{\rm~M}_{\odot}$ . This is neither supported nor excluded by available electromagnetic BH mass measurements .

In the following we present calculation of the gravitational radiation signal. The output of StarTrack is a binary merger at a given time. We then calculate the gravitational waveform associated with this merger, and determine whether this binary would have been observable by LIGO in the O1 configuration .

We model the full inspiral-merger-ringdown waveform of the binaries using the IMRPhenomD gravitational waveform template family . This is a simple and fast waveform family which neglects the effects of spin (which are not relevant for GW150914). We consider a detection to be given by a threshold SNR $>8$ in a single detector, and we use the fiducial O1 noise curve . We calculate the face-on, overhead SNR for each binary directly from Eq. 2 of . We then calculate the luminosity distance at which this binary would be detected with $\mbox{SNR}=8$. Note that as the distance to the binary changes, the observer frame (redshifted) mass also changes, and therefore calculating the horizon redshift requires an iterative process. Once this has been calculated, we then determine the predicted detection rates using Eq. 9 of . The effects of the antenna power pattern are incorporated in the $p_{\rm det}$ term in this equation.

Estimate of fiducial aLIGO sensitivity during the 16-day GW150914 analysis is shown in Figure 3. We estimate the sensitivity to coalescing compact binaries using a reference O1 noise curve. We assume that both detectors operate with the fiducial O1 noise curve, which is the same sensitivity we adopted to calculate compact binary detection rates. For comparison, this model agrees reasonably well with the “early-high” sensitivity model . Our expression is a 50th percentile upper limit, assuming no detections. The critical application of this expression is not related to its overall normalization; we are instead interested in its shape, which characterizes the strongly mass-dependent selection biases of LIGO searches.

Using these inputs, our fiducial estimate of the advanced LIGO sensitivity during the first 16 days of O1 for a specific mass bin, $\Delta M_{\rm i}$, is given by

where $T$ is 16 days, corresponding to the analysis of GW150914 , and the volume

is the sensitive volume averaged over the mass bin $\Delta M_{\rm i}$, and $p_{\rm det}(<w,M)$ is the orientation-averaged detection probability . The function $p_{\rm det}(<w,M)$ depends on the coalescing binary redshifted mass through the maximum luminosity distance (“horizon distance”) at which a source could produce a response of SNR$>$8 in a single detector. To calculate this distance, we adopt the same model for inspiral, merger, and ringdown used in the text to estimate compact binary detection rates. Extended Data Figure 7 shows our estimated horizon redshift as a function of the total redshifted binary merger mass for equal mass mergers.

Code availability. We have opted not to release the population synthesis code StarTrack used to generate binary populations for this study.

References and Notes