Stellar Transient Populations 🔍

In this notebook, we will cover:

Selecting binaries and combining metallicities
Exploring the TransientPopulation class
Calculating cosmic rates

We will do this in the context of the merging binary black hole population.

Creating multi-metallicity Populations

Reprocessed POSYDON v1 dataset

Please note that with the reprocessed POSYDON v1 data, only solar metallicity is available. You will not be able to follow along with the full tutorial! If you would still like to explore a population at solar metallicity, you can follow the One metallicity tutorial.

In the previous tutorial, you generated 8 population files with 1000 binaries. Since this is a small population of only 8.000 binaries, you could explore the complete population. But the larger your populations get, the better it is to select specific binaries.

For each population file, we can export the binaries we’re interested in into a new file. For this, we need to find the indices of the merging binaries. The relevant properties are that both S1_state and S2_state equal ‘BH’, while the binary has the event == 'CO_contact.

We will load one of the population files to build our merging binaries selection.

No indices?

It might be that the population you generated does not contain any merging BBHs, they’re rare after all.

As such, we have provided an example population run at $PATH_TO_POSYDON_DATA/population-synthesis/example/, this simulation contains 10.000 binaries at the 8 metallicity.

[ ]:

from posydon.popsyn.synthetic_population import Population

pop = Population('1e+00_Zsun_population.h5')
tmp_data = pop.history.select(columns=['S1_state', 'S2_state', 'event'])

pop.history.select is a read operation on the population file and can take quite a lot of time if the population is large.

If you have sufficient memory, it is more efficient to select several columns at the same time, as done here, instead of selecting a single column each time.

[ ]:

# Selection of S1 being a BH
S1_state = tmp_data['S1_state'] == 'BH'
# Selection of S2 being a BH
S2_state = tmp_data['S2_state'] == 'BH'
# Selection of the binary system being in contact during the double CO phase.
state = tmp_data['event'] == 'CO_contact'

# get the indices of all systems
indices = tmp_data.index

# delete the temporary data
del tmp_data

# get a mask for the indices that satisfy all the conditions
mask = S1_state & S2_state & state

# get the indices that satisfy all the conditions
selected_indices = indices[mask].to_list()

The selected_indices has to be a list for the selection to work correctly.

You can test you selection by doing the following:

[ ]:

print(selected_indices)
pop.history[selected_indices]

If your selected indices are empty, no BBH mergers occurred in your population. They’re quite rare after all!

Let’s use an example set of populations to make sure there will be BBH in your population.

This can only be used with the multi-metallicity populations!

[ ]:

from posydon.popsyn.synthetic_population import Population

data_path = f'PATH/TO/TUTORIAL/DATA/'

# Let's load one of the populations
pop = Population(data_path+'1e-01_Zsun_population.h5')

[ ]:

pop.oneline.columns

We don’t just want to write these binaries for a single metallicity. We want to export them for each metallicity file. As such, we loop over each Population file and export the binaries to the same file. It’s important to have the same columns in the populations, otherwise it’s not possible to add them together.

The indices of the binaries will be reset, when being exported, such that the new file contains unique indices for each binary.

[ ]:

from posydon.popsyn.synthetic_population import Population

# The names of the populations, if you followed along.
files = ['1e-01_Zsun_population.h5',
         '1e-02_Zsun_population.h5',
         '1e-03_Zsun_population.h5',
         '1e-04_Zsun_population.h5',
         '1e+00_Zsun_population.h5',
         '2e-01_Zsun_population.h5',
         '2e+00_Zsun_population.h5',
         '4.5e-01_Zsun_population.h5']


for file in files:
    pop = Population(data_path+file)
    # read the relevant data in one go
    # (faster than reading it in chunks, but requires more memory)
    tmp_data = pop.history.select(columns=['S1_state', 'S2_state', 'event'])
    # Selection of S1 being a BH
    S1_state = tmp_data['S1_state'] == 'BH'
    # Selection of S2 being a BH
    S2_state = tmp_data['S2_state'] == 'BH'
    # Selection of the binary system being in contact during the double CO phase.
    state = tmp_data['event'] == 'CO_contact'
    indices = tmp_data.index
    del tmp_data
    mask = S1_state & S2_state & state
    selected_indices = indices[mask].to_list()
    print(f'File: {file}, Number of systems: {len(selected_indices)}')

    # set overwrite to False to add to the file
    pop.export_selection(selected_indices, data_path+'BBH_contact.h5', append=True)

