Your First Binary Simulations with POSYDON 🌠

Tutorial goal:

• Create your first binary population

New concepts:

• Default population synthesis model

• Synthetic Population class

If you haven’t done so yet, export the path POSYDON environment variables. Set these parameters in your .bash_profile or .zshrc if you use POSYDON regularly.

%env PATH_TO_POSYDON=/YOUR/POSYDON/PATH/
%env PATH_TO_POSYDON_DATA=/YOUR/POSYDON_DATA/PATH/

env: PATH_TO_POSYDON=/Users/simone/Google Drive/github/POSYDON-public/
env: PATH_TO_POSYDON_DATA=/Volumes/T7/

Creating the Initialisation File

---

To run population synthesis with POSYDON, a population_params.ini file is required. This file described how the stellar population is created and what prescriptions are implemented.

POSYDON comes with a default population_params_detault.ini file found at PATH_TO_POSYDON/posydon/popsyn or have a look here.

Here we will run the notebook with the default parameters, which runs 10 binaries at \(10^{-4}Z_\odot\).

First, let’s copy the default population synthesis model to your working directory.

import os
import shutil
from posydon.config import PATH_TO_POSYDON

path_to_params = os.path.join(PATH_TO_POSYDON, "posydon/popsyn/population_params_default.ini")
shutil.copyfile(path_to_params, './population_params.ini')

'./population_params.ini'

Running the Population Synthesis Model

SyntheticPopulation() samples the given initial distributions given in the population_params.ini file and generates the initial conditions for our 10 binaries.

import pandas as pd
from posydon.popsyn.synthetic_population import SyntheticPopulation

synth_pop = SyntheticPopulation('./population_params.ini')
synth_pop.evolve()

Z=1.00e-04 Z_sun

/Users/simone/Google Drive/github/POSYDON-public/posydon/popsyn/binarypopulation.py:571: FutureWarning: The behavior of DataFrame concatenation with empty or all-NA entries is deprecated. In a future version, this will no longer exclude empty or all-NA columns when determining the result dtypes. To retain the old behavior, exclude the relevant entries before the concat operation.
  return pd.concat(holder, axis=0, ignore_index=False)

In the previous step, a synthetic binary population is created with the POSYDON default initial parameters and stored in synth_pop. With the last line of code the synthetic population will be evolved.

Congratulations! You run your first population synthesis model with POSYDON! 🎉

Exploring the Simulation Output

Here are the population synthesis files you just created, in this example we only visualize the first ten lines of the output table:

df = pd.read_hdf('./1.00e-04_Zsun_population.h5', key='history')
df.head(10)

	state	event	time	orbital_period	eccentricity	lg_mtransfer_rate	step_names	step_times	S1_state	S1_mass	...	S2_he_core_mass	S2_he_core_radius	S2_co_core_mass	S2_co_core_radius	S2_center_h1	S2_center_he4	S2_surface_h1	S2_surface_he4	S2_surf_avg_omega_div_omega_crit	S2_spin
