Computing Long-Duration Gamma-Ray Bursts from Double Compact Object Populations 🌌

Creating the GRB TransientPopulation

We will use the same BBH population as in the Stellar Transient Populations tutorial to compute the long gamma-ray bursts (LGRB) rate associated with the formation of merging BBHs.

Make sure the BBH_contact.h5 file is present to continue with this tutorial. You can find a copy of it in the example dataset.

from posydon.popsyn.synthetic_population import Population

file = 'BBH_contact.h5'

BBH_population = Population(file)
BBH_population.mass_per_metallicity

BBH_population.calculate_formation_channels(mt_history=True)

We had already selected only binaries that result in a BBH merger in the BBH_contact.h5 file, so we don’t need to perform any selection on the binaries.

Thus, our next step is to create a selection function for LGRBs based on the properties in the history and onelines of the binaries.

POSYDON already includes a basic selection function for GRBs named GRB_selection in posydon.popsyn.transient_select_func, which identifies if a GRB occured based on the presence of m_disk_radiated for either star. GRB_selection also outputs some pre and post supernova properties.

However, as you can see below, the GRB_selection requires an additional input parameter of S1_S2. The Population.create_transient_population function cannot pass these additional arguments to the sub-function, GRB_selection. Thus, we need to wrap the GRB_selection function.

This is useful, if you would like to create a similar population, but use different model parameters.

from posydon.popsyn.transient_select_funcs import GRB_selection

GRB_selection(BBH_population.history[10],
              BBH_population.oneline[10],
              BBH_population.formation_channels.loc[[10]],
              S1_S2='S2')

import pandas as pd

def GRB_wrapper(transient_chunk, oneline_chunk, formation_channels_chunk):
    """calculate the GRBs from star 1 and star 2 for a given chunk of the population"""
    # select potential GRBs from star 1
    df_1 = GRB_selection(transient_chunk, oneline_chunk, formation_channels_chunk, S1_S2='S1')
    # select potential GRBs from star 2
    df_2 = GRB_selection(transient_chunk, oneline_chunk, formation_channels_chunk, S1_S2='S2')

    # combine the two dataframes
    if df_1 is not None and df_2 is not None:
        GRB_df = pd.concat([df_1, df_2])
    elif df_1 is not None:
        GRB_df = df_1
    elif df_2 is not None:
        GRB_df = df_2
    else:
        return None

    return GRB_df

LGRB = BBH_population.create_transient_population(GRB_wrapper, 'LGRB')

You’ve now created a LGRB population, where either the first and/or second star has some amount of radiative disk.

This is, of course, just an example and the actual amount of disk required to power a LGRB will be higher than some of the values included in our “LGRB” rate. Moreover, the number of binaries in this population is insufficient to actually see the all unique events.

For now, let’s continue and calculate the metallicity bias function for the LGRBs. We don’t need to calculate the underlying_mass of the population, because we already did this in the BBH tutorial.

LGRB.get_efficiency_over_metallicity()

LGRB.plot_efficiency_over_metallicity(channels=True)

After this selection, we follow similar steps to the BBH analysis for creating and plotting a cosmic star formation history weighted rate.

If you’ve followed the previous tutorial on the BBH analaysis, you will still have access to those populations too.

LGRB_rates = LGRB.calculate_cosmic_weights('IllustrisTNG')

LGRB_rates.calculate_intrinsic_rate_density(mt_channels=True)

from posydon.popsyn.synthetic_population import Rates
BBH_rates  = Rates('BBH_contact.h5', 'BBH', 'IllustrisTNG')

BBH_rates.intrinsic_rate_density

Below we plot the BBH rate and the un-normalised LGRB rate together, showing how you can access multiple event rates, and transient populations from the same file.

import matplotlib.pyplot as plt
plt.plot(BBH_rates.centers_redshift_bins, BBH_rates.intrinsic_rate_density['total'], label='BBH')
plt.plot(LGRB_rates.centers_redshift_bins, LGRB_rates.intrinsic_rate_density['total'], label='LGRB')

plt.yscale("log")
plt.xlim(0,10)
plt.legend()
plt.show()

LGRB_rates.population

fig, ax = plt.subplots(1,1)
LGRB_rates.plot_hist_properties('S1_spin_postSN', intrinsic=True, label='S1', color='blue', ax=ax, show=False)
LGRB_rates.plot_hist_properties('S2_spin_postSN', intrinsic=True, label='S2', color='orange', ax=ax, show=False)
ax.legend()
ax.set_xlabel('BH spin post LGRB')
ax.set_ylabel('\# Events in bin')
plt.show()

It’s possible to calculate more properties for the LGRBs using the posydon.popsyn.GRB module.

It contains the get_GRB_properties that can provide more detailed GRB properties based on empirical or simulation relations for the LGRB energy and beaming factor.

Enjoy exploring this function and creating a more detailed GRB population! Keep in mind that it’s not possible to add additional columns to the transient population without overwritting the population.

The next tutorial will focus on rerunning specific binaries in a population or setting up an unique binary yourself.