If your population did not contain any BBH mergers, all the “Number of Systems” will be 0, and you won’t be able to open the file in the next cell.

If you ran the population with the tutorial populations, you should have a total of 636 merging BBHs.

We can confirm this by adding up the values above or by opening the new file. Population.number_of_systems gives us the total number of binaries in the file.

You now see that the mass_per_metallicity property contains information about all the metallicities in the file.

You can even combine multiple runs at the same metallicity together, if you like. This will combine their simulated mass.

[ ]:

from posydon.popsyn.synthetic_population import Population
BBH_pop = Population(data_path+'BBH_contact.h5', chunksize=10000)
print(BBH_pop.number_of_systems)

BBH_pop.mass_per_metallicity

Underlying mass

Before we move on, we have to discuss the occurrence rate of each binary in the population.

We have sampled several distributions at the start of our simulations. These include:

Initial mass function
Period
mass ratio
binary fraction

The default parameters do not sample the complete distributions. For example, we sample the binary between 7 and 150 $M_\odot$, and a binary fraction of 1. While you can use the population as is, you should consider the weight of the distributions not sampled.

You can calculate the underlying mass of these distributions, which is the actual weight of the binary in our population. Since we’ve stored essential information about the population run, we’re able to calculate this even after our selection. This will return the underlying_mass of the sampled population and add it to the population file.

We will need these weights to be able to continue with the rest of the calculations.

Caveats

Currently, only changing the binary fraction (f_bin) is supported in the calculation of the underlying mass! The initially sampled population can have any binary fraction.

Additionally, only parts of the period, IMF, and mass ratio distribution are included. We allow for different parameter limits to be set for sampling, but these are not included in this normalisation!

[ ]:

# This will calculate the underlying mass for the population
# with a binary fraction of 0.7 (default value)
BBH_pop.calculate_underlying_mass(f_bin=0.7, overwrite=True)

Transient population

Although we now have selected all binaries with a BBH merger, we don’t have the exact moment of the merger yet. This will be required to calculate merger rates across cosmic time.

For this we will create a TransientPopulation. This class is used to hold information about a specific event/moment in time. In our case, this is the moment of “CO_contact”, the moment of the BBH merger. However, we might want to store and calculate some additional values, such as the $M_\mathrm{chirp}$ or $\chi_\mathrm{eff}$.

The Population class has a function create_transient_population (see here for more details). In short, it takes a selection_function and a transient_name. The transient_name is a string identifying the transient population in the file, while selection_function extracts the TransientPopulation for us. This can be any custom function you want it to be, as long it outputs a pandas DataFrame with a ‘time’ and ‘metallicity’ column.

Several selection functions are provided in posydon.popsyn.transient_select_funcs:

BBH_selection_function
GRB_selection_function

We have copied part of the BBH_selection_function below to explain how a selections_function works.

[ ]:

import pandas as pd
def BBH_selection_function(history_chunk, oneline_chunk, formation_channels_chunk=None):
    '''A BBH selection function to create a transient population of BBHs mergers.'''

    indices = oneline_chunk.index.to_numpy()
    df_transients = pd.DataFrame(index = indices)

    df_transients['time'] = history_chunk[history_chunk['event'] == 'CO_contact']['time'] * 1e-6 #Myr
    df_transients['metallicity'] = oneline_chunk['metallicity']

    return df_transients

A selection_function always has as an input a chunk of the history, oneline, and formation_pathways (optional). For example, you set your chunksize=10000 when you initialise the Population, like we’ve done above for BBH_pop, each chunk will contain 10000 binaries. This means that the complete history, oneline and formation_pathways of those 10 binaries are passed to this function.

The BBH_selection_function selects the moment the binary reaches CO_contact as the moment of merger, and stores it in the time columns in Myr. The metallicity is also outputted, since this will be essential when combining it with the star formation history.

