\section*{The Methods}

Our Monte Carlo evolutionary modeling is performed with the {\tt StarTrack} binary 
population synthesis code~\cite{Belczynski2002}. In particular, we incorporate a 
calibrated treatment of tidal interactions in close binaries~\cite{Belczynski2008a}, 
a physical measure of the common envelope (CE) binding energy~\cite{Dominik2012,Xu2010}, 
and a rapid explosion supernova model that reproduces the observed mass gap between 
neutron stars and black holes~\cite{Fryer2012,Belczynski2012}. Our updated mass spectrum 
of black holes shows a strong dependence on the metallicity of the progenitor stars 
(Extended Data Fig.~\ref{fig:BHMassRel}). 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~\cite{Hurley2000}. Massive star winds are adopted from 
detailed studies of radiation driven mass loss~\cite{Vink2011}. 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$~\cite{Podsiadlowski1992}). 
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~\cite{Xu2010}. 
During CE compact objects accrete at $10\%$ Bondi-Hoyle rate as estimated by recent 
hydrodynamical simulations~\cite{Ricker2008,Macleod2015}. 
Our CE is done ``instantaneously'', so the time at the beginning
and end of CE is exactly the same (see Fig.~\ref{fig:ExamplEvol}); 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~\cite{Ivanova2013} 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~\cite{Pavlovskii2015}, nor survive as a binary if CE does 
happen~\cite{Belczynski2007}.
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~\cite{Hobbs2005}. 
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~\cite{Belczynski2015}.  
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~\cite{Fryer2012}.
We have also calculated a series of models with intermediate BH kicks (see Extended Data 
Fig.~\ref{fig:IntrRate}):
\mbox{$\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~\cite{Szecsi2015}.

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~\cite{Sana2012}.
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~\cite{Bastian2010}. 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~\cite{Dominik2013,Belczynski2015} as there are more BH-BH merger 
progenitors in our simulations~\cite{deMink2015}.

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~\cite{Sana2012,Duchene2013}).

We adopt an extinction corrected cosmic star formation rate based on numerous 
multi-wavelength observations~\cite{Madau2014}:
\begin{equation}
\mbox{SFR}(z)=0.015 {(1+z)^{2.7} \over 1+[(1+z)/2.9]^{5.6}} \,{\rm ~M}_{\odot} \, {\rm Mpc}^{-3} 
\, {\rm yr}^{-1}.
\label{sfr}
\end{equation} 
This SFR declines rapidly at high redshifts ($z>2$). This may be contrasted with some 
SFR models that we have used in the past~\cite{Strolger2004} which generated a greater 
number of stars at high redshifts. This revision will thus reduce the BH-BH merger 
rate densities at {\em 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.~\ref{fig:ExamplEvol} and Extended Data 
Fig.~\ref{fig:Emergence}).

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
\begin{equation}
\log(Z_{\rm mean}(z))=0.5+\log \left(  {y \, (1-R) \over \rho_{\rm b}} 
\int_z^{20} {97.8\times10^{10} \, sfr(z') \over H_0 \, E(z') \, (1+z')} dz' \right)
\label{Zmean}
\end{equation}
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.~\ref{sfr}, 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~\cite{Madau2014}, although the level was increased by $0.5$ dex to better 
fit observational data~\cite{Vangioni2015}.
At each redshift we assume a log--normal distribution of metallicity around the 
mean, with $\sigma=0.5$ dex~\cite{Dvorkin2015}. 
Our prescription (Extended Data Fig.~\ref{fig:MetEvol}) produces more low-metallicity 
stars than previously~\cite{Dominik2013}. Since BH-BH formation is enhanced at
low-metallicity~\cite{Belczynski2010a}, 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~\cite{Rodriguez2016}, 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~\cite{Almeida2015}. 

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 {\em and}\/ an optimistic common envelope phase). 
Albeit with low statistics and limited scope, precursor versions of these 
computationally demanding studies have already been performed~\cite{O'Shaughnessy2007}; 
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~\cite{AstroPaper} 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~\cite{Spera2015}; 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 {\em single}\/ $30{\rm ~M}_{\odot}$ 
BH requires $Z<25\%{\rm ~Z}_{\odot}$ (Extended Data Fig.~\ref{fig:BHMassRel}), while formation 
of two such BHs in a {\em binary} requires $Z<10\%{\rm ~Z}_{\odot}$ (Extended Data 
Fig.~\ref{fig:BHBHMass}). 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~\cite{Marchant2016}.
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}$~\cite{Spera2015}. The highest threshold obtained 
with binary evolution was reported at the level of $50\%{\rm ~Z}_{\odot}$~\cite{Eldridge2016}. 
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}$~\cite{Eldridge2016}. This is neither supported nor excluded by available 
electromagnetic BH mass measurements~\cite{fn2}. 

In the following we present calculation of the gravitational radiation signal.  
The output of {\tt 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~\cite{Dominik2015,Belczynski2015}. 

We model the full inspiral-merger-ringdown waveform of the binaries using the 
IMRPhenomD gravitational waveform template family~\cite{Khan2015,Husa2015}. 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~\cite{fn1}. 
We calculate the face-on, overhead SNR for each binary directly from Eq.~2 of
\cite{Dominik2015}. 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~\cite{Dominik2015}. 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~\ref{fig:RateComp}. 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~\cite{LVC2013}. 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
\begin{eqnarray}
R_{{\rm D},\Delta M_{\rm i}, UL} = \frac{0.7}{V_{\Delta M_{\rm i}}T} ,
\end{eqnarray}
where $T$ is 16 days, corresponding to the analysis of GW150914~\cite{DiscoveryPaper}, 
and the volume 
$$
V_{\Delta M_{\rm i}} = (\Delta M_{\rm i})^{-1}\int dM \int \frac{dz}{1+z} \frac{dV}{dz} 
p_{\rm det}(<w,M)
$$
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~\cite{Dominik2015,Belczynski2015}.  
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~\cite{Husa2015,Khan2015} used 
in the text to estimate compact binary detection rates. Extended Data 
Figure~\ref{fig:Horizon} shows our estimated horizon redshift as a function of the 
total redshifted binary merger mass for equal mass mergers.


\noindent
{\bf Code availability.}
We have opted not to release the population synthesis code {\tt StarTrack} 
used to generate binary populations for this study.


