Short GRBS

In Ghirlanda et al. 2016 a fitting algorithm was used to determine the redshift and luminosity of short GRBS. We can use the parameters to reproduce the population and the observed GBM survey.

from popsynth import SFRDistribution, BPLDistribution, PopulationSynth, NormalAuxSampler, AuxiliarySampler, HardFluxSelection
from popsynth import update_logging_level
update_logging_level("INFO")

%matplotlib inline

import matplotlib.pyplot as plt
from jupyterthemes import jtplot

jtplot.style(context="notebook", fscale=1, grid=False)
purple = "#B833FF"
yellow = "#F6EF5B"

import networkx as nx
import numpy as np
import warnings

warnings.simplefilter("ignore")

In the work, the luminosity function of short GRBs is model as a broken power law.

bpl = BPLDistribution()

bpl.alpha = -0.53
bpl.beta = -3.4
bpl.Lmin = 1e47 # erg/s
bpl.Lbreak = 2.8e52
bpl.Lmax = 1e55

To model the redshift distribution, an empirical form from Cole et al 2001 is used. In popsynth we call this the SFRDistribution (but perhaps a better name is needed).

sfr = SFRDistribution()

sfr.r0 = 5.
sfr.a = 1
sfr.rise = 2.8
sfr.decay = 3.5
sfr.peak = 2.3

We can checkout how the rate changes with redshift

fig, ax = plt.subplots()

z = np.linspace(0,5,100)

ax.plot(z, sfr.dNdV(z), color=purple)
ax.set_xlabel("z")
ax.set_ylabel(r"$\frac{\mathrm{d}N}{\mathrm{d}V}$")

Text(0, 0.5, '$\\frac{\\mathrm{d}N}{\\mathrm{d}V}$')

In their model, the authors also have some secondary parameters that are connected to the luminosity. These are the parameters for the spectrum of the GRB. It is proposed that the spectra peak energy (Ep) is linked to the luminosity by a power law relation:

\[\log E_{\mathrm{p}} \propto a + b \log L\]

We can build an auxiliary sample to simulate this as well. But we will also add a bit of scatter to the intercept of the relation.

intercept =NormalAuxSampler(name="intercept", observed=False)

intercept.mu = 0.034
intercept.sigma = .005

class EpSampler(AuxiliarySampler):

    _auxiliary_sampler_name = "EpSampler"

    def __init__(self):

        # pass up to the super class
        super(EpSampler, self).__init__("Ep", observed=True, uses_luminosity = True)

    def true_sampler(self, size):

        # we will get the intercept's latent (true) value
        # from its sampler

        intercept = self._secondary_samplers["intercept"].true_values

        slope = 0.84

        self._true_values = np.power(10., intercept + slope * np.log10(self._luminosity/1e52) + np.log10(670.))

    def observation_sampler(self, size):

        # we will also add some measurement error to Ep
        self._obs_values = self._true_values + np.random.normal(0., 10, size=size)

Now we can put it all together.

pop_synth = PopulationSynth(spatial_distribution=sfr, luminosity_distribution=bpl)

We will have a hard flux selection which is Fermi-GBM’s fluz limit of ~ 1e-7 erg/s/cm2

selection = HardFluxSelection()
selection.boundary = 1e-7

pop_synth.set_flux_selection(selection)

We need to add the Ep sampler. Once we set the intercept sampler as a secondary it will automatically be added to the population synth.

ep = EpSampler()

ep.set_secondary_sampler(intercept)

pop_synth.add_auxiliary_sampler(ep)

 INFO     |  registering auxilary sampler: Ep 

We are ready to sample our population. We will add some measurement uncertainty to the fluxes as well.

population = pop_synth.draw_survey(flux_sigma=0.2)

 INFO     |  The volume integral is 797.5506658438255 

 INFO     |  Expecting 770 total objects 
 INFO     |  Sampling: Ep 
 INFO     |  Ep is sampling its secondary quantities 
 INFO     |  Sampling: intercept 
 INFO     |  applying selection to fluxes 
 INFO     |  Detected 499 distances 
 INFO     |  Detected 499 objects out to a distance of 5.31 

population.display_fluxes(true_color=purple, obs_color=yellow, with_arrows=False, s= 5);

Let’s look at our distribution of Ep

fig, ax = plt.subplots()

ax.hist(np.log10(population.Ep.selected), histtype="step", color=yellow, lw=3, label="Ep observed")
ax.hist(np.log10(population.Ep.non_selected), histtype="step", color=purple, lw=3,  label="Ep hidden")
ax.set_xlabel("log Ep")

ax.legend()

<matplotlib.legend.Legend at 0x7f14b9a12190>

fig, ax = plt.subplots()


ax.scatter(population.fluxes.non_selected,
           population.Ep.non_selected,c=purple, alpha=0.5)
ax.scatter(population.fluxes.selected,
           population.Ep.selected,c=yellow, alpha=0.5)

ax.set_xscale("log")
ax.set_yscale("log")

ax.set_xlabel("log Ep")
ax.set_xlabel("log Flux")

Text(0.5, 0, 'log Flux')

Does this look like the observed catalogs?