We can test this function by inputting a single binary into it, as done below. If you’ve used the example populations, you might not have calculated the formation_channels yet. We will set that input to None.

[ ]:

BBH_selection_function(BBH_pop.history[0], BBH_pop.oneline[0], None)

This gives us a dataframe containing our 1 binary, its index (0), the moment of merger, and its metallicity.

Of course, we could like to know a bit more about our merger than just when it occurs. Let’s expand our output with the BH masses, their spin, their tilt, and the orbital period at DCO formation.

[ ]:

def BBH_selection_function(history_chunk, oneline_chunk, formation_channels_chunk):
    '''A BBH selection function to create a transient population of BBHs mergers.'''

    indices = oneline_chunk.index.to_numpy()
    df_transients = pd.DataFrame(index = indices)

    df_transients['time'] = history_chunk[history_chunk['event'] == 'CO_contact']['time'] * 1e-6 #Myr
    df_transients['metallicity'] = oneline_chunk['metallicity']

    # Added properties
    mask = (history_chunk['S1_state'] == 'BH') & (history_chunk['S2_state'] == 'BH') & (history_chunk['step_names'] == 'step_SN') & (history_chunk['state'] == 'detached')
    df_transients['t_inspiral'] = df_transients['time'] - history_chunk[mask]['time']*1e-6
    df_transients['S1_state']  = history_chunk[mask]['S1_state']
    df_transients['S2_state']  = history_chunk[mask]['S2_state']
    df_transients['S1_mass'] = history_chunk[mask]['S1_mass']
    df_transients['S2_mass'] = history_chunk[mask]['S2_mass']
    df_transients['S1_spin'] = history_chunk[mask]['S1_spin']
    df_transients['S2_spin'] = history_chunk[mask]['S2_spin']
    # we distinguish the tilt of the spin to the orbit after the first and second SN.
    df_transients['S1_spin_orbit_tilt_merger'] = oneline_chunk['S1_spin_orbit_tilt_second_SN']
    df_transients['S2_spin_orbit_tilt_merger'] = oneline_chunk['S2_spin_orbit_tilt_second_SN']
    df_transients['orbital_period'] = history_chunk[mask]['orbital_period']

    return df_transients

[ ]:

BBH_selection_function(BBH_pop.history[0], BBH_pop.oneline[0], None)

With this new function, we have a lot more information available in the TransientPopulation.

You can further customise this to your liking, if you want to store specific information. It’s also possible to calculate additional information based on any value in the history, oneline or formation_channels.

We will import some functions to calculate $\chi_\mathrm{eff}$, $q$ and $\mathcal{M}_\mathrm{chirp}$ and also add them to the dataframe.

[ ]:

from posydon.popsyn.transient_select_funcs import chi_eff, mass_ratio, m_chirp

def BBH_selection_function(history_chunk, oneline_chunk, formation_channels_chunk):
    '''A BBH selection function to create a transient population of BBHs mergers.'''

    indices = oneline_chunk.index.to_numpy()
    df_transients = pd.DataFrame(index = indices)

    df_transients['time'] = history_chunk[history_chunk['event'] == 'CO_contact']['time'] * 1e-6 #Myr
    mask = (history_chunk['S1_state'] == 'BH') & (history_chunk['S2_state'] == 'BH') & (history_chunk['step_names'] == 'step_SN') & (history_chunk['state'] == 'detached')
    df_transients['metallicity'] = oneline_chunk['metallicity']
    df_transients['t_inspiral'] = df_transients['time'] - history_chunk[mask]['time']*1e-6

    df_transients['S1_state']  = history_chunk[mask]['S1_state']
    df_transients['S2_state']  = history_chunk[mask]['S2_state']
    df_transients['S1_mass'] = history_chunk[mask]['S1_mass']
    df_transients['S2_mass'] = history_chunk[mask]['S2_mass']
    df_transients['S1_spin'] = history_chunk[mask]['S1_spin']
    df_transients['S2_spin'] = history_chunk[mask]['S2_spin']
    # we distinguish the tilt of the spin to the orbit after the first and second SN.
    df_transients['S1_spin_orbit_tilt_merger'] = oneline_chunk['S1_spin_orbit_tilt_second_SN']
    df_transients['S2_spin_orbit_tilt_merger'] = oneline_chunk['S2_spin_orbit_tilt_second_SN']
    df_transients['orbital_period'] = history_chunk[mask]['orbital_period']
    df_transients['eccentricity'] = history_chunk[mask]['eccentricity']

    # Added
    df_transients['chirp_mass'] = m_chirp(history_chunk[mask]['S1_mass'], history_chunk[mask]['S2_mass'])
    df_transients['mass_ratio'] = mass_ratio(history_chunk[mask]['S1_mass'], history_chunk[mask]['S2_mass'])
    df_transients['chi_eff'] = chi_eff(history_chunk[mask]['S1_mass'],
                                       history_chunk[mask]['S2_mass'],
                                       history_chunk[mask]['S1_spin'],
                                       history_chunk[mask]['S2_spin'],
                                       oneline_chunk['S1_spin_orbit_tilt_second_SN'],
                                       oneline_chunk['S2_spin_orbit_tilt_second_SN'])

    return df_transients

