Stellar Transient Populations 🔍

In this notebook, we will cover: 1. Selecting binaries and combining metallicities 2. Exploring the TransientPopulation class

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

If you haven’t done so yet, export the path POSYDON environment variables. For example:

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

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'])

import pandas as pd

pop = Population('BBH_contact.h5')

pop.oneline

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 is empty, no BBH mergers were in your population.

Lets 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

PATH_TO_POSYDON_DATA = '/Users/max/Documents/POSYDON_data/240305/POSYDON_data/tutorials/population-synthesis/example/'

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']


pop = Population(PATH_TO_POSYDON_DATA+files[0])

pop.history

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(PATH_TO_POSYDON_DATA+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, 'BBH_contact.h5', append=True)

If you 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_met property contains information about all the metallicities in the file.

You can even combined multiple runs at the same metallicity together, if you like. This will combine their simualted and underlying masses.

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

BBH_pop.mass_per_metallicity

Transient population

Although we now have selected all binaries with a BBH merger, we don’t have the extact moment of the merger yet.

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”. 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: 1. BBH_selection_function 2. 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):
    '''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=10 when you initialise the Population, like we’ve done above for BBH_pop, each chunk will contain 10 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.

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.

Lets 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']
    df_transients['S1_spin_orbit_tilt'] = oneline_chunk['S1_spin_orbit_tilt']
    df_transients['S2_spin_orbit_tilt'] = oneline_chunk['S2_spin_orbit_tilt']
    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 infomration 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
import pandas as pd

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']
    df_transients['S1_spin_orbit_tilt'] = oneline_chunk['S1_spin_orbit_tilt']
    df_transients['S2_spin_orbit_tilt'] = oneline_chunk['S2_spin_orbit_tilt']
    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'], oneline_chunk['S2_spin_orbit_tilt'])

    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. First, we calculate them for the population:

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']
    df_transients['S1_spin_orbit_tilt'] = oneline_chunk['S1_spin_orbit_tilt']
    df_transients['S2_spin_orbit_tilt'] = oneline_chunk['S2_spin_orbit_tilt']
    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'], oneline_chunk['S2_spin_orbit_tilt'])

    # 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
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 somehting 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 TransienPopulation is stored in file, you can continue your analysis by opening the file using the TransientPopulation class. You only need to do your selection once.

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.

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 matplotlib.pyplot as plt
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 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)

The rates can also be accessed using the Rates class. You will need to provide your transient_name and SFH_identifier. You can keep as many transients and SFHs in your file as you like.

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 SNe 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 1. TransientPopulation 2. The z_events 3. The weights

Using these, new weights are calculated. The BBH analysis comes with a detection function for several different 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 <>`__.