binary_index
0	detached	ZAMS	1.511066e+09	479.288083	0.0	NaN	initial_cond	0.000000	H-rich_Core_H_burning	10.635159	...	NaN	NaN	NaN	NaN	0.750999	2.490000e-01	NaN	NaN	NaN	NaN
0	RLO1	CC1	1.535382e+09	1220.670243	0.0	-4.442681	step_HMS_HMS	2.176263	H-rich_Central_C_depletion	4.387056	...	1.411542e-14	2.053381e-16	1.365520e-14	1.239204e-16	0.423320	5.766784e-01	0.388622	0.611377	0.907345	15.586620
0	disrupted	NaN	1.535382e+09	NaN	NaN	NaN	step_SN	0.001554	NS	1.346469	...	1.411542e-14	2.053381e-16	1.365520e-14	1.239204e-16	0.423320	5.766784e-01	0.388622	0.611377	0.907345	15.586620
0	disrupted	CC2	1.548703e+09	NaN	NaN	NaN	step_disrupted	5.963153	NS	1.346469	...	2.970098e+00	NaN	2.062403e+00	5.067320e-02	0.000000	7.232097e-14	0.719862	0.280118	0.024038	13.729559
0	disrupted	NaN	1.548703e+09	NaN	NaN	NaN	step_SN	0.001684	NS	1.346469	...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	0.000000
0	disrupted	END	1.548703e+09	NaN	NaN	NaN	step_end	0.000107	NS	1.346469	...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	0.000000
1	detached	ZAMS	3.316347e+09	0.911049	0.0	NaN	initial_cond	0.000000	H-rich_Core_H_burning	15.575811	...	NaN	NaN	NaN	NaN	0.750999	2.490000e-01	NaN	NaN	NaN	NaN
1	contact	oCE1	3.327664e+09	0.506718	0.0	-1.876899	step_HMS_HMS	0.347768	H-rich_Core_H_burning	14.908113	...	1.411542e-14	2.053381e-16	1.365520e-14	1.239204e-16	0.722955	2.770430e-01	0.750621	0.249361	0.988417	10.460393
1	merged	oMerging1	3.327664e+09	0.506718	0.0	-1.876899	step_CE	0.000091	H-rich_Core_H_burning	14.908113	...	1.411542e-14	2.053381e-16	1.365520e-14	1.239204e-16	0.722955	2.770430e-01	0.750621	0.249361	0.988417	10.460393
1	merged	CC1	3.331178e+09	NaN	NaN	NaN	step_merged	6.616714	H-rich_Central_C_depletion	16.152142	...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN

10 rows × 38 columns

It is also possible to collapse the evolution output on a single row per binary system as shown in the next example.

df_oneline = pd.read_hdf('./1.00e-04_Zsun_population.h5', key='oneline')
df_oneline.head(10)

	state_i	event_i	time_i	orbital_period_i	eccentricity_i	lg_mtransfer_rate_i	step_names_i	step_times_i	state_f	event_f	...	interp_class_HMS_HMS	interp_class_CO_HMS_RLO	interp_class_CO_HeMS	interp_class_CO_HeMS_RLO	mt_history_HMS_HMS	mt_history_CO_HMS_RLO	mt_history_CO_HeMS	mt_history_CO_HeMS_RLO	FAILED	WARNING
binary_index
0	detached	ZAMS	1.511066e+09	479.288083	0.0	NaN	initial_cond	0.0	disrupted	END	...	stable_MT	NaN	NaN	NaN	Stable RLOF during postMS	NaN	NaN	NaN	0	1
1	detached	ZAMS	3.316347e+09	0.911049	0.0	NaN	initial_cond	0.0	merged	END	...	unstable_MT	NaN	NaN	NaN	Unstable contact phase	NaN	NaN	NaN	0	1
2	detached	ZAMS	7.997196e+08	0.964342	0.0	NaN	initial_cond	0.0	merged	END	...	unstable_MT	NaN	NaN	NaN	Unstable contact phase	NaN	NaN	NaN	0	1
3	detached	ZAMS	6.368205e+08	1.807444	0.0	NaN	initial_cond	0.0	disrupted	END	...	stable_MT	NaN	NaN	NaN	Stable RLOF during postMS	NaN	NaN	NaN	0	1
4	detached	ZAMS	5.028601e+09	1.688541	0.0	NaN	initial_cond	0.0	initial_RLOF	END	...	stable_MT	NaN	NaN	NaN	Stable RLOF during postMS	NaN	NaN	NaN	0	0
5	detached	ZAMS	5.384771e+09	59.604185	0.0	NaN	initial_cond	0.0	disrupted	END	...	stable_MT	NaN	NaN	NaN	Stable RLOF during postMS	NaN	NaN	NaN	0	1
6	detached	ZAMS	1.242091e+10	3622.599836	0.0	NaN	initial_cond	0.0	disrupted	END	...	no_MT	NaN	NaN	NaN	no RLOF	NaN	NaN	NaN	0	1
7	detached	ZAMS	1.278701e+10	2.785114	0.0	NaN	initial_cond	0.0	merged	END	...	unstable_MT	NaN	NaN	NaN	Unstable RLOF during postMS	NaN	NaN	NaN	0	1
8	detached	ZAMS	3.362532e+09	2.522910	0.0	NaN	initial_cond	0.0	detached	END	...	stable_MT	stable_MT	NaN	NaN	Stable contact phase	Stable RLOF during postMS	NaN	NaN	0	1
9	detached	ZAMS	7.402367e+09	4518.548242	0.0	NaN	initial_cond	0.0	disrupted	END	...	no_MT	NaN	NaN	NaN	no RLOF	NaN	NaN	NaN	0	1