[ ]:

BBH_selection_function(BBH_pop.history[0], BBH_pop.oneline[0], None)

Before running this selection on the whole population, we would like to include the formation_channels too. These are calculated for the systems in the population file and add a good overview of the evolution of binaries.

[ ]:

# mt_history=True adds detailed information about the mass transfer in the HMS-HMS grid
BBH_pop.calculate_formation_channels(mt_history=True)

Now we can include the formation_channels in df_transients.

[ ]:

def BBH_selection_function(history_chunk, oneline_chunk, formation_channels_chunk):
    '''A BBH selection function to create a transient population of BBHs mergers.'''

    indices = oneline_chunk.index.to_numpy()
    df_transients = pd.DataFrame(index = indices)

    df_transients['time'] = history_chunk[history_chunk['event'] == 'CO_contact']['time'] * 1e-6 #Myr
    mask = (history_chunk['S1_state'] == 'BH') & (history_chunk['S2_state'] == 'BH') & (history_chunk['step_names'] == 'step_SN') & (history_chunk['state'] == 'detached')
    df_transients['metallicity'] = oneline_chunk['metallicity']
    df_transients['t_inspiral'] = df_transients['time'] - history_chunk[mask]['time']*1e-6

    df_transients['S1_state']  = history_chunk[mask]['S1_state']
    df_transients['S2_state']  = history_chunk[mask]['S2_state']
    df_transients['S1_mass'] = history_chunk[mask]['S1_mass']
    df_transients['S2_mass'] = history_chunk[mask]['S2_mass']
    df_transients['S1_spin'] = history_chunk[mask]['S1_spin']
    df_transients['S2_spin'] = history_chunk[mask]['S2_spin']
    # we distinguish the tilt of the spin to the orbit after the first and second SN.
    df_transients['S1_spin_orbit_tilt_merger'] = oneline_chunk['S1_spin_orbit_tilt_second_SN']
    df_transients['S2_spin_orbit_tilt_merger'] = oneline_chunk['S2_spin_orbit_tilt_second_SN']
    df_transients['orbital_period'] = history_chunk[mask]['orbital_period']
    df_transients['eccentricity'] = history_chunk[mask]['eccentricity']

    df_transients['chirp_mass'] = m_chirp(history_chunk[mask]['S1_mass'], history_chunk[mask]['S2_mass'])
    df_transients['mass_ratio'] = mass_ratio(history_chunk[mask]['S1_mass'], history_chunk[mask]['S2_mass'])
    df_transients['chi_eff'] = chi_eff(history_chunk[mask]['S1_mass'],
                                       history_chunk[mask]['S2_mass'],
                                       history_chunk[mask]['S1_spin'],
                                       history_chunk[mask]['S2_spin'],
                                       oneline_chunk['S1_spin_orbit_tilt_second_SN'],
                                       oneline_chunk['S2_spin_orbit_tilt_second_SN'])

    # added
    df_transients = pd.concat([df_transients, formation_channels_chunk[['channel']]], axis=1)

    return df_transients

