Customize Simulation

In this notebook, we demonstrate how to customize a simulation of black hole binary populations by defining your own distributions and domains. This allows you to explore a wide range of astrophysical scenarios and test the sensitivity of your inference results to different assumptions about the underlying population.

The notebook run should take less than 5 minutes to compute.

In [1]:

Copied!





from archeo.constants.enum import Fits
from archeo.data_structures.distribution import Uniform
from archeo.data_structures.math import Domain
from archeo.data_structures.physics.binary import BinaryGenerator
from archeo.data_structures.physics.black_hole import BlackHoleGenerator
from archeo.data_structures.physics.mahapatra import MahapatraMassFunction
from archeo.postprocessing.dataframe import convert_simulated_binaries_to_dataframe
from archeo.simulation.simulate_merger import simulate_black_hole_mergers
from archeo.constants.enum import Fits
from archeo.data_structures.distribution import Uniform
from archeo.data_structures.math import Domain
from archeo.data_structures.physics.binary import BinaryGenerator
from archeo.data_structures.physics.black_hole import BlackHoleGenerator
from archeo.data_structures.physics.mahapatra import MahapatraMassFunction
from archeo.postprocessing.dataframe import convert_simulated_binaries_to_dataframe
from archeo.simulation.simulate_merger import simulate_black_hole_mergers

In [2]:

Copied!





import warnings

warnings.filterwarnings("ignore", category=UserWarning)
warnings.filterwarnings("ignore", category=FutureWarning)
warnings.filterwarnings("ignore", category=RuntimeWarning)
warnings.filterwarnings("ignore", category=SyntaxWarning)
import warnings

warnings.filterwarnings("ignore", category=UserWarning)
warnings.filterwarnings("ignore", category=FutureWarning)
warnings.filterwarnings("ignore", category=RuntimeWarning)
warnings.filterwarnings("ignore", category=SyntaxWarning)

Here we first define a generator of black hole, where we specify the distribution of parameters such as mass and spin. In this example, we use the MahapatraMassFunction for the mass distribution to illustrate the use of an astrophysics motivated mass function, and a uniform distribution for the spin magnitude. You can easily swap in different distribution or customize the parameters to explore various population models.

In [3]:

Copied!





bh_generator = BlackHoleGenerator(
    mass_distribution=MahapatraMassFunction(mass=Domain(low=5, high=65)),
    spin_magnitude_distribution=Uniform(low=0, high=1),
)
bh_generator = BlackHoleGenerator(
    mass_distribution=MahapatraMassFunction(mass=Domain(low=5, high=65)),
    spin_magnitude_distribution=Uniform(low=0, high=1),
)

We then define a generator of binary systems, where we specify the range of the mass ratio and whether the spins are aligned or isotropic. For simplicity, we assume both black holes in the system come from the same distribution (generator), but you can also specify different distributions for the primary and secondary black holes if desired.

Note that we use mass ratio in [1, 6] because surfinbh's models restrict to this range.

In [4]:

Copied!





binary_generator = BinaryGenerator(
    primary_black_hole_source=bh_generator,
    secondary_black_hole_source=bh_generator,
    mass_ratio=Domain(low=1.0, high=6.0),
    is_aligned_spin=True,
)
binary_generator = BinaryGenerator(
    primary_black_hole_source=bh_generator,
    secondary_black_hole_source=bh_generator,
    mass_ratio=Domain(low=1.0, high=6.0),
    is_aligned_spin=True,
)

Then we run the simulation using the simulate_black_hole_mergers function, which takes in the binary generator and other parameters such as the number of binaries. n_workers specifies the number of parallel workers to use for the simulation, which can speed up the computation significantly when simulating a large number of binaries.

In [5]:

Copied!





size = 5000
n_workers = 1
random_state = 42
fits = Fits.NRSUR3DQ8REMNANT
# If you want to simulate a population in precession setting,
# please use Fits.NRSUR7DQ4REMNANT instead.

black_hole_mergers = simulate_black_hole_mergers(
    binary_generator, fits, size, n_workers, random_state
)
size = 5000
n_workers = 1
random_state = 42
fits = Fits.NRSUR3DQ8REMNANT
# If you want to simulate a population in precession setting,
# please use Fits.NRSUR7DQ4REMNANT instead.

black_hole_mergers = simulate_black_hole_mergers(
    binary_generator, fits, size, n_workers, random_state
)

Loaded NRSur3dq8Remnant fit.

100%|██████████| 5000/5000 [00:00<00:00, 1783899.29it/s]

Here we show how a black hole merger sample actually looks like. Each entry is a tuple of Binary and BlackHole (remnant) data structures, which contain all the parameters of the binary and the remnant black hole, respectively. You can access any parameter of interest from these data structures for further analysis or visualization.

In [6]:

Copied!

binary, remnant = black_hole_mergers[0]
binary, remnant = black_hole_mergers[0]

In [7]:

Copied!

binary
binary

Out[7]:

Binary(primary_black_hole=BlackHole(mass=17.893000000004307, spin_magnitude=0.09237083927589662, spin_vector=(0.0, 0.0, -0.09237083927589662), speed=0.0), secondary_black_hole=BlackHole(mass=10.223000000001743, spin_magnitude=0.3936355202774421, spin_vector=(0.0, 0.0, -0.3936355202774421), speed=0.0))

In [8]:

Copied!

remnant
remnant

Out[8]:

BlackHole(mass=27.047814394642238, spin_magnitude=0.5891627875731139, spin_vector=(0.0, 0.0, 0.5891627875731139), speed=111.98845275921877)

Finally we can also convert the simulation results into a pandas data fram for easier manipulation and analysis, which is especially useful when dealing with large datasets, and visualizing the dataset with packages like matplotlib and seaborn.

In [9]:

Copied!

df_binaries = convert_simulated_binaries_to_dataframe(black_hole_mergers)

df_binaries.head(10)
df_binaries = convert_simulated_binaries_to_dataframe(black_hole_mergers)

df_binaries.head(10)

Out[9]:

	m_1	a_1	a_1z	m_2	a_2	a_2z	m_f	a_f	k_f	chi_eff	q
0	17.893	0.092371	-0.092371	10.223	0.393636	-0.393636	27.047814	0.589163	111.988453	-0.201911	1.750269
1	41.465	0.473436	-0.473436	8.446	0.060754	0.060754	49.205834	0.140479	176.289219	-0.383039	4.909425
2	17.976	0.854547	-0.854547	9.915	0.604192	0.604192	26.985044	0.425651	407.402142	-0.335979	1.813011
3	15.542	0.966116	0.966116	13.847	0.340004	0.340004	27.063866	0.881890	75.158232	0.671116	1.122409
4	11.652	0.502721	-0.502721	8.207	0.869650	-0.869650	19.191158	0.460602	94.414534	-0.654359	1.419764
5	18.257	0.051515	0.051515	8.207	0.088134	-0.088134	25.491424	0.617676	125.350637	0.008207	2.224564
6	26.853	0.847548	-0.847548	7.861	0.244793	-0.244793	34.103435	0.084981	278.162089	-0.711054	3.415978
7	17.090	0.337857	0.337857	13.903	0.181818	0.181818	29.254119	0.765575	26.198156	0.267860	1.229231
8	17.034	0.430347	0.430347	8.826	0.043603	0.043603	24.572559	0.761794	71.287979	0.298351	1.929980
9	10.704	0.904893	-0.904893	6.706	0.165500	0.165500	16.835611	0.429193	357.283266	-0.492598	1.596183