10 rows × 100 columns

For the last example, only a selection of columns will be selected to visualize the evolution of the ninth binary (with a corresponding index value of eight) in the population:

# let's check the BBH system at index 8
cols = ['time', 'step_names', 'state', 'event', 'orbital_period', 'eccentricity', 'S1_state', 'S2_state', 'S1_mass', 'S2_mass']
df.loc[8, cols]

	time	step_names	state	event	orbital_period	eccentricity	S1_state	S2_state	S1_mass	S2_mass
binary_index
8	3.362532e+09	initial_cond	detached	ZAMS	2.522910	0.000000	H-rich_Core_H_burning	H-rich_Core_H_burning	93.546567	93.023072
8	3.366016e+09	step_HMS_HMS	detached	CC1	3.909472	0.000000	H-rich_Central_C_depletion	H-rich_Core_H_burning	58.365944	119.598219
8	3.366016e+09	step_SN	detached	NaN	4.648114	0.083735	BH	H-rich_Core_H_burning	43.861220	119.598219
8	3.366233e+09	step_detached	RLO2	oRLO2	4.902223	0.000000	BH	H-rich_Core_H_burning	43.861220	101.734326
8	3.366846e+09	step_CO_HMS_RLO	detached	CC2	2.846054	0.000000	BH	H-rich_Central_C_depletion	45.175531	67.246779
8	3.366846e+09	step_SN	detached	NaN	8.712412	0.458613	BH	BH	45.175531	30.278034
8	1.380000e+10	step_dco	detached	maxtime	5.781819	0.355638	BH	BH	45.175531	30.278034
8	1.380000e+10	step_end	detached	END	5.781819	0.355638	BH	BH	45.175531	30.278034

Feel free to explore the dataset!

If you want to learn more about population synthesis and how to build more complex models it is advised to continue with the remaining tutorials and consult the POSYDON documentation.

Local MPI runs

To speed up population synthesis runs, you can run on a computing cluster, as described in HPC Facilities, or you can distribute the population synthesis across multiple cores on your local machine using MPI.

To enable local MPI runs, go into the population_params_detault.ini and change use_MPI to True.

It’s important to note that you cannot run have this option enabled for cluster runs!

We create a binary population simulation script to run the population:

%%writefile script.py
from posydon.popsyn.synthetic_population import SyntheticPopulation

if __name__ == "__main__":
    synth_pop = SyntheticPopulation("./population_params.ini")
    synth_pop.evolve()

This script can be initiated using a local where NR_processors is the number of processors you would like to us.

mpiexec -n ${NR_processors} python script.py

This will create a folder for each metallicity in the population and store output of the parallel runs in it. You will have to concatenate these runs into a single population file per metallicity, which can be achieved using the following code:

from posydon.popsyn.synthetic_population import SyntheticPopulation
from posydon.config import PATH_TO_POSYDON_DATA
import os


synth_pop = SyntheticPopulation("./population_params.ini")
# Get the path to the batches in the current folder
x = os.listdir('.')
path_to_batches = [i for i in x if i.endswith('_batches')]
synth_pop.merge_parallel_runs(path_to_batches)