[ ]:

# NOTE: the formation channels have to be given as a dataframe, as such we call it with loc[[0]] for testing
# Make sure you use double brackets for loc to return a dataframe
BBH_selection_function(BBH_pop.history[0],
                       BBH_pop.oneline[0],
                       BBH_pop.formation_channels.loc[[0]])

The example binary gives us a dataframe with all the information we’re interested in. But if something is missing feel free to customize the selection further!

We now create the transient population with the BBH_selection_function, we’ve created and create_transient_population.

This will create a TransientPopulation instance and write the transient population to the current file.

BBH_mergers.population will load the complete dataframe we’ve just created. Each of the indices in this population refer back to the original binaries, which are still accessible through BBH_mergers.history or BBH_mergers.oneline.

[ ]:

BBH_mergers = BBH_pop.create_transient_population(BBH_selection_function, 'BBH')

[ ]:

BBH_mergers.population

Since the Transient Population is stored in the same file, you can continue your analysis by opening the file using the TransientPopulation class. You only need to do your selection once and you can then continue from this stage.

[ ]:

from posydon.popsyn.synthetic_population import TransientPopulation

BBH_mergers = TransientPopulation(filename='BBH_contact.h5', transient_name='BBH')

[ ]:

BBH_mergers.mass_per_metallicity

With TransientPopulation class you can compute the efficiency of the transient population at a given metallicity, i.e. the number of merging DCO per unit star formation and visualize the results.

Underlying mass

Make sure you’ve calculated the underlying mass for the population, as done earlier in this tutorial.

Applying the underlying mass normalisation assumes that no events occur outside the sampled parameter space.

[ ]:

BBH_mergers.get_efficiency_over_metallicity()

[ ]:

# if you want to plot the channels, you can set channels=True
BBH_mergers.plot_efficiency_over_metallicity(channels=True)

[ ]:

BBH_mergers.efficiency

You can also plot the delay time distribution based on the time columns in the TransientPopulation.

[ ]:

import numpy as np
bins = np.logspace(6.5,10.5)
BBH_mergers.plot_delay_time_distribution(bins=bins)
# you can also set a specific metallicity. However, this normalises the rate with only the mass of the selected metallicity.
BBH_mergers.plot_delay_time_distribution(metallicity=0.1, bins=bins)

You might also be interested in seeing the evolution of the binary in more detail.

Therefore, you can plot the TransientPopulation over a grid slice. For example, below we plot binaries over the $q=0.7$ slice at $Z=10^{-4} Z_\odot$ on the HMS-HMS grid.

If no slice is give, all mass ratios are plotted.

You can also plot a property from the TransientPopulation.population DataFrame as a colourmap. And specific formation_channels if they’re in the DataFrame.

[ ]:

BBH_mergers.plot_popsyn_over_grid_slice('HMS-HMS', 1e-4, slices=[0.7], save_fig=False)

[ ]:

BBH_mergers.plot_popsyn_over_grid_slice('HMS-HMS', 1e-4, slices=[0.7], prop='S1_spin', prop_range=[0,0.3], save_fig=False, channel='ZAMS_oRLO1_CC1_oRLO2_CC2_END') # SMT channel

Cosmic Rate

So far these events do not consider the metallicity or star formation rate evolution of the Universe.

POSYDON comes with several built-in star formation histories and metallicity evolutions, a few examples are:

IllustrisTNG
Neijssel2019
Madau+Fragos2017

We can apply these to our population with the calculate_cosmic_weights function. Here, we will the metallicity and SFR evolution of the IllustrisTNG, which is the default model used, if no MODEL_in is given.

The function returns an instance of the Rates class, which gives us access to some new variables:

z_birth: the redshift and age of the universe at which we probe the star formation
z_events: the redshift at which an event takes place
weights: the weight of the event based on the SFR, its metallicity, and its weight in the population.

[ ]:

from posydon.popsyn.synthetic_population import TransientPopulation
BBH_mergers = TransientPopulation(filename='BBH_contact.h5', transient_name='BBH')

[ ]:

MODEL = {
    'delta_t' : 100, # Myr
    'SFR' : 'IllustrisTNG', # Neijssel2019, Madau+Fragos2017
    'sigma_SFR' : None,
    'Z_max' : 2.,
}

rates = BBH_mergers.calculate_cosmic_weights('IllustrisTNG', MODEL_in=MODEL)

[ ]:

rates.z_birth.head(2)

[ ]:

rates.z_events.head(2)

[ ]:

rates.weights.head(2)

Similar to the TransientPopulation, the rates can also be accessed using the Rates class. You will need to provide your transient_name and SFH_identifier. These rates are also stored in the file. You can have as many transients and star formation histories in your population file as you want.

[ ]:

from posydon.popsyn.synthetic_population import Rates

rates = Rates(filename='BBH_contact.h5', transient_name='BBH', SFH_identifier='IllustrisTNG')

You can directly work with the weights and the z_events to get events at specific redshifts.

Or you can calculate the rate density over redshift using the calculate_intrinsic_rate_density function. This will calculate the rate over redshift for you.

You can use the output of the function immediately, or get them from rates.intrinsic_rate_density.

[ ]:

out = rates.calculate_intrinsic_rate_density(mt_channels=True)

[ ]:

rates.intrinsic_rate_density

[ ]:

rates.plot_intrinsic_rate(channels=True, xlim=(0,10), ylim=(1e-2,1e3))

Properties of the systems

Sometime syou mgiht want some more details about the actual population, and its properties.

plot_hist_properties() allows you to plot these.

By default, the function will create its own figure and plot it. If you would like more control over the output, you can give it your own pyplot axis and set show=False. This allows you to adapt and change the plot to your liking after finishing adding the properties.

[ ]:

import numpy as np
import matplotlib.pyplot as plt

bins = np.linspace(0,100,101)
fig, ax = plt.subplots(1,1)
rates.plot_hist_properties('S1_mass', intrinsice=True, bins=bins, color='red', ax =ax, label='S1', show=False)
rates.plot_hist_properties('S2_mass', intrinsice=True, bins=bins, ax=ax, label='S2', show=False)

ax.set_ylabel('Rate density [Gpc$^-3$ yr$^-1$]')
ax.set_xlabel('Mass [Msun]')
ax.legend()
plt.show()

Observable population

Sometimes, however, you’re not just interested in the intrinsic population, and want an observable population. This could be a supernova detection fraction for a telescope survey, or the LVK detection efficiency for GW mergers.

You can apply these using the calculate_observable_population. Similar to the create_transient_population, this function takes a observable_func, which described the observability of a transient.

An observable_func takes chunks of

TransientPopulation.population
The TransientPopulation.z_events
The TransientPopulation.weights

Using these, new weights are calculated and stored as a separate observable population of your TransientPopulation. The BBH analysis framework comes with a detection function for several detector sensitivities and configurations.

[ ]:

from posydon.popsyn.transient_select_funcs import DCO_detactability

[ ]:

def DCO_wrapper(transient_chunk, z_events_chunk, weights_chunk):
    sensitivity = 'design_H1L1V1'
    return DCO_detactability(sensitivity, transient_chunk, z_events_chunk, weights_chunk, verbose=False)

[ ]:

# We also give it a name, which is used as an identifier in the file
rates.calculate_observable_population(DCO_wrapper, 'design_H1L1V1')

[ ]:

# We can now access this observable population
rates.observable_population('design_H1L1V1')

You can also plot the intrinsic and observable population together.

[ ]:

bins = np.linspace(0,100,101)
fig, ax = plt.subplots(1,1)

rates.plot_hist_properties('S1_mass', intrinsice=True, observable='design_H1L1V1', bins=bins, ax = ax, label='S1', show=False)
ax.set_ylabel('Rate density [Gpc$^-3$ yr$^-1$]')
ax.set_xlabel('Mass [Msun]')
ax.legend()
plt.show()

Cogratulations, you are now ready to analyze any DCO population data you generated with POSYDON. Feel free to further explore the BBH model or to use this tutorial to study other populations. The next tutorials show you how to select GRBs, perform a SFH calculation at a single metallicity, and how to run an individual binary.