The Cosmic History of Long Gamma Ray Bursts

The cosmic formation rate of long Gamma Ray Bursts (LGRBs) encodes the evolution, across cosmic times, of their progenitors' properties and of their environment. The LGRB formation rate and the luminosity function, with its redshift evolution, are derived by reproducing the largest set of observations collected in the last four decades, namely the observer-frame prompt emission properties of GRB samples detected by the Fermi and Compton Gamma Ray Observatory (CGRO) satellites and the redshift, luminosity and energy distributions of flux-limited, redshift complete, samples of GRBs detected by Swift. The model that best reproduces all these constraints consists of a GRB formation rate increasing with redshift $\propto (1+z)^{3.2}$, i.e. steeper than the star formation rate, up to $z\sim3$ followed by a decrease $\propto(1+z)^{-3}$. On top of this, our model predicts also a moderate evolution of the characteristic luminosity function break $\propto(1+z)^{0.6}$. Models with only luminosity or rate evolution are excluded at $>5\sigma$ significance. The cosmic rate evolution of LGRBs is interpreted as their preference to occur in environments with metallicity $12+\log(\rm O/H)<8.6$, consistently with theoretical models and host galaxy observations. The LGRB rate at $z=0$, accounting for their collimation, is $\rho_0=79^{+57}_{-33}$ Gpc$^{-3}$ yr$^{-1}$ (68% confidence interval). This corresponds to $\sim$1\% of broad-line Ibc supernovae producing a successful jet in the local Universe. This fraction increases up to $\sim$7% at $z\ge3$. Finally, we estimate that at least $\approx0.2-0.7$ yr$^{-1}$ of Swift and Fermi detected bursts at $z<0.5$ are jets observed slightly off-axis.


INTRODUCTION
Distinctive features of Long durantion Gamma Ray Bursts (LGRBs) such as their detection up to high redshifts (z=8.2-9.2 - Salvaterra et al. 2009, Tanvir et al. 2009, Cucchiara et al. 2011, their association with faint blue star forming galaxies, their positional coincidence with the brightest UV regions within their hosts (Fruchter et al. 2006) and their association with core-collapase supernovae (Hjorth & Bloom 2011) point to a massive star progenitor thus establishing a direct link with the formation and relatively rapid evolution of massive stars across cosmic times. The study of increasingly larger samples of hosts revealed the preference of LGRBs to form in sub-solar metallicity environments (Vergani et al. 2015;Perley et al. 2016;Palmerio et al. 2019). This is in agreement with progenitors models (MacFadyen & Woosley 1999) requiring a low metallicity to keep high specific angular momentum at the time of the core-collapse in order to produce the GRB jet (Yoon et al. 2006). Therefore, the physical conditions for the production of a relativistic jet determine the intrinsic properties of the LGRB population such as their characteristic luminosities and cosmic rate. GRB population studies aims at reconstructing the population luminosity and formation rate functions. The free parameters of these functions are often constrained by reproducing a set of rest frame and observer frame properties (Daigne et al. 2006;Wanderman & Piran 2010;Robertson & Ellis 2012) of detected GRB samples and by properly accounting for sample incompleteness Ghirlanda et al. 2015;Palmerio & Daigne 2020).
In this work we build a GRB population model ( §2) and derive the cosmic formation rate of LGRBs and their luminosity function together with its evolution with redshift. To this aim we exploit, for the first time, the largest set of observational constraints collected by different GRB detectors in the last four decades ( §3). We account for the presence of relativistic jets and for the orientation-dependence of the observed GRB properties. We derive the evolution of the ratio of GRB-to-cosmic-star formation rate and test its consistency with the host metallicity bias ( §4). We discuss implications of our population model in §5. Throughout the paper we assume standard flat cosmology Ω M = 0.3, h = 0.7.

