Binary black hole merger rate constraints using GWTC-3 and full-O3 stochastic background constraints

<h1>README</h1> <p>This dataset contains posterior measurements of the redshift-dependent merger&nbsp;rate, mass distribution, and spin distribution of binary black holes following&nbsp;the O3b observing run of the LIGO-Virgo-KAGRA network, including both direct compact binary detections and constraints on the astrophysical&nbsp;gravitational-wave background.</p> <p>In particular, the goal of this work is to constrain a more complex model for the black hole merger rate, with the comoving rate density evolving as</p> <p>$$<br>R(z) \propto \frac{(1+z)^\alpha}{1 + \left(\frac{1+z}{1+z_p}\right)^{\alpha + \beta}}.<br>$$</p> <p>At redshifts \(z &lt; z_p\), the merger rate grows approximately as \(R(z) \propto (1+z)^\alpha\), whereas at \(z&gt;z_p\) it falls as \(R(z) \propto (1+z)^{-\beta}\).</p> <p>The analysis was performed as described in <a href="https://iopscience.iop.org/article/10.3847/2041-8213/ab9743">Callister <em>et al</em> (2020)</a> and <a href="https://link.aps.org/doi/10.1103/PhysRevD.104.022004">Abbott&nbsp;<em>et al</em> (2021)</a>, now using binary black&nbsp;hole detections from the GWTC-3 catalog (<a href="https://link.aps.org/doi/10.1103/PhysRevX.13.041039">Abbott <em>et al</em> 2023a</a>, <a href="https://link.aps.org/doi/10.1103/PhysRevX.13.011048">2023b</a>).</p> <ul> <li>The parameter estimation samples used are those provided by the LIGO-Virgo-KAGRA collaboration at https://zenodo.org/records/8177023</li> <li>Selection effects are calculated and mitigated using the suite of pipeline&nbsp;injections available at https://zenodo.org/records/7890398</li> <li>Cross-correlation measurements of the stochastic gravitational-wave background&nbsp;are available at https://dcc.ligo.org/LIGO-G2001287, and correspond to&nbsp;results presented in <a href="https://link.aps.org/doi/10.1103/PhysRevD.104.022004">Abbott <em>et al</em> (2021).</a></li> </ul> <p>As discussed in <a href="https://link.aps.org/doi/10.1103/PhysRevX.13.011048">Abbott <em>et al</em> (2023b)</a>, the results of this combined BBH + stochastic analysis are categorically&nbsp;unchanged relative to results previously obtained using GWTC-2 events (<a href="https://link.aps.org/doi/10.1103/PhysRevD.104.022004">Abbott <em>et al</em> 2021</a>); sensitivities&nbsp;are not yet sufficient to resolve the redshift at which the black hole merger&nbsp;rate peaks and turns over.</p> <h1>Contents</h1> <ul> <li><code><strong>processed_emcee_samples_together_r00r01.npy</strong></code>: File containing posterior samples when jointly analyzing BBH detections and stochastic background upper limits.</li> <li><code><strong>processed_emcee_samples_noStochastic_r00r01.npy</strong></code>: File containing posterior samples analyzing only direct BBH detections.</li> <li><code><strong>run_emcee_plPeak.py</strong></code>: Script performing joint hierarchical inference using BBH detections and stochastic background limits; used to generate posterior samples in <code>processed_emcee_samples_together_r00r01.npy</code></li> <li><code><strong>run_emcee_plPeak_noStochastic.py</strong></code>: Script performing joint hierarchical inference using BBH detections and stochastic background limits; used to generate posterior samples in <code>processed_emcee_samples_noStochastic_r00r01.npy</code></li> </ul> <h1>Accessing posterior samples</h1> <p>Posterior samples are contained in the files <code>processed_emcee_samples_together_r00r01.npy</code> and <code>processed_emcee_samples_noStochastic_r00r01.npy</code>. This is loaded via python as, e.g.</p> <blockquote> <p>&gt;&gt;&gt; import numpy as np</p> <p>&gt;&gt;&gt; dataset = np.load('processed_emcee_samples_together_r00r01.npy')</p> </blockquote> <p>Contained in this file is a single <code>numpy</code> array of size <code>(# of posterior samples, # of hyperparameters)</code>:</p> <blockquote> <p>&gt;&gt;&gt; dataset.shape</p> <p>(1152, 13)</p> </blockquote> <p>&nbsp;</p> <p>The 13 hyperparameters are defined as follows:</p> <table> <tbody> <tr> <td>Column</td> <td>Name</td> <td>Definition</td> </tr> <tr> <td><code>dataset[:, 0]</code></td> <td><code>xeff_mu</code></td> <td>Mean effective inspiral spin</td> </tr> <tr> <td><code>dataset[:, 1]</code></td> <td><code>xeff_sig</code></td> <td>Standard deviation of effective inspiral spin</td> </tr> <tr> <td><code>dataset[:, 2]</code></td> <td><code>R0</code></td> <td>Total BBH merger rate at redshift \(z=0\)</td> </tr> <tr> <td><code>dataset[:, 3]</code></td> <td><code>mMin</code></td> <td>Minimum black hole mass</td> </tr> <tr> <td><code>dataset[:, 4]</code></td> <td><code>mMax</code></td> <td>Maximum black hole mass</td> </tr> <tr> <td><code>dataset[:, 5]</code></td> <td><code>lmbda</code></td> <td>Power-law index on primary mass distribution</td> </tr> <tr> <td><code>dataset[:, 6]</code></td> <td><code>mu_peak</code></td> <td>Mean of Gaussian peak in primary mass distribution</td> </tr> <tr> <td><code>dataset[:, 7]</code></td> <td><code>sig_peak</code></td> <td>Standard deviation of Gaussian peak</td> </tr> <tr> <td><code>dataset[:, 8]</code></td> <td><code>frac_peak</code></td> <td>Fraction of events occupying Gaussian peak</td> </tr> <tr> <td><code>dataset[:, 9]</code></td> <td><code>bq</code></td> <td>Power-law index on mass ratio distribution \(p(q\|m_1)\)</td> </tr> <tr> <td><code>dataset[:, 10]</code></td> <td><code>alpha</code></td> <td>Slope of \(R(z) \propto (1+z)^\alpha \) at low redshifts</td> </tr> <tr> <td><code>dataset[:, 11]</code></td> <td><code>beta</code></td> <td>Slope of \(R(z) \propto (1+z)^{-\beta}\) at high redshifts</td> </tr> <tr> <td><code>dataset[:, 12]</code></td> <td><code>zpeak</code></td> <td>Redshift at which \(R(z)\) peaks and turns over</td> </tr> </tbody> </table> <p>&nbsp;</p> <p>The exact usage of the above parameters can be seen in the included scripts <code>run_emcee_plPeak.py</code> and <code>run_emcee_plPeak_noStochastic.py</code>, with which the inference was performed.&nbsp;</p>