Advanced - GW190521 BH1 as 3rd gen BH
In this notebook, we perform a quick demonstration of using archeo to analyze the GW190521 event, assuming the primary black hole is a third generation black hole (yield from a previous merge of one second generation black hole and one first generation black hole). Again, similar to gw190521_quick.ipynb, we will use a subsampled version of the LVK posterior samples for GW190521 as our input data.
The full dataset can be found dcc.ligo.
The notebook run should take less than 5 minute to complete.
import pandas as pd
import archeo
import warnings
warnings.filterwarnings("ignore", category=UserWarning)
warnings.filterwarnings("ignore", category=FutureWarning)
warnings.filterwarnings("ignore", category=RuntimeWarning)
warnings.filterwarnings("ignore", category=SyntaxWarning)
Prepare Prior¶
For the 2g1g analysis, we will assume all hierarchical mergers undergo a generice spin setting including precession, similar to PPq6 in Araújo Álvarez et al., 2024. Note that we assume the first generation black hole mass range is in [5, 65] solar masses based on PISN gap in [65, 130] solar masses. You may feel free to modify the upper limit based on your own assumptions.
pipeline_1g1g = archeo.get_binary_generation_pipeline("2g_aligned_spin")
df_1g1g_binaries, _ = pipeline_1g1g(size=5000, n_workers=1)
pipeline_2g1g = archeo.get_binary_generation_pipeline("ng_aligned_spin")
df_2g1g_binaries, binary_generator = pipeline_2g1g(df_bh1_binaries=df_1g1g_binaries, size=5000, n_workers=1)
Loaded NRSur3dq8Remnant fit.
100%|██████████| 5000/5000 [00:06<00:00, 829.08it/s]
Loaded NRSur3dq8Remnant fit.
100%|██████████| 5000/5000 [00:06<00:00, 777.50it/s]
Verify Distribution in Priors¶
Here we quickly visualize the remnant mass and spin distributions of the 2g prior and the primary black hole mass and spin distributions of the 1g prior to ensure they are as expected.
We will also visualize the remnant mass and spin distributions of the 2g1g prior to make a sense of what is happening.
import matplotlib.pyplot as plt
from archeo.visualization.base import plot_pdf
# Verify Distribution in Priors
# Here the distributions should be fairly close to each other
# since we are using the 2g prior samples
# for the primary black hole in the 2g1g prior.
# However, the sample size is small so there will be some noise.
fig, axes = plt.subplots(1, 2, figsize=(10, 5))
plot_pdf(axes[0], df_1g1g_binaries["m_f"], "blue", name="2g Remnant Mass", unit="M_\odot")
plot_pdf(axes[0], df_2g1g_binaries["m_1"], "lightblue", name="2g1g BH1 Mass", unit="M_\odot")
plot_pdf(axes[1], df_1g1g_binaries["a_f"], "blue", name="2g Remnant Spin")
plot_pdf(axes[1], df_2g1g_binaries["a_1"], "lightblue", name="2g1g BH1 Spin")
for i in range(2):
axes[i].set_title(
"{param} Distribution Comparison".format(param="Mass" if i == 0 else "Spin")
)
axes[i].set_ylabel("Probability Density")
axes[i].set_xlabel("Mass (M_\odot)" if i == 0 else "Spin")
axes[i].grid()
axes[i].legend()
# Understand the remnant mass and spin distributions of the 2g1g prior
# Here we should at least see the remnant mass distribution shifted to higher masses
# compared to the 2g prior since the primary black hole is more massive on average.
fig, axes = plt.subplots(1, 2, figsize=(10, 5))
plot_pdf(axes[0], df_1g1g_binaries["m_f"], "blue", name="2g Remnant Mass", unit="M_\odot")
plot_pdf(axes[0], df_2g1g_binaries["m_f"], "lightblue", name="2g1g Remnant Mass", unit="M_\odot")
plot_pdf(axes[1], df_1g1g_binaries["a_f"], "blue", name="2g Remnant Spin")
plot_pdf(axes[1], df_2g1g_binaries["a_f"], "lightblue", name="2g1g Remnant Spin")
for i in range(2):
axes[i].set_title(
"{param} Distribution Comparison".format(param="Mass" if i == 0 else "Spin")
)
axes[i].set_ylabel("Probability Density")
axes[i].set_xlabel("Mass (M_\odot)" if i == 0 else "Spin")
axes[i].grid()
axes[i].legend()
Load Parameter Estimation Samples¶
Here the steps are identical to those in gw190521_quick.ipynb.
# This data file is a subset of the LVK posterior samples for GW190521
# See https://github.com/wyhwong/archeo/tree/main/src/tests/test_data/gw190521_lvk_subsampled.json
# For the full dataset, please refer to: https://dcc.ligo.org/P2000158/public
gw190521_pe_samples = pd.read_json("./gw190521_lvk_subsampled.json")
gw190521_pe_samples.head(10)
| mass_1_source | mass_2_source | a_1 | a_2 | |
|---|---|---|---|---|
| 0 | 84.513663 | 76.297304 | 0.687007 | 0.899445 |
| 1 | 81.059729 | 69.663127 | 0.769092 | 0.878550 |
| 2 | 92.370429 | 86.105830 | 0.540142 | 0.949552 |
| 3 | 77.363091 | 67.721241 | 0.088123 | 0.973257 |
| 4 | 99.454443 | 41.802952 | 0.710768 | 0.976182 |
| 5 | 89.774109 | 72.749785 | 0.096520 | 0.859641 |
| 6 | 65.693613 | 63.273635 | 0.664623 | 0.867052 |
| 7 | 85.314123 | 58.547913 | 0.604760 | 0.867737 |
| 8 | 88.031480 | 85.329605 | 0.741590 | 0.659539 |
| 9 | 71.147886 | 65.844700 | 0.914395 | 0.073000 |
ancestors_bh1 = archeo.infer_ancestral_posterior_distribution(
df_binaries=df_2g1g_binaries,
mass_posterior_samples=gw190521_pe_samples["mass_1_source"].values.tolist(),
spin_posterior_samples=gw190521_pe_samples["a_1"].values.tolist(),
n_workers=1,
)
ancestors_bh1.head(10)
100%|██████████| 5000/5000 [00:05<00:00, 974.11it/s]
| m_1 | a_1 | v_1 | m_2 | a_2 | v_2 | m_f | a_f | k_f | chi_eff | chi_p | q | logL | spin_measure | mass_measure | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 65.831050 | 0.353322 | 365.279865 | 22.124244 | 0.609223 | 0.000000 | 84.633698 | 0.712134 | 103.448526 | 0.417691 | 0.0 | 2.975516 | -5.744604 | 0.687007 | 84.513663 |
| 1 | 53.119709 | 0.802843 | 140.728853 | 32.413284 | 0.984165 | 0.000000 | 81.036127 | 0.783052 | 251.992394 | 0.125645 | 0.0 | 1.638825 | -5.472671 | 0.769092 | 81.059729 |
| 2 | 58.394803 | 0.220110 | 0.000000 | 36.909588 | 0.645634 | 24.305592 | 91.823105 | 0.544498 | 93.547576 | -0.384907 | 0.0 | 1.582104 | -5.683980 | 0.540142 | 92.370429 |
| 3 | 63.264025 | 0.750997 | 124.797344 | 15.707040 | 0.426390 | 0.000000 | 77.766324 | 0.053967 | 233.438467 | -0.686434 | 0.0 | 4.027750 | -7.418581 | 0.088123 | 77.363091 |
| 4 | 55.884501 | 0.444157 | 61.081178 | 48.053683 | 0.214754 | 0.000000 | 98.456368 | 0.736262 | 117.553609 | 0.139523 | 0.0 | 1.162960 | -5.259097 | 0.710768 | 99.454443 |
| 5 | 75.477138 | 0.898486 | 22.410608 | 14.883065 | 0.345397 | 0.000000 | 89.288024 | 0.129097 | 207.617190 | -0.693609 | 0.0 | 5.071344 | -7.824046 | 0.096520 | 89.774109 |
| 6 | 48.383650 | 0.276768 | 147.457452 | 20.130373 | 0.443253 | 0.000000 | 65.641566 | 0.709175 | 115.004725 | 0.325684 | 0.0 | 2.403515 | -6.907755 | 0.664623 | 65.693613 |
| 7 | 54.973616 | 0.016539 | 0.000000 | 34.110196 | 0.456059 | 182.718164 | 85.512721 | 0.613566 | 83.928041 | -0.184832 | 0.0 | 1.611648 | -6.119298 | 0.604760 | 85.314123 |
| 8 | 55.529923 | 0.008108 | 0.000000 | 37.793997 | 0.879787 | 29.857368 | 88.257144 | 0.745623 | 185.899460 | 0.361118 | 0.0 | 1.469279 | -5.521461 | 0.741590 | 88.031480 |
| 9 | 54.706320 | 0.748035 | 0.000000 | 21.524488 | 0.544730 | 303.951564 | 71.691424 | 0.873448 | 26.528114 | 0.690630 | 0.0 | 2.541585 | -5.878136 | 0.914395 | 71.147886 |
Visualize Inferred Posteriors of BH1's Ancestor¶
Here we just visualize the summary for convenience.
## Visualize Inferred Posteriors of BH1's Ancestor
from archeo.visualization.estimation import table_estimates
table_estimates(
dfs={
"Prior (Aligned Spin)": df_2g1g_binaries,
"GW190521 BH1 Ancestor": ancestors_bh1,
}
)
| Recovery Rate | p2g_GC | p2g_MW | p2g_NSC | p2g_EG | $m_1$ | $m_2$ | $q$ | $m_f$ | $a_f$ | $v_f$ | $\chi_{p}$ | $\chi_{eff}$ | ||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | Prior (Aligned Spin) | 1.0000 | 7.88 | 46.70 | 46.70 | 46.70 | $67.04_{-32.64}^{+37.66}$ | $34.68_{-21.86}^{+25.94}$ | $1.86_{-0.79}^{+3.05}$ | $99.57_{-47.87}^{+46.07}$ | $0.66_{-0.54}^{+0.23}$ | $149.91_{-127.55}^{+209.84}$ | $0.00_{-0.00}^{+0.00}$ | $0.02_{-0.75}^{+0.71}$ |
| 1 | GW190521 BH1 Ancestor | 0.9838 | 14.60 | 61.74 | 61.74 | 61.74 | $60.98_{-17.23}^{+20.26}$ | $27.89_{-14.57}^{+17.91}$ | $2.20_{-1.12}^{+3.20}$ | $84.95_{-14.03}^{+20.35}$ | $0.69_{-0.63}^{+0.24}$ | $128.81_{-115.08}^{+205.26}$ | $0.00_{-0.00}^{+0.00}$ | $0.12_{-0.89}^{+0.74}$ |