THE MODEL
We build a synthetic LGRB population model based on two distribution functions. The GRB formation rate per unit comoving volume ρ(z) is modelled with the parametric function (Cole et al. 2001;Madau & Dickinson 2014): where ρ 0 , in Gpc −3 yr −1 , is the LGRB local density rate. The choice of this functional form is motivated by its relatively small number of free parameters and, for the ease of comparison, by its systematic use to fit the Cosmic Star Formation Rate (CSFR - Madau & Dickinson 2014;Madau & Fragos 2017;Hopkins & Beacom 2006). Flat priors are considered for the the free parameters (p z,1 , p z,2 , p z,3 ) which can vary over wide ranges including also the values corresponding to the CSFR. This ensures that we can recover, if present, any significant deviation of the GRB formation rate from the CSFR. The luminosity function of LGRBs is often described in the literature by a broken power-law (e.g. Salvaterra et al. 2009;Salvaterra et al. 2012): where a 1 and a 2 are the indices of the power laws respectively below and above the break luminosity L iso,b which is left free to evolve with redshift ∝ (1 + z) δ . Therefore, L iso,0 represents the break of the luminosity function at z = 0. Φ(L iso , z) is normalised to its integral and it is defined for L iso ≥ 10 47 erg/s. As alternative working hypothesis sometimes considered in the literature is a Schechter (Schechter 1976) luminosity function. We discuss in §A.3 the results obtained with under this different assumption. In summary, our model allows for the evolution of the GRB population both in terms of rate density (Eq1) and luminosity (Eq.2).
LGRBs follow the E p − L iso (Yonetoku et al. 2004) correlation, where E p is the peak energy of the νF ν spectrum, and the E p − E iso correlation (Amati & Frontera 2002), linking the peak energy with the isotropic energy E iso . It has been proved (Nava et al. 2012;Dai et al. 2021) that these correlations do not evolve with redshift. Following Ghirlanda et al. 2015, in our simulation we sample Eq.1 extracting a redshift and independently we assign an E p value sampling a broken powerlaw. The values of L iso and E iso corresponding to E p are derived from the E p − L iso and E p − E iso correlations, respectively, accounting for their gaussian scatters (Nava et al. 2012). For the ease of comparison with the literature, once we have derived the posterior parameter distributions of the model, we map, through the E p − L iso correlation, the peak energy distribution function back to the luminosity function of which we provide the resulting parameter values in Tab.1. GRB rest frame durations are estimated as f · E iso /L iso , where f is a scaling factor sampled from a free gaussian distribution G[f |µ, σ].
It is known that GRBs are highly collimated (e.g. Ghirlanda et al. 2007). Therefore, we compute their intrinsic rate ρ(z) (Eq.1) by accounting for the elusive bulk of the population consisting of events with jets oriented off-side our line of sight. We assume uniform conical jets distributed with random orientations θ v with respect to the line of sight. Jet opening angle values θ jet appear correlated with E iso (Ghirlanda et al. 2007) such that we assume a free parametric correlation θ jet -E iso . Even off-axis jets (θ v >θ jet ) may be detected depending on their orientation and bulk   Lorentz factor Γ, characterising the prompt emission phase. We account, for the first time within a population study of LGRBs, for the relativistic effects on the observed flux, fluence and peak energy induced by the jet orientation and beaming through the doppler factor δ = 1/Γ(1 − β cos θ v ) (see e.g. Ghisellini et al. 2006). To this aim we sample the Γ-E iso correlation reported by Ghirlanda et al. 2018 to assign Γ values to simulated events. The GRB prompt emission spectrum is described with the empirical Band function (Band et al. 1993) consisting of two smoothly joined power laws. To each simulated GRB we assign a low and high energy spectral index extracted from gaussian distributions G[α|µ, σ] and G[β|µ, σ] respectively, whose parameters (e.g. Goldstein et al. 2012) are given in Tab.1.
The population model has 12 free parameters (Tab.1): [ρ 0 , p z,1 , p z,2 , p z,3 ] describe the GRB cosmic formation rate (Eq.1) and [a 1 , a 2 , L iso,0 , δ] define the luminosity function (Eq.2); the remaining four parameters (µ, σ) and (K, S) describe the stretching factor of the duration distribution and the E iso -θ jet correlation respectively. All other parameters reported in Tab.1 are fixed to their presently known values and are derived from the corresponding references given above.
The twelve free parameters are constrained by matching a set of observed distributions C J ( §3) and by maximising a binned likelihood. The probability that the model distribution has m i objects in the i bin of the distribution, given that there are r i real events in this bin, follows a Poisson distribution. The total log-likelihood can be expressed as: where the internal sum is performed over the n bins of each distribution of the constraint C J (with J = 1, ..., 14). The observational constraint distributions C J (described in §3) are build from the observed properties of LGRB samples detected by different instruments and are represented by different simbols in the panels of Fig.1. The posterior distributions of the 12 free parameters are derived through a Markov Chain Monte Carlo (MCMC) method which samples from wide uniform prior distributions of the free parameters:   Fig.1), the observer-frame peak energy (E p /1 + z) and duration (T 90 ) (panels b and d in Fig.1). In order to minimize the sample incompleteness, particularly severe in the faint end of the peak flux distribution 1 , we consider LGRBs with F P ≥4 ph cm −2 s −1 , integrated over the 10 keV -1 MeV energy range for Fermi /GBM and 20 keV -2 MeV for GCRO/BATSE. This selection results in 904 GRBs (out of 1906) for Fermi /GBM (red symbols in Fig.1) and 370 (out of 1529) for CGRO/BATSE (blue symbols in Fig.1). The inset in Fig.1-(a) shows the peak flux (50-300 keV) of 2810 CGRO/BATSE bursts in the Stern et al. (2001) For the redshift distribution (panel e) we consider two well selected Swift /BAT flux-limited samples for which the large majority of bursts have a redshift measurement. The BAT6 sample Pescalli et al. 2016) consists of Swift LGRBs with 15-150 keV peak flux ≥ 2.6 ph cm −2 s −1 . The SHOALS sample (Perley et al. 2016) comprising Swift GRBs with a 15-150 keV integrated fluence 10 −6 erg cm −2 . Both these samples have 80-90% redshift measurements either from afterglow or host spectroscopy. The BAT6 sample also provides the constraints on E iso and L iso distributions (panel f). Finally, we consider the jet opening angle distribution (panel g) of GRBs collected in Ghirlanda et al. 2007.
Further we also verified the consistency of our best model (Fig.5 in §A.1) with the distributions of the peak flux, fluence and duration (T 90 ) of the samples of LGRBs observed by the Swift /BAT 2 and by the Gamma Ray Burst Monitor (GRBM) on board the Beppo/SAX satellite (Frontera et al. 2008).

RESULTS
The model which best matches the set of observational constraints is shown by the solid histograms in Fig.1 (and  Fig.5). The shaded/dashed regions around the model line represent the 68% model uncertainty. This is computed by The luminosity function has a faint-end slope ∼ −1 and a bright-end slope ∼ −2.2 with a break luminosity at z=0 L b,0 = 10 52 erg s −1 . These values are consistent with those reported by previous studies Wanderman & Piran 2010). We find a mild evolution of this break which increases with redshift as L b = L b,0 (1+z) 0.64 .
The cosmic rate of LGRBs, ρ(z), is shown in Fig.3-(a) (solid red line). The shaded orange region represents the 1σ uncertainty as obtained by resampling the posterior distributions of the free parameters. ρ(z) increases with redshift ∝ (1 + z) 3.2 , it peaks at z ∼ 3 and declines at larger redshifts ∝ (1 + z) −3 . These slopes are respectively steeper and similar to those of the CSFR 3 of Madau & Dickinson 2014 (solid blue line in Fig.3-(a)).
The local LGRB formation rate is ρ 0 = 79 +57 −33 Gpc −3 yr −1 (68% confidence interval). This is derived by normalizing the model population: we consider the ∼ 170±13 bright GRBs with peak flux in the 15-150 keV energy range ≥2.6 ph cm 2 s −1 detected by Swift /BAT in 8 years. We assume for Swift/BAT an average duty cycle of 0.75 and the fully coded field of view of 1.4 sr. We verified that our results are unchanged if we renormalize the population to the bright GRB sample detected by Fermi /GBM, namely 902 LGRBs with P(10-1000)≥ 5 ph cm 2 s −1 detected in 11.8 years within its 8.8 sr field of view, also accounting for a 0.55 duty cycle. Since our population accounts for the random orientation of GRB jets in the sky, the derived local rate represents the intrinsic one.
We performed the same analysis also under the assumption of a Schechter luminosity function and the results are reported in §A.2. We find that the shape of the LGRB rate at low redshifts is unchanged (Fig.6). The Schechter scenario produces a slightly shallower GRB formation rate decrease at high redshifts. This is however consistent with the CSFR given the uncertainties on high redshift measurements.

DISCUSSION AND CONCLUSIONS
The LGRB efficiency as a function of redshift, defined as the ratio between the GRB rate and the cosmic star formation rate, is shown in Fig.3-b by the solid red line. Its 1σ uncertainty (shaded violet region in Fig.3) is computed by resampling the posterior distributions of the free parameters shown in Fig.2 10 5 times. The uncertainty For the CSFR we have considered the function of Madau & Dickinson 2014 (solid blue line in Fig.3-a). In the local Universe this ratio is ∼ 5 × 10 −5 M −1 ⊙ and it increases up to a factor of 10 at z ∼ 3. In Fig.3-b we show with the blue and red symbols respectively the fraction of star formation leading to Supernovae Ibc and broad line (BL) Ic (Perley et al. 2020) assuming the core-collapse supernova rate derived by Strolger et al. 2015. We find that, in the local Universe, 1.3 +1.0 −0.6 % of SNIc BL produce a jet and the associated GRB event. This fraction increases with redshift up to ∼7% at z > 3.
The general trend of the GRB cosmic rate at low redshifts shows a decrease of the efficiency of producing GRBs with decreasing redshifts (Fig.3). This result argues against the claims of a increase of the LGRB efficiency towards lower redshifts (so called low redshift excess - Yu et al. 2015;Petrosian et al. 2015;Tsvetkova et al. 2017;Lloyd-Ronning et al. 2019, but see Bryant et al. 2021;Pescalli et al. 2016;Le et al. 2020) obtained through non-parametric methods. The low redshift GRB uptick could be partly absorbed if the jet opening angle evolves with redshift (Lloyd-Ronning et al. 2020).
The GRB efficiency redshift evolution can be interpreted as due to the increase of the metallicity with cosmic time which prevents a larger fraction of massive stars from producing a GRB (Yoon et al. 2006). In order to verify this possibility we compute the fraction of cosmic star formation in galaxies with average metallicity below a given threshold following the method of Robertson & Ellis 2012. We combine ( §A.3) the mass-metallicity relation and its redshift evolution (Maiolino et al. 2008) with the galaxy mass function (McLeod et al. 2021) and the galaxy main sequence (Tomczak et al. 2014). We find that a metallicity threshold 12 + log(O/H) ≤ 8.6 , represented by the the dashed blue models in Fig.3, accounts for the derived GRB efficiency evolution. The metallicity thresholds we derive are consistent with the observed GRB host galaxy metallicity (Vergani et al. 2015;Palmerio et al. 2019). To further verify this consistency we computed the stellar masses of galaxies with metallicity below the threshold value found and compare with the observed masses of a complete sample of GRB hosts (Vergani et al. 2015;Palmerio et al. 2019 -black symbols in Fig.4 left panel). The maximum and average stellar masses (green and blue curves in Fig.4 Previous LGRB population studies (e.g. Firmani et al. 2004;Salvaterra et al. 2012;Palmerio & Daigne 2020) showed that a pure density evolution (PDE) of the GRB formation rate, where the luminosity function is invariant with redshift, and pure luminosity evolution (PLE), where the GRB formation rate is proportional to the CSFR, cannot be distinguished on the basis of the GRB peak flux and luminosity-redshift distributions. Our best model is intermediate between these two extreme scenarios and, indeed, we find both the GRB formation rate evolving with redshift ( Fig.3-b) and a mild evolution of the characteristic GRB luminosity ∝ (1 + z) 0.64 (Tab.1). We at more than 5σ significance the extreme cases of PLE and PDE. Moreover, the PLE is also excluded by the host galaxy mass distribution (cf red line in Fig.4). The finding of an evolution of the LGRB characteristic luminosity ∝ (1 + z) 0.64 could hide the change of some properties of the progenitors with redshifts possibly ruled by their metallicity (Yoon et al. 2006) which could induce a flattening of the initial mass function (Fryer et al. 2021).
This possibility can be tested by an increasing fraction of GRBs detected at high redshifts as envisaged by the next generation missions optimised to this purpose (e.g. the Transient and High Energy Sky and Early Universe Surveyor -THESEUS Amati et al. 2021or Gamow Explorer -White et al. 2020).
Our population model, which implements the jet opening angle of GRBs and accounts for the observed properties as a function of the viewing angle, allows us to estimate the fraction of GRBs detected by Fermi /GBM and Swift /BAT which are observed off-axis, i.e. θ v >θ jet . Fig.4 (right panel) shows the expected cumulative rate as a function of redshifts of Fermi /GBM (blue filled stripe) and Swift /BAT (orange filled stripe). The hatched curves show the detection rates of off-axis events according to our population. We estimate that ∼10% of LGRBs detected by Fermi and Swift at z < 0.5, corresponding to a rate 0.2-0.7 yr −1 respectively, should be observed off-axis. This fraction reduces to less than 3% at z < 2. In this work we adopted a uniform jet structure consistently with what suggested by the study of the luminosity function of long GRBs . Therefore, our results should be considered as a conservative lower limit since a universal structured jet should produce a larger fractions (e.g. Salafia et al. 2015) of detected off-axis events.
We acknowledge support by ASI-INAF agreement n. 2018-29-HH-0. We acknowledge PRIN-MIUR 2017 (grant 20179ZF5KS)  In order to verify the consistency of our model with GRB samples detected by other instruments we considered the full Swift /BAT and BeppoSAX/GRBM (Frontera et al. 2008) samples. We applied a peak flux cut extracting 536 Swift /BAT GRBs with a 15-150 keV peak flux F p ≥ 0.5 phot cm −2 s −1 . Their observer frame distributions are shown by the green asterisks in the top panels of Fig.5. For Beppo/SAX we considered the 118 GRBs with peak flux, integrated in the 40-700 keV energy range, F p ≥ 10 −6 erg cm −2 s −1 . The BeppoSAX distributions are shown by the pink asterisks in the bottom panels of Fig.5. The population model is shown by the solid histogram and shaded regions.

A.2. Schechter luminosity function
We considered the possibility that the GRB luminosity function is described by a Schechter function (Schechter 1976), namely a powerlaw with a high-luminosity exponential cutoff. This model has two free parameters, the slope of the powerlaw α and the cutoff luminosity L c . We fix the minimum luminosity of this function to 10 47 erg s −1 . With this model we can reproduce all the observational constraints described in §3. We obtain α = 1.15 ± 0.15 and log(L c ) = 53.21 ± 0.28. The resulting GRB formation rate is shown in Fig.6  (2) of Maiolino et al. 2008. The values of K(z) and M 0 (z) are interpolated from those reported in Tab.5 of their paper. We also corrected (through the 0.04 term) for the different assumption of a Chabrier initial mass function in host galaxy mass estimates shown in Fig.4. Figure 5. Consistency checks of the population model with the distributions of the peak flux, fluence and duration of GRBs detected by the Swift/BAT satellite (top row) and Beppo-SAX (bottom row). Swift/BAT GRBs (green symbols) are obtained selecting GRBs with a peak flux (15-150 keV) ≥0.5 ph cm −2 s −1 while Beppo-SAX GRBs (pink symbols) are selected with peak flux (40-700 keV) ≥ 10 −6 erg cm −2 s −1 . Model is shown by the solid lines and its 1σ uncertainty (described as in Fig.1) by the solid filled histogram.