© J. T. Palmerio & F. Daigne 2021

1. Introduction

Gamma-ray bursts (GRBs) are the most powerful electromagnetic phenomena and are associated with the relativistic ejection following the birth of a stellar mass compact object (black hole or magnetar; see e.g., Piran 2005; Kumar & Zhang 2015). Gamma-ray bursts exhibit two main radiative phases: (i) the prompt emission, which is commonly detected in the hard X-rays and soft γ-rays (1 keV–10 MeV) and usually lasts a few hundreds of milliseconds to a few hundreds of seconds; and (ii) the afterglow emission, which is initially bright, decays rapidly, and peaks successively in X-rays, optical and radio (see e.g., Gehrels et al. 2009; Gehrels & Razzaque 2013). Gamma-ray bursts have proven themselves to be powerful probes of the Universe owing to some unique properties. They form up to very high redshift (spectroscopically confirmed at z ∼ 8.2 by Salvaterra et al. 2009 and Tanvir et al. 2009, and using photometry only at z ∼ 9 by Cucchiara et al. 2011) and can be detected thanks to their γ-ray light, which is largely unaffected by dust and hydrogen absorption. The transient, fading nature of their afterglow benefits from time dilation due to cosmological expansion (Lamb & Reichart 2000); 1 day after the prompt emission on Earth corresponds to 6 h in the source frame at z = 3 and 2 h at z = 10, essentially catching the afterglow earlier in its light curve (and thus brighter) as redshift increases.

Among the various classes of GRBs, the case of long GRBs (LGRBs, which have a duration of the prompt emission longer than ∼2 s; Mazets et al. 1981; Kouveliotou et al. 1993) is the most promising for the study of the distant Universe. These are the most frequent GRBs, and both theoretical progenitor models (‘collapsar’ model; Woosley 1993; Paczynski 1998) and observations have firmly associated them with the collapse of certain massive stars. Indeed, they generally occur in faint, blue, low-mass star-forming galaxies (Le Floc’h et al. 2003; Savaglio et al. 2009; Palmerio et al. 2019) and in the brighter regions of their hosts (Fruchter et al. 2006; Svensson et al. 2010; Blanchard et al. 2016; Lyman et al. 2017). In addition, the majority of low-redshift LGRBs are found in association with core collapse supernovae (Bloom et al. 1999; Hjorth et al. 2003; Krühler et al. 2017; Izzo et al. 2017, 2020; Melandri et al. 2019). Due to the short-lived nature of their massive star progenitors, LGRBs are expected to occur up to very high redshift, possibly in association with the first generation of stars (Bromm & Loeb 2006), and thus can be used as lighthouses to study galaxies (e.g., Le Floc’h et al. 2006; Perley et al. 2013; Vergani et al. 2017) and the intergalactic medium at high redshift (e.g., Fynbo et al. 2006; Hartoog et al. 2015; Selsing et al. 2019; Tanvir et al. 2019; Vielfaure et al. 2020).

From their association with certain massive stars, LGRBs as a population are considered potential tracers of star formation (e.g., Lamb & Reichart 2000; Porciani & Madau 2001; Kistler et al. 2008; Robertson & Ellis 2012; Vangioni et al. 2015). However, the precise link between star formation and LGRBs remains shrouded in uncertainty for two main reasons. The first reason is that this link depends on many factors that are sometimes poorly constrained and which may evolve with redshift, such as the stellar initial mass function (IMF), the properties of the progenitor star (e.g., mass range, metallicity distribution, initial rotation distribution), the fraction of binary progenitors (see e.g., Chrimes et al. 2020), or the LGRB luminosity, spectrum, and jet opening angle. The second reason is the difficulty in obtaining large, unbiased, and redshift-complete observational samples of LGRBs as they require a rapid response from ground follow-up and significant telescope time to get meaningful statistics. In practice, we often have to choose between biased, incomplete or smaller samples (see e.g., Hjorth et al. 2012; Salvaterra et al. 2012; Krühler et al. 2015; Perley et al. 2016a, for the case of smaller but unbiased and highly complete samples).

The goal of population models is to overcome the limitations of biased or incomplete samples by modelling the underlying intrinsic population and fitting it to carefully selected observational samples. This allows one to then constrain characteristics of the intrinsic population, such as the redshift distribution, the rate, the evolution of the energetics with redshift, or the correlations between physical parameters. Over the last 20 years, population models have become an effective approach to study GRBs due to growing sample sizes and increasingly affordable computing power. Porciani & Madau (2001) studied the luminosity function of LGRBs using the peak flux distribution of the Burst And Transient Source Experiment on board Compton Gamma Ray Observatory (CGRO/BATSE) from Kommers et al. (2000) and estimated the rate of LGRBs under the hypothesis that their redshift distribution follows the cosmic star formation rate density (CSFRD). They used a fixed Band function (Band et al. 1993; see Appendix A) with parameters α = −1, β = −2.25, and E_p = 511 keV for the LGRB spectrum. Assuming a constant LGRB production efficiency from stars, they found (1 − 2) × 10⁻⁶ LGRBs pointing towards us per core collapse. Daigne et al. (2006) used a similar approach, but in addition to including observations of the peak energy distribution of CGRO/BATSE LGRBs, they allowed for three scenarios regarding the redshift distribution of LGRBs. These scenarios correspond to three different hypotheses about the LGRB production efficiency: constant, mildly increasing, and strongly increasing with redshift. They also used a Band spectrum but added realistic distributions for the values of α and β as well as two different scenarios for the source-frame peak energy distribution: log-normal or correlated to the peak luminosity. Using the observed redshift distribution from the early Swift results of Jakobsson et al. (2006) as a crosscheck, they found evidence for an LGRB production efficiency increasing with redshift. With the advent of Swift and thanks to dedicated follow-up campaigns, samples of GRBs with a significant fraction of redshift measurements were obtained, allowing for two new approaches to emerge. The first approach is to model the redshift recovery efficiency as in Wanderman & Piran (2010). These authors directly inverted the observed luminosity-redshift distribution of a sample of 120 Swift LGRBs, assuming a fixed Band spectrum as in Porciani & Madau (2001), to obtain the intrinsic luminosity function and redshift distribution. They found some evidence that the LGRB rate does not follow the CSFRD, in agreement with Daigne et al. (2006). Using a different approach, Salvaterra et al. (2012) designed a complete and unbiased sample of 58 bright Swift LGRBs named BAT6 (later extended to 99 LGRBs by Pescalli et al. 2016), essentially paying the cost of sample size in order to achieve high redshift completeness and avoid the uncertainties caused by modelling the redshift recovery efficiency. They found evidence for strong redshift evolution of the luminosity function, although they noted the degeneracy with the redshift evolution of the LGRB efficiency and suggested both might be at play.

Building on these works, we present in this paper our own population model, which combines the latest, largest, and most diverse datasets of LGRBs available to date with currently relevant scenarios in order to investigate some of the questions raised above. We fit our population model parameters using an intensity constraint based on a large sample of CGRO/BATSE LGRBs, including faint events, a spectral constraint based on a sample of bright Fermi/GBM LGRBs with well-measured spectral parameters, and a redshift constraint based on the extended BAT6 sample. Finally, we discuss the implications of our results in terms of intrinsic spectral-energetics correlations, soft bursts, cosmic evolution of the LGRB luminosity and/or rate, the use of LGRBs as tracers of star formation at high redshift, and strategies for detecting high redshift bursts.

The paper is organised as follows. The assumptions and scenarios we use to model the intrinsic LGRB population are presented in Sect. 2, and the observations we used to constrain our model are presented in Sect. 3. Results for the best fitting scenarios are given in Sect. 4, a discussion on the implications and predictions of our model is presented in Sect. 5, and our conclusions are compiled in Sect. 6. Finally, additional relevant material, including the statistical framework used and details about our method, is provided in the appendix. All errors are reported at the 1σ confidence level unless otherwise stated. We use a flat cold dark matter (ΛCDM) cosmology¹: Ω_m = 0.27, Ω_Λ = 0.73, and H₀ = 71 km s⁻¹ Mpc⁻¹.

2. Modelling the intrinsic LGRB population

We should specify to lift any ambiguity that throughout the rest of this paper the term ‘intrinsic’ is used to qualify the entire, true LGRB population, without any selection criteria (e.g., intrinsic redshift distribution versus redshift distribution of the Swift sample) while the term ‘source-frame’ will be reserved to refer to properties as measured without the effect of cosmological redshift dilation (e.g., source-frame E_p, versus observed E_p, noted E_p obs; when not specified the default is source-frame).

Each burst in our model is described by the following properties: a redshift, a peak luminosity and a spectrum. For the majority of the work presented in this paper, these properties alone are enough to compute all the observed quantities we are interested in, and in particular the peak flux. However, for some specific studies described in Sect. 3.5, other observed quantities such as duration or fluence are necessary. For this purpose, we add two properties to fully describe an LGRB: a source-frame duration and a variability coefficient that corresponds to the ratio of the mean luminosity (averaged over the duration of the burst) over the peak luminosity. One intrinsic population is thus described by the probability density function of all these physical properties. This allows us to generate a synthetic population with a large number of bursts, N_GRB, by Monte Carlo sampling of these intrinsic distributions. We then aim to constrain the parameters of these distributions by carefully comparing the generated LGRB population to real observed samples of LGRBs. For this purpose, we generate mock samples for different instruments, based on simple selection criteria (e.g., peak flux threshold significantly above the detector threshold) to avoid any modelling of a complex instrumental efficiency, and compare to the corresponding real observed sample. Since our LGRB population is generated by Monte Carlo sampling, its accuracy depends on N_GRB; we tested different values and settled on N_GRB = 10⁶ as a compromise between population accuracy and computational time.

2.1. Properties and quantities describing one LGRB

More specifically, each LGRB in our model is described by a redshift z, four peak properties, and two additional time-integrated properties that are used only for some crosschecks in Sect. 3.5. The four peak properties are defined as occurring when the LGRB is at peak brightness and are independent of the LGRB duration; they are: (i) L, the isotropic-equivalent bolometric peak luminosity in units of erg s⁻¹, defined as $L={\int }_0^{\infty }{L}_{\text{E}}\text{d}E$, where L_E is the source-frame power density in units of [erg s⁻¹ keV⁻¹], (ii) E_p, the energy at which the E L_E spectrum peaks in the source-frame, in units of keV, (iii) α, the low-energy slope of the photon spectrum L_E/E, and (iv) β, the high-energy slope of the photon spectrum L_E/E. In addition to these four quantities, the intrinsic shape of the photon spectrum L_E/E is assumed to be given by a Band function (Band et al. 1993).

The peak brightness of an LGRB depends on the timescale with which its light curve is sampled; in the following, unless stated otherwise, this timescale is assumed to be 1.024 s in the observer frame (the typical timescale for measuring peak fluxes). This assumption is ubiquitous for these types of studies, but in reality it introduces a small error on the peak flux estimations because of time dilation; 1.024 s in the observer frame corresponds to different timescales in the source frame depending on the redshift. However, this does not significantly impact the calculated values of the peak flux; a rough estimate puts this correction between 1, 7, 13, and 18% at redshifts 0.1, 1, 3, and 6 respectively (Heussaff 2015).

The two time-integrated properties that depend on the light curve of the LGRB are: (i) T₉₀, the duration over which 5–95% of photons are emitted in the source frame, and (ii) C_var, the variability coefficient defined as the mean luminosity $\overline{L}$ divided by the peak luminosity L.

Using these aforementioned properties, it is possible to calculate the peak photon flux N^pk for any instrument observing between E_min, obs and E_max, obs (in units of [ph s⁻¹ cm⁻²]) with:

$$\begin{array}{cc}\hfill N^{\mathrm{pk}}& ={\displaystyle {\int }_{E_{\min ,\mathrm{obs}}}^{E_{\max ,\mathrm{obs}}}{N}_{E_{\mathrm{obs}}}(E_{\mathrm{obs}})\mathrm{d}{E}_{\mathrm{obs}},}\hfill \end{array}$$(1)

$$\begin{array}{cc}& =\frac{(1+z)}{4\pi D_{\mathrm{L}}^2(z)}{\displaystyle {\int }_{E_{\min ,\mathrm{obs}}(1+z)}^{E_{\max ,\mathrm{obs}}(1+z)}\frac{L_{\mathrm{E}}}{E}\mathrm{d}E,}\hfill \end{array}$$(2)

where N_Eobs is the observed photon flux density (in units of [ph s⁻¹ cm⁻² keV⁻¹]) and D_L is the luminosity distance.

Similarly we can calculate the peak energy flux F^pk for any instrument observing between E_min, obs and E_max, obs (in units of [erg s⁻¹ cm⁻²]) with:

$$\begin{array}{c}\hfill F^{\mathrm{pk}}=\frac{1}{4\pi D_{\mathrm{L}}^2(z)}{\int }_{E_{\min ,\mathrm{obs}}(1+z)}^{E_{\max ,\mathrm{obs}}(1+z)}{L}_{\mathrm{E}}\mathrm{d}E.\end{array}$$(3)

A summary of all the properties and calculated quantities is presented in Table 1. The paragraphs below are devoted to presenting in detail the assumed probability density distributions for each property, and the parameters of these distributions that we want to constrain from observations.

2.2. Intrinsic redshift distribution

The intrinsic redshift distribution is perhaps one of the most important ingredients of an LGRB population model because any constraint on this distribution has critical implications for the identification of LGRB progenitors and the use of LGRBs as tracers of star formation. However, it remains poorly constrained due to the small fraction of GRBs with measured redshifts (∼30% for Swift GRBs) leading to significant uncertainties regarding its shape; for these reasons many authors have used various functional forms to represent it (see e.g., Porciani & Madau 2001; Daigne et al. 2006; Wanderman & Piran 2010; Salvaterra et al. 2012; Amaral-Rogers et al. 2016).

In this work we chose to use a simple functional form, with only 3 free parameters, which could adequately represent the CSFRD while also having the freedom to deviate from it at high or low redshift independently. The comoving rate density of LGRBs ṅ_GRB is thus parametrised in our model as:

$$\begin{array}{c}\hfill {\dot{n}}_{\mathrm{GRB}}(z)={\dot{n}}_{\mathrm{GRB}}^{0}f(z;a,b,z_m)[\mathrm{y}{r}^{-1}Mp{c}^{-3}],\end{array}$$

with f(z) a broken exponential:

$$\begin{array}{c}\hfill f(z;a,b,z_m)=\{\begin{array}{cc}{e}^{\mathit{az}}\hfill & z<z_m\hfill \\ e^{\mathit{bz}}{e}^{(a-b)z_m}\hfill & z\ge z_m\hfill \end{array},\end{array}$$(4)

where z_m is the redshift of the break, a and b are the low- and high-redshift slopes, respectively, and ${\dot{n}}_{\mathrm{GRB}}^0$ is a normalisation given by our model (see Sect. 3.4).

We can relate ṅ_GRB to the CSFRD by introducing the LGRB efficiency η(z) defined as:

$$\begin{array}{c}\hfill {\dot{n}}_{\mathrm{GRB}}(z)=\eta (z){\dot{n}}_{\mathrm{cc}}(z),\end{array}$$

where ṅ_cc is the comoving rate density of core collapses (in units of [yr⁻¹ Mpc⁻³]). We then assume that ṅ_cc is related to the CSFRD by:

$$\begin{array}{c}\hfill {\dot{n}}_{\mathrm{cc}}(z)=\frac{p_{\mathrm{cc}}(z)}{\overline{m}(z)}{\dot{\rho }}_{\ast }(z),\end{array}$$

where ${\dot{\rho }}_{\ast }(z)$ is the CSFRD in units of [M_⊙ yr⁻¹ Mpc⁻³] and $\overline{m}$ is the mean mass deduced from the stellar IMF:

$$\begin{array}{c}\hfill \overline{m}(z)={\displaystyle {\int }_{m_{\inf }}^{m_{\sup }}mI(m,z)\mathrm{d}m[M_{\odot }].}\end{array}$$

Here, p_cc(z) is the fraction of formed stars ending with a core collapse, given by:

$$\begin{array}{c}\hfill p_{\mathrm{cc}}(z)={\displaystyle {\int }_{m_{\mathrm{cc}}}^{m_{\sup }}I(m,z)\mathrm{d}m,}\end{array}$$

where m_cc is the minimum mass at which a core collapse can occur. In our case we used a Salpeter (Salpeter 1955) stellar IMF with a slope of 1.35 for I(m, z), m_cc = 8 M_⊙, m_inf = 0.1 M_⊙, and m_sup = 100 M_⊙. Using the numeric values quoted previously, this finally leads to:

$$\begin{array}{c}\hfill {\dot{n}}_{\mathrm{GRB}}(z)=\frac{\eta (z){\dot{\rho }}_{\ast }(z)}{135M_{\odot }}·\end{array}$$(5)

In the following, we assume p_cc and $\overline{m}$ to be constant with redshift, meaning any discrepancy between the redshift distribution of LGRBs and the CSFRD is caused by the redshift evolution of the LGRB efficiency η(z). As discussed in Sect. 5.4, such an evolution should be directly linked to the physics of LGRB progenitors.

The CSFRD is then taken to have a similar shape as the broken exponential defined above:

$$\begin{array}{c}\hfill {\dot{\rho }}_{\ast }(z)={\dot{\rho }}_{\ast }^{0}f(z;a_{\ast },b_{\ast },z_{m,\ast })[M_{\odot }\mathrm{y}{r}^{-1}Mp{c}^{-3}],\end{array}$$

with ${\dot{\rho }}_{\ast }^0=2.8\times {10}^{-2}{M}_{\odot }{\mathrm{yr}}^{-1}{\mathrm{Mpc}}^{-3}$, a_* = 1.1, b_* = −0.57 and z_m, * = 1.9; these values are obtained by fitting the broken exponential form of Eq. (4) to the slightly more complex functional form of Springel & Hernquist (2003) with the parameter values of model 3 from Vangioni et al. (2015), which are deduced of the observed cosmic formation rate extracted by Behroozi et al. (2013) from galaxy luminosity functions measurements including high redshift measurements at z ≃ 8 − 10 (Oesch et al. 2014, 2015; Bouwens et al. 2015), and are compatible with the chemical enrichment history of the Universe and with the constraints on reionisation from the cosmic microwave background data (see Vangioni et al. 2015 for details).

The particular case in which η is constant with redshift corresponds to a = a_*, b = b_* and z_m = z_m, *. Then the efficiency η is equal to

$$\begin{array}{c}\hfill {\eta }_0=135M_{\odot }\frac{{\dot{n}}_{\mathrm{GRB}}^0}{{\dot{\rho }}_{\ast }^0}=\frac{{\dot{n}}_{\mathrm{GRB}}^0}{2.1\times {10}^{-4}\mathrm{y}{r}^{-1}Mp{c}^{-3}}·\end{array}$$

In the more general case where a, b, and z_m are different from the CSFRD values a_*, b_* and z_m, *, we can obtain the evolving LGRB efficiency η(z) by rearranging Eq. (5):

$$\begin{array}{c}\hfill \eta (z)={\eta }_0\frac{f(z;a,b,z_m)}{f(z;a_{\ast },b_{\ast },z_{m,\ast })}·\end{array}$$(6)

2.3. Intrinsic luminosity function

The luminosity function is another widely studied distribution of LGRB populations with a variety of different functional forms used (e.g., Salvaterra et al. 2012; Wanderman & Piran 2010; Pescalli et al. 2016). In our model, due to its ubiquitous use throughout astronomy and its flexibility, we described it using a Schechter function (Schechter 1976) as:

$$\begin{array}{c}\hfill \varphi (L)=A\{\begin{array}{cc}{\left(\frac{L}{L_{\ast }}\right)}^{-p}\times exp(-\frac{L}{L_{\ast }})\hfill & L>L_{\min }\hfill \\ 0\hfill & L\le L_{\min }\hfill \end{array},\end{array}$$(7)

where A is a normalisation factor given by $1={\int }_0^{\infty }\varphi (L)\mathrm{d}L$. In practice this means we are actually working with a probability density but we will continue referring to it as a luminosity function out of convenience. This functional form has the advantage of having only three parameters: the minimum luminosity L_min, the break luminosity L_*, and the slope p. It is worth noting that this form is very similar to a simple power law, with the difference that the break at the high-luminosity end is smooth rather than sudden, which is slightly more realistic (see the discussion in Atteia et al. 2017). We also tested luminosity models with a broken power law but found little success in deriving meaningful constraints on all parameters simultaneously². For this reason, and in order to have a lower number of free parameters, we chose to use the Schechter function. Furthermore, after some exploration we decided to fix L_min since it cannot be constrained from current observations as it requires seeing the turnover at low peak fluxes in the log N − log P diagram, which is to date unobserved (Stern et al. 2001). The value chosen was L_min = 5 × 10⁴⁹ erg s⁻¹, which corresponds to the lowest-luminosity burst in the Swift/eBAT6 sample (see Sect. 5.2 for a discussion on the low-luminosity population in our model); this value is similar to the ones assumed by other studies (typically taken between 10⁴⁸ and 10⁵⁰ erg s⁻¹). It should be noted that this parameter strongly affects the normalisation of our model and in particular the global LGRB rate (which we will discuss in Sect. 5); however, it does not affect the observed rate of LGRBs with current instruments (i.e., the number of LGRBs in the various mock samples presented in Sect. 2.6) since the majority of their bursts have luminosities above L_min.

We explored two scenarios for the intrinsic luminosity function: the first in which it is constant, and the second in which it is allowed to evolve with redshift. In the second case we assume for simplicity that the slope of the luminosity function p is constant and only L_min and L_* evolve with redshift as (1 + z)^k_evol as in Salvaterra et al. (2012). The redshift-evolving luminosity function is given by:

$$\begin{array}{c}\hfill \varphi (L,z)=\frac{1}{{(1+z)}^{k_{\mathrm{evol}}}}\varphi (\frac{L}{{(1+z)}^{k_{\mathrm{evol}}}},z=0),\end{array}$$(8)

where the 1/(1 + z)^k_evol pre-factor comes from the condition that the probability density remain normalised to unity. In this case, the values of the parameters quoted are always given for the de-evolved luminosity function, at z = 0.

2.4. Intrinsic photon spectrum

There are a number of different functional forms used to fit GRB high energy spectra (typically 10 keV–1 MeV) throughout the literature (see e.g., Preece et al. 2000; Kaneko et al. 2008; Goldstein et al. 2013; Gruber et al. 2014; Bhat et al. 2016). We decided to use the Band function (Band et al. 1993), which provides the best fits for GRBs with high signal-to-noise ratio in the GBM catalogue (see left panels of Fig. A.1) and whose shape has only three parameters: the low-energy slope α, the high-energy slope β, and the peak energy E_p. The Band function is presented in more detail in Appendix A.

2.4.1. Intrinsic peak energy distribution

In the past, various authors have found relations between the spectral and energetic properties of GRBs and in particular between the peak energy E_p of the GRB spectrum and the isotropic-equivalent energy E_iso (Amati et al. 2002; Amati 2006; Lu et al. 2012) or luminosity L_iso (Yonetoku et al. 2004, 2010; Frontera et al. 2012). Some authors have suggested that these correlations are caused by strong selection effects (Nakar & Piran 2005; Band & Preece 2005; Butler et al. 2007; Shahmoradi & Nemiroff 2011; Heussaff et al. 2013), while others have argued that selection effects do not suffice to explain the observed correlation (Ghirlanda et al. 2008, 2012; Nava et al. 2008). For our models we tested two different scenarios regarding the E_p distribution of LGRBs: a scenario with an intrinsic correlation between E_p and L_iso and a scenario where the intrinsic E_p distribution follows a log-normal distribution independently of L_iso.

Several population models (e.g., Salvaterra et al. 2012; Pescalli et al. 2016; Lan et al. 2019) assume that the observed correlation between the peak luminosity of a burst and its peak energy is intrinsic. The correlation has three free parameters and can be written as:

$$\begin{array}{c}\hfill E_{\mathrm{p}}=E_{p0}{\left(\frac{L}{L_0}\right)}^{{\alpha }_{\mathrm{A}}},\end{array}$$(9)

where E_p0 is drawn from a log-normal distribution with scatter σ_Ep, and L₀ is a constant fixed at 1.6 × 10⁵² erg s⁻¹. The population models cited above use the values obtained by fitting observed samples of LGRBs and assumed to be valid for the intrinsic LGRB population; the most recent observed values from the extended BAT6 sample are E_p0 = 309 ± 6 keV, α_A = 0.54 ± 0.05, σ_Ep = 0.28 (Pescalli et al. 2016). In the intrinsic E_p − L correlation scenario studied in the present paper, we decided to relax this condition and leave the parameters of this relation free to vary to allow for an intrinsic E_p − L correlation to be affected by additional selection effects to shape the observed correlation.

The independent log-normal E_p scenario assumes E_p is independent of L and follows a log-normal distribution with two free parameters: the mean E_p0 [keV] and the spread σ_Ep [dex]. It should be noted that the correlated scenario becomes the independent log-normal scenario when α_A goes to 0. The motivation behind this scenario is to see whether we can obtain an observed E_p − L correlation without imposing any intrinsic correlation.

2.4.2. Intrinsic spectral slopes distribution

The spectral slopes α and β are drawn from the catalogue of Gruber et al. (2014), Bhat et al. (2016). In particular, we chose the sample fulfilling their ‘good’ criteria (i.e., with small errors, Gruber et al. 2014) adding the conditions that 2 < β and α < 2. In this way we have a diversity of realistic values for these parameters as opposed to fixing them as is commonly done in GRB population models (e.g., Wanderman & Piran 2010; Salvaterra et al. 2012).

2.5. Intrinsic duration and variability coefficient distributions

The intrinsic duration and variability coefficients are used solely for crosschecks and are drawn after the best fit parameters for a given scenario have already been found by MCMC (see Sect. 4). This means that all mock samples with a selection based only on the peak flux have already been built, including the Fermi/GBM^bright sample of bright GBM bursts (cut at 0.9 ph s⁻¹ cm⁻² to avoid faint flux incompleteness, see Sect. 3.2). This allows us to compute the probability density function of the intrinsic duration T₉₀ by first fitting the observed distribution of T_90 obs from the GBM catalogue, cutting out LGRBs below 0.9 ph s⁻¹ cm⁻², with a log-normal distribution, which yields a mean μ = 1.45 and standard deviation σ = 0.47. Then we assume the intrinsic distribution of T₉₀ is log-normal with the same spread and correct for cosmological time dilation as $\mathrm{\mu }-log(1+{\overline{z}}_{\mathrm{GBM}})$, where ${\overline{z}}_{\mathrm{GBM}}$ is the median redshift of our mock Fermi/GBM^bright sample.

The variability coefficient C_var is defined as the ratio of the mean luminosity $\overline{L}$ to the peak luminosity L. It can be estimated from the GBM catalogue with:

$$\begin{array}{c}\hfill C_{\mathrm{var}}=\frac{\overline{L}}{L}=\frac{\overline{N}}{N^{\mathrm{pk}}}=\frac{\mathcal{N}}{T_{90\mathrm{obs}}{N}^{\mathrm{pk}}},\end{array}$$(10)

where 𝒩 is the photon fluence (in units of [ph cm⁻²]) and N^pk is the peak photon flux (in units of [ph s⁻¹ cm⁻²]) defined in Eq. (1). The C_var − T_90 obs plane of the GBM catalogue with $N_{50-300\mathrm{keV}}^{\mathrm{pk}}\ge 0.9\mathrm{ph}{\mathrm{s}}^{-1}{\mathrm{cm}}^{-2}$ is shown in Fig. 1 in the case of a Band spectral fit using the flux and fluence values in the 50–300 keV band. We fit the correlation between C_var and T_90 obs with a linear regression and find:

$$\begin{array}{c}\hfill log{C}_{\mathrm{var}}=-0.413log{T}_{90\mathrm{obs}}.\end{array}$$(11)

Fig. 1. Distribution of Cvar and T90 obs for the GBM catalogue of Bhat et al. (2016) for burst with $N_{50-300\mathrm{keV}}^{\mathrm{pk}}\ge 0.9\mathrm{ph}{\mathrm{s}}^{-1}{\mathrm{cm}}^{-2}$. The data are binned in a 2D histogram whose shading represents the density of points; contours are also shown to guide the eye. The black dashed line corresponds to the linear fit of Eq. (11). The black filled lines in the side histograms are the 1D Gaussian kernel density estimations of the data.

We crosschecked these values on different catalogues (BATSE 5B, Goldstein et al. 2013) and different energy bands (GBM 10–1000 keV, Bhat et al. 2016) and find that the slope is within 10% around −0.4. We fit the decorrelated C_var with a log-normal distribution yielding μ = 0.04 and σ = 0.22. Thus to obtain C_var, we: (i) draw z and T₉₀ as mentioned above and calculate T_90 obs, (ii) draw log $C_{\mathrm{var}}^{\mathrm{decorr}}$ from a log-normal distribution with μ = 0.04 and σ = 0.22, and (iii) calculate C_var from Eq. (11). We note that there can be some nonphysical values of C_var > 1, which we set to 1 in our implementation.

2.6. Mock samples

From the distributions presented above, and for a given set of parameters in the case of the adjusted distributions, we can randomly draw and calculate the properties and quantities given in Table 1 for a large number of LGRBs N_GRB. From this intrinsic population we can create several mock samples by applying different peak flux, or fluence cuts; these are summarised in Table 2. The first three mock samples of Table 2 are only based on a peak flux cut and therefore do not depend on T₉₀ and C_var. They are the ones used for comparison with the observational constraints presented in Sect. 3 in order to constrain the parameters of the population’s distributions. The SHOALS sample relies on a fluence criteria and thus also requires T₉₀ and C_var; it is used as a crosscheck (see Sect. 3.5.2) rather than a constraint for this reason.

3. Observational constraints

In order to constrain our LGRB population model, we used carefully selected observational constraints covering all key aspects of the observed LGRB population. We used different missions and instruments to achieve this, taking advantage of the wide variety of observational data publicly available as of October 2018. This use of samples from three different instruments is one new aspect of our population model with respect to others and requires that special attention be paid to the selection process. Additionally, this implicitly assumes that the LGRBs from these various missions all come from the same underlying LGRB population, an assumption that is validated given the good quality fits described in Sect. 4.

3.1. Intensity constraint

One of the most important constraints, in particular for the luminosity function of LGRBs, is one based on the intensity of the bursts. In this instance, the intensity of the bursts can be assimilated to the peak flux, thus it becomes of interest to constrain the number of bursts at each peak flux, traditionally known as the log N − log P diagram³, represented in panel a of Fig. 2.

Fig. 2. Observational constraints used in our population model. Panel a: intensity constraint; log N − log P diagram built from the offline search of CGRO/BATSE LGRBs of Stern et al. (2001), corrected for efficiency of detection at low fluxes. Panel b: spectral constraint; Ep number density distribution from the GBM spectral catalogue (Gruber et al. 2014) for long GRBs with $N_{50-300\mathrm{keV}}^{\mathrm{pk}}>0.9\mathrm{ph}{\mathrm{s}}^{-1}{\mathrm{cm}}^{-2}$. The actual histogram used in the fitting procedure is shown as black crosses while a finer-grained histogram is shown in light grey to guide the eye. The black line is a Gaussian kernel density estimation of the data. Panel c: redshift constraint; redshift number density distribution of the Swift/eBAT6 sample. The black line is a Gaussian kernel density estimation of the data.

This diagram is a good way to estimate the peak isotropic-equivalent luminosity function of GRBs; however, there is a difficulty residing in the fact that while peak fluxes are proportional to the luminosity, they also depend on redshift. This means a burst with a high peak flux could be low luminosity at low redshift, or high luminosity at high redshift. Fruitful studies (Kommers et al. 2000; Stern et al. 2001) have focused on the turnover at low peak flux, trying to determine if it is real (i.e., due to a minimum luminosity of GRBs) or if it is caused by the lower efficiency of detectors at these fluxes. In the case of Stern et al. (2001), who went down to the lowest peak flux limit by conducting an off-line search of all CGRO/BATSE records, they conclude that there is no turnover down to a peak flux of about ∼0.1 ph s⁻¹ cm⁻². Using the catalogue⁴ from Stern et al. (2001), we reconstructed a modified version of their original log N − log P diagram, using wider bins towards high peak fluxes to insure at least ten objects in each bin for proper Gaussian statistics to be applicable.

3.2. Spectrum constraint

In order to constrain the spectral properties of our intrinsic population we focused on a quantity that is fundamental in defining the GRB prompt spectrum: the peak energy, E_p. Similarly to the log N − log P distribution, the E_p obs distribution is the result of the intrinsic E_p distribution, convolved with redshift, which raises the same aforementioned problems. In addition to this, there is the issue of properly measuring E_p, which is difficult for instruments with a narrow energy band (e.g., Swift/BAT 15–150 keV). We decided to use data from the third Fermi/GBM catalogue⁵ (Gruber et al. 2014; von Kienlin et al. 2014; Bhat et al. 2016) since GBM is an instrument with a large sample (≥1500 GRBs) and a broad spectral bandwidth (10 keV–30 MeV). In order to have a clean comparison sample we created the log N − log P diagram for the GBM sample and compared it to the one from Stern et al. (2001), shown in Fig. 3, using the GBM peak flux in the 50–300 keV band. Since a thorough study on the efficiency of the GBM detectors is not yet available, we normalised the log N − log P from GBM to the corrected one from Stern et al. (2001) by multiplying the GBM histogram by a constant and performing χ² minimisation on the bright bins (shown within the grey shaded area in Fig. 3) to find the best value. This normalisation yields an average exposure factor of f_exp ≃ 0.57 for the GBM sample (fraction of observing time multiplied by the fraction of sky observed, see Eq. (13)) slightly higher than the value for CGRO/BATSE of 0.46 (Salvaterra et al. 2012; Goldstein et al. 2013). We then cut the GBM sample at 0.9 ph s⁻¹ cm⁻², indicated by the black dashed vertical line in Fig. 3, below which the log N − log P diagrams of GBM and Stern et al. (2001) start to diverge. By doing this procedure we are ensuring an E_p obs distribution that is unbiased from faint flux incompleteness, at the price of sample size. The resulting E_p obs distribution is shown in panel b of Fig. 2.

Fig. 3. Efficiency-corrected log N − log P diagram for CGRO/BATSE from Stern et al. (2001) in black. The log N − log P of GBM is shown in green, adjusted to the one of CGRO/BATSE by multiplying by a constant whose value is obtained by minimising the χ2 between the two histograms over the grey shaded region. The black dashed vertical line represents the Npk cut to avoid biases due to faint flux incompleteness for our spectral constraint.

3.3. Redshift constraint

As illustrated by the previous sections, the distance of GRBs is inherently intertwined with any observable and unfortunately the majority of GRBs do not have a measured redshift. To date, there are over 500⁶ GRBs with measured redshift; however, it is not possible to simply use all GRBs with a redshift since redshift distributions are often plagued with strong selection effects and biases that are difficult to model (see, however, Wanderman & Piran 2010). For instance, the ability to measure a redshift for GRBs relies fundamentally on the capacity to locate it, which biases this distribution against so-called ‘dark’ bursts (Greiner et al. 2011; Melandri et al. 2012) that exhibit highly extinguished optical afterglows. Another selection effect, called the redshift desert, is due to the fact that most emission and common absorption lines are shifted outside the window of optical spectrographs around z ∼ 2, although the advent of newer spectrographs such as X-shooter (Vernet et al. 2011) has mostly remedied this. It is thus crucial to use a well-controlled redshift distribution to avoid biasing the intrinsic LGRB population, even at the cost of sample size; we therefore used the redshift distribution from the extended BAT6 sample.

The Swift/BAT6 sample (Salvaterra et al. 2012) is a complete sample of Swift LGRBs with a selection based on the peak flux $N_{15-150\mathrm{keV}}^{\mathrm{pk}}>2.6\mathrm{ph}{\mathrm{s}}^{-1}{\mathrm{cm}}^{-2}$ and favourable observing conditions (Jakobsson et al. 2006). These conditions are chosen so that they increase the chance of redshift recovery without biasing the redshift distribution. They are based on criteria that do not depend on the redshift of the LGRB: (i) the burst must be well localised by Swift/XRT and the information was distributed quickly; (ii) there is low galactic foreground extinction (A_V < 0.5); (iii) the burst declination is between −70° and +70°; (iv) the burst’s angular distance to the sun is greater that 55°; (v) there are no nearby bright stars. This results in 58 LGRBs for the original BAT6, later extended to 99 LGRBs by Pescalli et al. (2016). This extended BAT6 sample (Swift/eBAT6) has 82 bursts with redshifts, yielding an 83% redshift completeness, and its redshift distribution is shown in panel c of Fig. 2. It should be noted that while the eBAT6 bursts are statistically representative of the population of Swift bursts with $N_{15-150\mathrm{keV}}^{\mathrm{pk}}>2.6\mathrm{ph}{\mathrm{s}}^{-1}{\mathrm{cm}}^{-2}$, this population constitutes only the brightest 25% of Swift bursts.

3.4. Event rates for our mock samples

One of the advantages of using the CGRO/BATSE log N − log P of Stern et al. (2001) is that they derived a detection efficiency for all bursts within the field of view. Applying this efficiency correction while also taking into account the mean observed solid angle of the sky ⟨Ω^BATSE⟩ = 8.17 sr and the live time of the search $T_{\mathrm{live}}^{\mathrm{BATSE}}=6.54\mathrm{yr}$ (Goldstein et al. 2013), we can use this corrected log N − log P to obtain the total observed all-sky LGRB rate:

$$\begin{array}{c}\hfill R_{\mathrm{BATSE}}^{\mathrm{tot}}=1170(\pm 40)\mathrm{LGRB}{\mathrm{yr}}^{-1}\mathrm{in}4\pi ,\end{array}$$

for BATSE LGRBs (i.e., brighter than 0.067 ph s⁻¹ cm⁻² in the 50–300 keV range). We can use $R_{\mathrm{BATSE}}^{\mathrm{tot}}$ to normalise our population with:

$$\begin{array}{c}\hfill T_{\mathrm{sim}}=\frac{N_{\mathrm{BATSE}}}{R_{\mathrm{BATSE}}^{\mathrm{tot}}}[\mathrm{y}r],\end{array}$$

where N_BATSE is the number of LGRBs in our mock CGRO/BATSE sample and T_sim represents the time it would take to generate N_BATSE LGRBs with an LGRB rate of $R_{\mathrm{BATSE}}^{\mathrm{tot}}$. Defining the intrinsic rate of all LGRBs above L_min as:

$$\begin{array}{c}\hfill R_{\mathrm{intr}}=\frac{N_{\mathrm{GRB}}}{T_{\mathrm{sim}}}[{\mathrm{yr}}^{-1}],\end{array}$$

where N_GRB is the number of LGRBs in our simulated population, we can obtain ${\dot{n}}_{\mathrm{GRB}}^0$ with:

$$\begin{array}{c}\hfill R_{\mathrm{intr}}={\dot{n}}_{\mathrm{GRB}}^0{\displaystyle {\int }_0^{z_{\max }}\frac{f(z;a,b,z_m)}{1+z}\frac{\mathrm{d}V}{\mathrm{d}z}\mathrm{d}z}\end{array}$$(12)

where z_max is a maximum redshift, set to 20 in our implementation. We can also calculate the total rate of LGRBs predicted by our population for a given sample:

$$\begin{array}{c}\hfill R_{\mathrm{sample}}^{\mathrm{tot}}=\frac{N_{\mathrm{sample}}}{N_{\mathrm{GRB}}}{R}_{\mathrm{intr}}[{\mathrm{yr}}^{-1}].\end{array}$$

The observed rate of LGRBs for a given sample is then simply the total rate times the average exposure factor, f_exp:

$$\begin{array}{cc}\hfill R_{\mathrm{sample}}^{\mathrm{obs}}& =f_{\exp }{R}_{\mathrm{sample}}^{\mathrm{tot}}[{\mathrm{yr}}^{-1}]\hfill \\ \hfill & =\frac{⟨\mathrm{\Omega }⟩}{4\pi }\frac{T_{\mathrm{live}}}{T_{\mathrm{total}}}{R}_{\mathrm{sample}}^{\mathrm{tot}},\hfill \end{array}$$(13)

where ⟨Ω⟩ is the average angle of the sky observed by the instrument and T_live and T_total are the live time and total time of observation of the instrument, respectively.

3.5. Additional crosschecks

In addition to the aforementioned constraints used for fitting, we control that the distributions of the spectral slopes α and β of the GBM sample generated by our model are consistent with the observed distributions and include two other observables to help discriminate scenarios with similar likelihood. A brief description is given below, and more details are given in Sects. 5.1 and 5.3.

3.5.1. The Swift/eBAT6 E_p − L plane

The Swift/eBAT6 sample is used as a redshift constraint, but there is more information to be extracted from this complete sample. More specifically, we use the E_p − L plane since it contains information about the observed correlation between the isotropic-equivalent luminosity of the bursts and their peak energy. This is relevant, in particular when trying to distinguish between scenarios with intrinsic correlated or independent log-normal E_p distributions (see Sect. 3.2). Figure 4 shows the E_p − L plane for the Swift/eBAT6 sample (data are from Pescalli et al. 2016). Fitting the E_p − L plane of the extended BAT6 sample, Pescalli et al. (2016) derived E_p0 = 309 ± 6 keV, α_A = 0.54 ± 0.05, σ_Ep = 0.28, their fit is represented by the red line in Fig. 4. It should be noted that ∼25% of the original BAT6 sample have peak energies only determined by Swift/BAT, which can be inaccurate for determining spectral parameters due to its small bandwidth (15–150 keV). This could potentially impact the E_p distribution of the sample, although a precise quantification of this effect has not yet been determined (but see Nava et al. 2012).

Fig. 4. Ep − L plane for the Swift/eBAT6 sample. The individual data points are colour-coded by redshift; the filled red line represents the Ep − Liso relation found by Pescalli et al. (2016) for this sample and the dashed red lines represents the scatter. A 2D Gaussian kernel density estimation of the data is shown as grey shaded contours. The black dashed, dot-dashed, and dotted lines represent the detection threshold for $N_{15-150\mathrm{keV}}^{\mathrm{pk}}=2.6\mathrm{ph}{\mathrm{s}}^{-1}{\mathrm{cm}}^{-2}$ and a fixed Band spectrum (α = 0.6, β = 2.5) at redshifts 0.3, 2, and 5, respectively. The side histograms represent the binned data and the colour of each bin represents the median redshift in that bin following the colour-code of the central panel; the black curve is the 1D Gaussian kernel density estimation.

3.5.2. The SHOALS redshift distribution

The Swift GRB Host Galaxy Legacy Survey (SHOALS, Perley et al. 2016a,b) is the largest unbiased sample of LGRB host galaxies to date with 119 objects and a 92% redshift completeness. The selection is based on a fluence cut, a quantity that is not straightforwardly calculated by our model since we do not simulate light curves for our LGRBs. It can however be calculated from the two additional properties T₉₀ and C_var as shown in Eq. (15). Because of the additional assumptions needed to create the distributions of the duration and the variability indicator, we decided to use only Swift/eBAT6 as our primary redshift constraint and to perform a crosscheck a posteriori to make sure the best fit populations also adequately represent SHOALS.

4. Results

In order to find the best fit parameters, we used a Bayesian framework with an indirect likelihood and a Markov chain Monte Carlo (MCMC) approach to explore the parameter space; the details are presented in Appendix B. Due to computational considerations we tried to reduce the number of free parameters as much as possible in a justifiable way. For instance, in the independent log-normal E_p scenario, we noticed the value of the standard deviation σ_Ep did not move significantly from 0.45, we therefore decided to fix it at this value. Moreover, we found a strong degeneracy between the luminosity evolution parameter k_evol and the slopes of the redshift distribution; we thus fixed k_evol to 0, 0.5, 1.0, and 2.0 and left the redshift distribution free to vary. In the following, these four cases will be referred to as ‘no evolution’, ‘mild evolution’, ‘evolution’ and ‘strong evolution’ of the luminosity function, respectively. The full set of parameters explored and the ranges used for their flat priors are shown in Table B.1.

The first result is that for a population with a non-evolving luminosity function (k_evol = 0) and a constant LGRB efficiency (i.e., a redshift distribution that follows the shape of the CSFRD), we cannot find a set of parameters that yield a good fit to the data. Indeed, as shown in red in Fig. 5 for the intrinsic E_p − L correlation scenario, a population with no evolution of the luminosity function and a constant LGRB efficiency over-predicts the number of bright and low-redshift LGRBs. This is an important result connected to the physics of LGRB progenitors and was already found by Daigne et al. (2006) after two years of Swift observations. It was then confirmed by following studies (Wanderman & Piran 2010; Salvaterra et al. 2012) and we confirm it again on larger, more precise samples and find that it is true both in the independent log-normal E_p scenario and the scenario with intrinsic E_p − L correlations. Our other models can then be used to better understand the underlying evolution (efficiency and/or luminosity) responsible for this discrepancy.

Fig. 5. Best fits to the observational constraints from models with intrinsic Ep − L correlations and a luminosity function that does not evolve with redshift (kevol = 0). The red curve corresponds to the case where the LGRB rate follows the shape of CSFRD and the blue curve corresponds to the case where the LGRB efficiency is free to evolve with redshift. The black curves and grey histograms are the observational constraints described in Fig. 2. This figure is discussed in Sect. 4.

Building in complexity and allowing either the redshift distribution parameters or k_evol to vary, good fits to the data are found; the best fit parameters for each case are reported in Table 3 for the independent log-normal E_p scenario and in Table 4 for the intrinsic E_p − L correlation scenario; the quality of the fits in each case is comparable.

The scenarios where the redshift distribution is fixed proportional to the CSFRD find strong evolution of the luminosity function (k_evol = 1.6 ± 0.1 for the independent log-normal E_p and k_evol = 2.1 ± 0.1 for the intrinsic E_p − L scenario). The value of k_evol = 2.1 ± 0.1 is in agreement with previous studies, which find k_evol = 2.1 ± 0.6 (Salvaterra et al. 2012) and k_evol = 2.5 (Pescalli et al. 2016) using similar E_p scenarios.

Results obtained with an LGRB rate strictly proportional to the cosmic star formation rate are consistent with results obtained by fixing k_evol and leaving the LGRB rate free to vary. In the intrinsic E_p − L case, this corresponds to a strong cosmic evolution of the luminosity function (k_evol = 2) and in the independent log-normal E_p case, this corresponds to an intermediary case between k_evol = 1 and k_evol = 2. For this reason in the following, the discussion will focus on the scenarios where the LGRB rate is left free to vary.

The intrinsic E_p − L correlation scenarios find an intrinsic slope of the correlation significantly shallower than observed (α_A ≃ 0.3 versus ∼0.5 observed) and a larger dispersion (σ_Ep ≃ 0.45 versus ∼0.28), close to the best value derived in the independent log-normal E_p scenario. The value of these best parameters are independent of the value of the scenario for the cosmic evolution of the luminosity and/or the rate.

The behaviour of the parameters of the luminosity function and the redshift distribution is similar in both scenarios; as k_evol increases, L_* and the slopes of the redshift distribution, a and b, decrease. Comparing the two scenarios for the peak energy E_p, we find very similar results regarding the luminosity function for the same value of k_evol: p differs by less than 0.1, L_* is typically bigger by a factor of ∼2 in the intrinsic E_p − L correlation scenario. We also find similar results regarding the redshift distribution for the same value of k_evol: the break redshift z_m ≃ 2.1 − 2.2 independently of k_evol and of the E_p scenario, the high-redshift slope b varies strongly with k_evol and is similar in both E_p scenarios and finally the low-redshift slope a also varies strongly with k_evol but is typically steeper by +0.2 in the intrinsic E_p − L correlation scenario.

Finally, the results are comparable for both scenarios regarding the peak energy distribution, which is expected given that σ_Ep is the same and that the correlation slope α_A is shallow.

Interestingly, the local LGRB rate ${\dot{n}}_{\mathrm{GRB}}^0$ increases with k_evol in the independent log-normal E_p scenarios while it stays roughly constant in the intrinsic E_p − L correlation scenarios.

We conclude that we cannot distinguish between an independent log-normal E_p scenario and an intrinsic E_p − L correlation scenario solely using the constraints presented in Sect. 3; we discuss this in more detail in Sect. 5.1 using additional crosschecks. Nonetheless, many results do not depend on the E_p scenario and can be robustly discussed equivalently; this is the case in particular regarding the redshift evolution of the luminosity function and/or of the LGRB efficiency. Using only the intensity, spectral and redshift constraints we cannot distinguish between luminosity-evolving and efficiency-evolving scenarios. This degeneracy has been an issue in LGRB population models based on smaller, complete samples of bright LGRBs from Swift (see e.g., Salvaterra et al. 2012; Pescalli et al. 2016). Despite using fainter and larger samples, we find that lifting this degeneracy remains difficult; we discuss this in more detail using additional crosschecks in Sect. 5.3.

5. Discussion

5.1. A possible intrinsic spectral–energetic correlation

The observed E_p − L plane of the extended BAT6 sample is a powerful tool to compare the predictions of the intrinsic E_p − L correlation and the independent log-normal E_p scenarios. There has been significant debate concerning the causes of this observed correlation, whether due to observational selection effects or not (e.g., Nakar & Piran 2005; Band & Preece 2005; Butler & Kocevski 2007; Ghirlanda et al. 2008; Nava et al. 2008; Heussaff et al. 2013). Two caveats should be noted while comparing our population to the observations. The first is that the observed Swift/eBAT6 sample is missing about ∼17% of measurements in the E_p − L plane; the second is that ∼25% of the original BAT6 sample have peak energies only determined by Swift/BAT, which can be inaccurate for determining spectral parameters due to its small bandwidth (15–150 keV). The quantification of these effects on the observed Swift/eBAT6 E_p − L plane is beyond the scope of this paper and we do not perform rigorous statistical tests for this reason; this comparison is meant more as an indication. The observed E_p − L plane is shown for our mock eBAT6 sample in Fig. 6 for our best fit models with the two extreme values k_evol = 0 and 2 (no/strong luminosity evolution) and the two scenarios regarding the peak energy E_p. Looking at this plane for the independent log-normal E_p scenarios with mild or no evolution, we find a dearth of high-L bursts compared to the observations (panel a of Fig. 6); in the case of strong evolution, this discrepancy is less pronounced (panel b of Fig. 6). On the other hand, the luminosity distributions of the intrinsic E_p − L correlation scenarios agree well with the observations for all values of k_evol (panels c and d of Fig. 6). For this reason this scenario is preferred and the following discussions will be focused on this scenario, although we checked that the general trend was the same for the independent log-normal E_p scenario. This suggests that the observed E_p − L correlation cannot be explained only from selection effects, in agreement with Ghirlanda et al. (2008, 2012), Nava et al. (2008).

Fig. 6. Swift/eBAT6 Ep − L plane for models with no redshift evolution of the luminosity function (a, c) and with strong redshift evolution of the luminosity function (b, d). Top panels (a, b): independent log-normal Ep scenarios; bottom panels (c, d): intrinsic Ep − L correlation scenarios. The predictions for the entire intrinsic LGRB population is shown as grey contours and the predictions for the mock eBAT6 sample are shown in coloured filled contours. The observed eBAT6 sample is shown as coloured scatter points, the 2D Gaussian kernel density estimation shown in Fig. 4 is omitted for clarity. A full description of the observed data points, the side histograms and the dashed, dot-dashed and dotted curves of the figure is presented in Fig. 4.

Nonetheless, looking at the best fit values in Table 4, we find that the slope of the intrinsic relation α_A is shallower than the observed relation (∼0.3 intrinsic versus 0.54 observed) and the intrinsic scatter is larger (∼0.44 intrinsic versus 0.28 observed). This implies that some selection effects do shape the observed correlation by making it more pronounced and narrower, as is suggested by the detection limits shown as dashed, dot-dashed, and dotted lines for redshift 0.3, 2, and 5, respectively in Fig. 4. This result is in good agreement with the current understanding of the GRB prompt emission physics. The three main possible emission mechanisms (internal shocks, Rees & Meszaros 1994; Kobayashi et al. 1997; Daigne & Mochkovitch 1998; magnetic reconnection, Usov 1994; Spruit et al. 2001; Giannios & Spruit 2005; Zhang & Yan 2011; McKinney & Uzdensky 2012; dissipative photosphere, Mészáros & Rees 2000; Beloborodov & Mészáros 2017) indeed predict that the peak energy increases with the luminosity, but depends also on several other parameters, usually linked both to the jet properties (e.g., Lorentz factor) and the microphysics of the emitting region. Therefore, theoretical models predict a general tendency for the peak energy to increase with luminosity, but never a strict and narrow E_p − L correlation unless several physical parameters of the model are strongly correlated. This has for instance been studied precisely in the case of the internal shock model by Barraud et al. (2005), Mochkovitch & Nava (2015). Our new result can now better specify which intrinsic property of the prompt GRB emission should be compared to theoretical models: a weak correlation with a large ∼0.4 dex scatter of the form

$$\begin{array}{c}\hfill E_{\mathrm{p}}\simeq (500-600\mathrm{keV}){\left(\frac{L}{1.6\times {10}^{52}\mathrm{e}rg{s}^{-1}}\right)}^{0.3}.\end{array}$$(14)

5.2. The population of soft bursts

It is remarkable that in all scenarios, the intrinsic distribution of LGRBs (shown as grey contours in Fig. 6) is dominated by largely undetected low-luminosity bursts, which are also mostly low-E_p bursts in the preferred scenario with a moderate intrinsic E_p − L correlation, in agreement with previous studies (Daigne et al. 2006). This could be linked to observed sub-classes of GRBs, such as low-luminosity GRBs (LL-GRBs, Liang et al. 2007; Virgili et al. 2009; Stanway et al. 2014), X-ray Rich GRBs or X-ray Flashes (XRR-GRBs and XRFs, respectively, Barraud et al. 2005; Sakamoto et al. 2005, 2008). With the spectral range and sensitivity of current instruments, their observational samples are still quite small (Virgili et al. 2009), making the precise description of their properties challenging. For XRFs, the γ-ray spectra are not very well characterised as the E_p is often below the threshold of detectors and furthermore, the afterglows of XRR/XRFs are rarely observed (Sakamoto et al. 2008).

If the intrinsic spectrum–luminosity correlation is indeed valid, one would expect a link between LL-GRBs and XRR-GRBs/XRFs. However, we note that the observed frequency of these events is larger than the predictions of our model. For instance, Sakamoto et al. (2005) proposed a classification based on the softness S defined as the ratio of energy fluences in the 2–30 keV and 30–400 keV bands. Among bursts detected both by the HETE2/WXM and FREGATE instruments, they find 22% of classical GRBs with S < 0.3, 42% of XRR-GRBs with 0.3 < S < 1 and 36% of XRFs with S < 1. For the corresponding mock sample (see Table 2) in our model, the fraction of XRFs is much lower (∼4%) and even XRR-GRBs are underrepresented (∼30%). This suggests that either the luminosity function of LGRBs has a different slope at low luminosity, which is unconstrained by our model, or that XRFs are not just the low-luminosity tail of the classical LGRB population but rather a full-fledged distinct population of their own. These two scenarios have different implications in terms of emission mechanisms and progenitor physics and distinguishing between them would require modelling the population of XRFs in details. Such a dedicated study remains difficult to this day due to the small size of samples of such soft events often without any afterglow or redshift measurements (but see Jenke et al. 2016; Katsukura et al. 2020). We can hope that the situation will improve in the future thanks to a new generation of satellites with a lower-energy threshold and good localisation capabilities such as SVOM/ECLAIRs (4–120 keV, Godet et al. 2014; Wei et al. 2016), THESEUS (2 keV–20 MeV, Amati et al. 2018), HiZ-GUNDAM (1–10 keV, Yonetoku et al. 2014) or the Gamow Explorer (0.5–10 keV, White 2020).

5.3. Cosmic evolution of LGRBs: rate or luminosity?

When discussing the cosmic evolution of LGRBs and especially the respective contribution of the evolution of the comoving rate or the luminosity function, a key property to address is that of metallicity. Long gamma-ray bursts have been shown to occur in predominantly lower metallicity environments with respect to what would be expected if LGRBs were pure tracers of star formation (Vergani et al. 2015; Japelj et al. 2016; Perley et al. 2016b; Graham & Fruchter 2017; Palmerio et al. 2019). At first order, assuming the effect of metallicity is a simple threshold on the LGRB production efficiency from stars, this could suggest that the cosmic evolution of the LGRB rate scenario is more plausible. One should keep in mind however that the effects of metallicity may be more subtle such as modifying the stellar IMF (La Barbera et al. 2013) or the fraction of binary progenitors (Fryer & Heger 2005; de Mink et al. 2009; Podsiadlowski et al. 2010; Chrimes et al. 2020). The potential dependence of LGRB progenitors properties on metallicity could very well lead not only to an evolution of the stellar efficiency to produce LGRB, η(z), but also to an evolution of the LGRB properties. For instance, cosmic evolution of the isotropic equivalent luminosity may be expected if the beaming angle evolves (Lloyd-Ronning et al. 2020; Salafia et al. 2020) since the beaming depends on properties of the progenitor, such as the stellar density profile, which are in turn affected by metallicity. Therefore, the reality may be an intermediary case where both the rate and the luminosity of LGRBs evolve with cosmic time. Unfortunately, our results presented in Sect. 4 show that, solely with the three redshift, spectral, and intensity constraints used by our population model, it is not possible to distinguish between an efficiency and a luminosity evolution. In the following, we attempt to lift this degeneracy by extending our analysis to a larger unbiased sample (SHOALS) and pushing further our comparison between the two extreme scenarios of k_evol = 0 and k_evol = 2.

Having found a number of good models that fit the three observational constraints described in Sect. 3, we can compare the predictions of each model for the redshift distribution of the SHOALS complete, unbiased sample (Perley et al. 2016a,b), which is ∼20% larger than the eBAT6 sample. This comparison is used as a crosscheck rather than a constraint since it involves an additional step, namely the calculation of a fluence, which introduces extra uncertainties. Calculating fluences is not straightforward in our model because we do not simulate a light curve for each LGRB drawn; instead we draw C_var and T_90 obs following the prescription presented in Sect. 2.5 and, using the peak photon flux N^pk following Eq. (1), we get the photon fluence 𝒩 (in units of [ph cm⁻²] in the same energy interval [E_min, obs;E_max, obs]) with:

$$\begin{array}{c}\hfill \mathcal{N}=C_{\mathrm{var}}{T}_{90}{N}^{\mathrm{pk}}.\end{array}$$(15)

This assumes that the time-integrated prompt spectrum is the same as the peak-brightness prompt spectrum which is a strong additional assumption. However, using bursts from the observed Fermi/GBM^bright sample (with $N_{50-300\mathrm{keV}}^{\mathrm{pk}}\ge 0.9\mathrm{ph}{\mathrm{s}}^{-1}{\mathrm{cm}}^{-2}$) we find that the peak energy of these two spectra are tightly correlated: $⟨log(E_{\mathrm{p}}^{\mathrm{Pk}}/E_{\mathrm{p}}^{\mathrm{Ti}})⟩=0.04\pm 0.2$. Additionally, we compared the peak flux–fluence planes of our mock samples with observations and found good agreement suggesting that our method is sound.

The SHOALS sample selection is based on the energy fluence ℱ (in units of [erg cm⁻²] between E_min, obs and E_max, obs) which can be consequently calculated from:

$$\begin{array}{c}\hfill \mathcal{F}=C_{\mathrm{var}}{T}_{90}{F}^{\mathrm{pk}},\end{array}$$(16)

where F^pk is the peak energy flux (in units of [erg s⁻¹ cm⁻²] in the same energy interval) given by Eq. (3). The observed SHOALS sample is defined as having ℱ_{15−150 keV} > 10⁻⁶ erg cm⁻² and favourable observing conditions (see Perley et al. 2016a for more details). We construct our simulated SHOALS sample using this energy fluence threshold and look at the predictions regarding the redshift distribution. We restrict ourselves to scenarios with k_evol = 0 (no evolution) and k_evol = 2 (strong evolution) as they are the ones that will exhibit the most differences, the resulting distributions are shown in Fig. 7. All populations that fit the observational constraints described in Sect. 3 predict a mock SHOALS redshift distribution that cannot be distinguished from the observed SHOALS at the 95% confidence level as is attested by performing K–S tests between the cumulative distributions (p-values range from 0.7 to 0.95).

Fig. 7. Cumulative redshift distribution of the observed SHOALS sample in black; the 95% confidence bound calculated from bootstrapping is shown in grey and the upper limits are shown as arrows at the bottom of the plot. The predictions for the intrinsic Ep − L correlation scenarios with kevol = 0 and kevol = 2 are shown in blue and dashed light blue, respectively.

This means that even with the sample size and fluence cut of SHOALS (117 objects, ℱ_{15−150 keV} > 10⁻⁶ erg cm⁻²), we cannot distinguish between the scenarios with no redshift evolution of the luminosity function (k_evol = 0) and scenarios with strong evolution (k_evol = 2). We investigated for what values of sample size and fluence cut we could distinguish between the k_evol = 0 and the k_evol = 2 scenarios using KS-tests; the results are shown in Fig. 8 while the detailed methodology is shown in Appendix C. This plot illustrates how difficult it is to lift the degeneracy between the cosmic evolution of the LGRB rate and of the LGRB luminosity. The sample size should be increased by a factor larger than 10 at fixed fluence cut in order to have a good chance to distinguish between the two types of cosmic evolution at the 95% confidence level; decreasing only the fluence cut does not really improve the test. In practice, a combination of lower fluence cut and larger sample size could be used, but even then a factor of ∼3 in both directions would be needed to get a decent probability of discriminating between the two types of cosmic evolution at the 95% confidence level. Lifting the degeneracy between luminosity and efficiency cosmic evolution remains difficult using the redshift distribution of reasonably sized samples.

Fig. 8. Fluence cut–sample size plane for the intrinsic Ep − L correlation scenario. The colour shading corresponds to the probability of discriminating between the kevol = 0 and the kevol = 2 scenarios at the 2σ (95%) confidence level. The star corresponds to the current SHOALS sample of Perley et al. (2016b).

Another approach to lifting the degeneracy is to look directly at the luminosity distribution in different redshift bins. Figure 9 shows the luminosity CDF in four redshift bins for the intrinsic E_p − L correlation scenario at different peak flux cuts. The scenario with k_evol = 0 is shown in dashed while k_evol = 2 is shown in full. As expected, if one had access to the entire population (leftmost panel), it would be easy to distinguish between no luminosity evolution and strong luminosity evolution since the CDF would be the same in all redshift bins for the case k_evol = 0. The difficulty resides in the fact that as we increase the peak flux cut (as shown in middle panel of Fig. 9), we remove mostly lower-luminosity bursts and shift the distributions towards higher luminosities, essentially creating the same effect as luminosity evolution. This is illustrated in the centre left and centre right panels of Fig. 9 where the same plots are shown but for bursts with $N_{15-150\mathrm{keV}}^{\mathrm{pk}}\ge 0.1\mathrm{ph}{\mathrm{s}}^{-1}{\mathrm{cm}}^{-2}$ and 2.6 ph s⁻¹ cm⁻², respectively. The rightmost panel of Fig. 9 shows the same curves as in the centre right panel but adding the 95% confidence area from the observed eBAT6 sample. In addition to the good agreement between our models and the observed data, we can see with current sample sizes (∼100 objects) that it is not possible to argue in favour of a specific value for k_evol. Both of these tests become even more difficult if we include the intermediate scenarios with k_evol = 0.5 and k_evol = 1.

Fig. 9. Luminosity cumulative distribution function in four different redshift bins for the intrinsic Ep − L correlation scenario. The full lines correspond to kevol = 2 and the dashed lines to kevol = 0. Leftmost panel: entire intrinsic LGRB population, centre left and centre right panels: bursts with $N_{15-150\mathrm{keV}}^{\mathrm{pk}}\ge 0.1\mathrm{ph}{\mathrm{s}}^{-1}{\mathrm{cm}}^{-2}$ and 2.6 ph s−1 cm−2 respectively. Rightmost panel: same as the centre right panel but adding the observed eBAT6 sample and where each curve was shifted by a constant C for clarity. The shaded area represents the 95% confidence region for the CDF of the observed eBAT6 sample, calculated using bootstrapping (without taking into account limits).

5.4. The LGRB rate and its connection to the cosmic star formation rate

5.4.1. The LGRB production efficiency by massive stars

The top panel of Fig. 10 shows the LGRB comoving rate density as a function of redshift for the intrinsic LGRB population in the intrinsic E_p − L correlation scenario for different values of k_evol. The case of k_evol = 2 follows closely the star formation rate density, while the high-redshift slope gets shallower as k_evol decreases. Models with lower values of k_evol have higher global rates of LGRBs. Using Eq. (6), we can derive the LGRB efficiency as a function of redshift, shown in the bottom panel of Fig. 10. Assuming p_cc and $\overline{m}$ are constant with z, models with k_evol < 2 imply some amount of evolution of the LGRB efficiency with redshift, with the strongest evolution found for models with k_evol = 0. There are many possible causes for this increasing LGRB efficiency with redshift, one of which is metallicity. This is in line with observational studies of complete, unbiased samples of LGRB host galaxies (Hjorth et al. 2012; Vergani et al. 2015, 2017; Japelj et al. 2016; Perley et al. 2016a; Palmerio et al. 2019) that find that metallicity is the driving factor behind the LGRB production efficiency. The fraction of star-formation happening below a given metallicity (from Eq. (5) of Langer & Norman 2006) is shown in the bottom panel of Fig. 10 as dashed lines, arbitrarily scaled. We can see some similarity between our derived efficiency and these curves, although the behaviour at high redshift is different. This might be an indication of an additional break at z ≥ 6 in the redshift distribution of LGRBs, although with the amount of data available at these redshifts to date, it would be hard to constrain. Another possibility is that the effect of metallicity is more complicated than a simple threshold. Indeed, metallicity can also play an indirect role by influencing the stellar IMF (e.g., La Barbera et al. 2013) or the fraction of binary progenitors (Chrimes et al. 2020).

Fig. 10. Cosmic evolution of the LGRB rate and efficiency. Top panel: comoving LGRB rate density as a function of redshift for the best fitting models of the intrinsic Ep − L correlation scenarios. Bottom panel: LGRB efficiency η(z) as as function of redshift obtained with Eq. (6). The estimation of the fraction of star formation below Zth from Eq. (5) of Langer & Norman (2006) is shown, rescaled for comparison, in dashed. The shaded area represents the 95% confidence interval.

5.4.2. The local LGRB rate

The local LGRB rate density for LGRBs pointing towards us we derived with our population model is ${\dot{n}}_{\mathrm{GRB}}^0=0.7-0.8{\mathrm{yr}}^{-1}{\mathrm{Gpc}}^{-3}$ for the intrinsic E_p − L correlation scenarios and ${\dot{n}}_{\mathrm{GRB}}^0=1.0-1.5{\mathrm{yr}}^{-1}{\mathrm{Gpc}}^{-3}$ for the independent log-normal E_p scenarios (see Tables 3 and 4 for details). This local LGRB rate density is highly dependent on the value of L_min = 5 × 10⁴⁹ erg s⁻¹; the following equation can be used to extrapolate to lower values of L_min for a Schechter or power-law function assuming that the additional LGRBs are not detected in the observed samples:

$$\begin{array}{c}\hfill {\dot{n}}_{\mathrm{GRB}}^0(L_{\min })={\dot{n}}_{\mathrm{GRB}}^0(L_{\min }^{\mathrm{ref}}){\left(\frac{L_{\min }}{L_{\min }^{\mathrm{ref}}}\right)}^{1-p},\end{array}$$(17)

where p is the slope of the luminosity function and $L_{\min }^{\mathrm{ref}}=5\times {10}^{49}$ erg s⁻¹ in our case. Taking into account the difference in L_min, our values are mostly consistent with those of previous authors (but see Appendix D). We should note however that the slope at low-luminosities may be different (see discussion in Sect. 5.2) and therefore any extrapolation to lower L_min should keep this in mind.

Using Eq. (5) we find η₀ = (3−6) × 10⁻⁶, as shown in the bottom panel of Fig. 10, or in other words 3 − 6 LGRBs pointing towards us for every million core collapse. This can be extended by dividing by the average jet opening angle to estimate the total number of LGRBs occurring per core collapse (including LGRBs pointing away from us). Using ⟨θ_jet⟩≃5° (Harrison et al. 1999; Ghirlanda et al. 2013; Gehrels & Razzaque 2013; Lloyd-Ronning et al. 2020), we get 0.8 − 1.6 LGRBs pointing in any direction per thousand core collapse at z = 0 and up to a factor of 10 higher at z = 6 (depending on the value of k_evol). This highlights the fact that LGRBs are rare events among transients, requiring specific conditions to form and that these conditions may be more readily met earlier in the Universe.

5.4.3. Measuring the CSFRD with LGRBs

Due to their association with massive stars, LGRBs have been used as probes of the CSFRD (Kistler et al. 2008; Robertson & Ellis 2012; Hao et al. 2020), in particular up to high redshift where estimations from rest-frame UV measurements are plagued by dust uncertainties (Bouwens et al. 2015, 2016). However, in order to estimate the CSFRD from the LGRB rate, one must make hypotheses on the redshift evolution (or lack thereof) of the LGRB luminosity function and the LGRB production efficiency. Typically what is done is to calibrate the relationship between the CSFRD and the rate of bright LGRBs (by using only LGRBs above a certain luminosity) at low redshift (z < 4) and extrapolate this relationship out to high redshift. In doing so, authors usually assume a non-evolving luminosity function; however, as we have shown in Sect. 5.3, there is a strong degeneracy between the cosmic evolution of the LGRB efficiency and the LGRB luminosity function, which cannot be lifted with the current sample sizes. This means that estimates of the CSFRD using LGRBs should take into account the uncertainty due to the possible redshift evolution of the luminosity function and reflect this in their error bars. This correction is not a small effect, looking at Fig. 10, we see that at z = 6, it amounts to a factor of at least 10 between the extreme scenarios.

5.5. Strategy for high redshift detection

Using our intrinsic LGRB populations, we can investigate which type of strategies are most effective at observing high-redshift LGRBs, a primary goal for using the full potential of LGRBs as a cosmological tool. In Fig. 11 we show the expected all-sky yearly rate of LGRBs at z ≥ 6 as a function of limiting peak flux for different energy bands. We see that, at a given peak flux limit, the rate of high-redshift LGRBs is higher for softer energy bands, and depends in particular on the lower-energy threshold of the band, in line with the results of Ghirlanda et al. (2015). This bodes well for SVOM, to be launched in 2022 (Wei et al. 2016), whose coded-mask telescope ECLAIRs is designed to trigger at these lower energies with a low-energy threshold of 4 keV (Godet et al. 2014). This holds even more for other future missions, such as THESEUS⁷, which is in discussion for the next decade and whose XGIS instrument will detect GRBs in the 2 keV–20 MeV energy band with a larger effective area (Amati et al. 2018), the Gamow Explorer (White 2020) which would detect GRBs with a lower-energy threshold (potentially as low as 0.5 keV), or the HiZ-GUNDAM (Yonetoku et al. 2014), with its 1–10 keV wide field X-ray imaging detector. Furthermore, we see that for models with k_evol = 2, the rate of LGRBs at z ≥ 6 reaches a plateau at a peak flux limit N^pk ∼ 10⁻² ph s⁻¹ cm⁻²; this implies that the whole of the population is seen and going to lower peak fluxes will not increase the proportion of high-redshift LGRBs. This could be an interesting test to distinguish between different luminosity evolution scenarios. However, the possibility of getting unbiased, redshift-complete samples of LGRBs down to such low peak fluxes remains to this day unlikely.

Fig. 11. All-sky rate of LGRBs at z ≥ 6, as a function of limiting peak flux for different energy bands in the intrinsic Ep − L correlation scenario. Each panel corresponds to a luminosity evolution scenario (from no evolution at kevol = 0 to strong evolution at kevol = 2). To obtain the actual observed rate for a given mission, one should take into account the field of view of the instrument and its duty cycle. The shaded area represents the 95% confidence interval.

6. Conclusion

We created a model for the intrinsic population of LGRBs, which we constrained by comparing to carefully selected observational constraints, specifically tailored to constrain all major aspects of the population. We explored two main scenarios concerning the distribution of the peak energy E_p: one where the peak energy follows a log-normal distribution independently of other properties (independent log-normal E_p) and one where the peak energy is correlated to the luminosity (intrinsic E_p − L correlation). In addition, in each case we explored four scenarios for the evolution of the luminosity function: no evolution (k_evol = 0), mild evolution (k_evol = 0.5), evolution (k_evol = 1) and strong evolution (k_evol = 2). We derived the best fit parameters in each scenario using a Bayesian Markov chain Monte Carlo exploration scheme with an indirect likelihood. Using two additional crosschecks we discussed the results of our intrinsic population in the context of the intrinsic spectral-energetics correlations, soft bursts, the cosmic evolution of the luminosity function, the LGRB rate and its connection to the cosmic star formation rate and the strategy for detection of high redshift bursts. Our conclusions can be summarised as follows:

• The observed E_p − L correlation cannot be explained only by selection effects; however, these do play a role in shaping the relation. Compared to the observed relation, we derive a shallower slope for the intrinsic correlation α_A ∼ 0.3 and a larger scatter of ∼0.4 dex.

• Our intrinsic population is not able to accurately describe the proportion of soft bursts (e.g., XRR-GRBs, XRFs), which are expected to be weak following the spectral-energetics correlation; this may imply a different slope of the luminosity function at low luminosities or a separate population all-together, either case is difficult to constrain with current samples.

• We confirm the degeneracy between cosmic evolution of the LGRB rate and of the luminosity function. Under a strong cosmic evolution of the luminosity function (k_evol = 2) we find that the LGRB rate follows the cosmic star formation rate; however under a non-evolving luminosity function we find a increasingly higher LGRB production efficiency from stars at high redshift. Naturally, mixed scenarios where both the efficiency and the luminosity function evolve can also reproduce as well our observational constraints.

• We show that this degeneracy cannot be lifted with current sample sizes and completeness. Compared to the current SHOALS sample, a factor of 3 larger sample size and lower fluence cut would be needed to distinguish between the two most extreme scenarios.

• Assuming that the cosmic star formation rate is well measured at high redshift, the scenario with a non-evolving luminosity function (k_evol = 0) implies an LGRB production efficiency from stars a factor of ∼10 higher at z = 6 than at z = 0; assuming it is not, this could be instead interpreted as evidence for the underestimation of the cosmic star formation rate at high redshift.

• The degeneracy between luminosity and redshift evolution is usually not taken into account in studies aiming at measuring high redshift star formation rate using LGRBs, we emphasise the need to do so as it impacts the estimation of the LGRB rate by up to a factor of 10 at z = 6 and 50 at z = 9.

• We derive a local LGRB rate density for LGRBs pointing towards us ${\dot{n}}_{\mathrm{GRB}}^0\sim 0.8{\mathrm{yr}}^{-1}{\mathrm{Gpc}}^{-3}$ in the preferred intrinsic E_p − L scenario which translates to ∼1 LGRB pointing in any direction per thousand core collapse at z = 0, assuming a mean jet opening angle of ⟨θ_jet⟩≃5°.

• We confirm that the best strategy for detecting high redshift bursts is to use detectors with a soft low-energy threshold of ∼ a few keV such as the ECLAIRs telescope on board the SVOM satellite to be launched in 2022 or the XGIS instrument of the future THESEUS mission. We find that the combination of such a low-energy threshold with a good sensitivity could help to partially lift the degeneracy between luminosity and rate evolution.

Acknowledgments

This work has been partially supported by the French Space Agency (CNES). JTP would like to thank Sébastien Carassou, Tom Charnock and Jens-Kristian Krogager for invaluable discussion regarding the statistical framework used in this paper. The corner plots used in this paper for representing multi-dimensional datasets were performed with the PYTHON corner.py code (Foreman-Mackey 2016). Parts of the results in this work make use of the colormaps in the CMasher package (van der Velden 2020).

References

Appendix A: Band function

The Band function (Band et al. 1993) is given by:

$$\begin{array}{c}\hfill f(E)=A\{\begin{array}{cc}{\left(\frac{E}{100\mathrm{keV}}\right)}^{\alpha }exp(-(\alpha +2)\frac{E}{E_{\mathrm{p}}})\hfill & E\le E_c,\hfill \\ {\left(\frac{E}{100\mathrm{keV}}\right)}^{\beta }\left(\frac{E_c}{100\mathrm{keV}}\right)exp(\beta -\alpha )\hfill & E>E_c,\hfill \end{array}\end{array}$$(A.1)

where α is the low-energy slope, β the high-energy slope, E_p the peak of the E L_E spectrum, $E_c=\frac{\alpha -\beta }{\alpha +2}{E}_{\mathrm{p}}$ and A is a normalisation that satisfies $1={\int }_0^{\infty }f(E)\mathrm{d}E$. The goodness of fit of this functional form, shown in the right panel of Fig. A.1 for the 3rd GBM catalogue above 0.9 ph s⁻¹ cm⁻², illustrates that even when Band is not the best fitting model it still provides a good fit to the data (i.e., the reduced χ² is around 1).

Fig. A.1. Quality of Band fits to GRB spectra from GBM catalogue (Bhat et al. 2016). Left: spectral slopes α (top) and β (bottom) as a function of peak flux for the GBM spectral catalogue (Bhat et al. 2016). The entire catalogue of long GRBs above 0.9 ph s−1 cm−2 is in light-grey, the bursts complying with the ‘good’ criteria (i.e., small errors on the parameters) are shown in blue, the bursts for which Band is the best-fit spectral model are shown in orange. Right: reduced χ2 distribution for Band models from the GBM spectral catalogue (Bhat et al. 2016).

Appendix B: Statistical framework

B.1. Likelihood

We decided to follow a methodology already used for galaxy catalogues from image extraction (Carassou et al. 2017), which is based on a Parametric Bayesian Indirect Likelihood (pBIL, Drovandi et al. 2015). The basic idea is to use an auxiliary likelihood when a complex problem renders its own likelihood intractable. To this end, we used the binned maximum likelihood method (Barlow & Beeston 1993), which assumes the number of objects in each bin follows a Poissonian law and we construct our auxiliary likelihood such that for each bin i, the probability of o_i given the model s_i is:

$$\begin{array}{c}\hfill l_i=\frac{e^{-s_i}{s}_i^{o_i}}{o_i!}.\end{array}$$(B.1)

The likelihood for the entire histogram becomes:

$$\begin{array}{c}\hfill \mathcal{L}={\displaystyle \prod _{i=1}^b\frac{e^{-s_i}{s}_i^{o_i}}{o_i!},}\end{array}$$(B.2)

where b is the total number of bins in the histogram and o_i and s_i are the observed and simulated number count in bin i, respectively. The log-likelihood is thus:

$$\begin{array}{c}\hfill ln\mathcal{L}={\displaystyle \sum _{i=1}^b{o}_{i}ln(s_i)-s_i,}\end{array}$$(B.3)

where the factor ln(o_i!) is neglected since it is a constant and our goal is to maximize lnℒ. This likelihood presents a problem if a single simulated bin is empty therefore we added an infinitesimal ϵ = 10⁻³. We checked the impact of different values of ϵ, which affected mostly models with bad likelihood (i.e., models with many empty bins); the effect on good models was negligible. Multiple histograms can be included by simply adding the logarithm of the likelihood (i.e., multiplying the likelihoods) however, the importance of a given histogram on the overall likelihood relies on its number of objects. For this reason, the intensity constraint (see Sect. 3.1) has the strongest weight of our constraints due to its large sample size (N ∼ 7000 once the efficiency correction is taken into account). In order to strengthen the impact of the redshift constraint (weak due to its small sample size, N = 82), we added a weight of 10 to its log-likelihood. We tested different weights and this one was chosen as a balance between having a notable impact and being unrealistically constraining with respect to the other constraints.

B.2. Parameter space exploration

We decided to use a MCMC scheme to explore the parameter space, with a modified version of the Metropolis–Hastings algorithm called simulated annealing (Kirkpatrick et al. 1983). This modification is based on the idea of cooling metals and allows for the Markov chain to initally accept worse likelihoods with a factor τ = τ₀ that decreases at a user-defined rate. As the effective temperature, τ, decreases, the chain is less and less likely to accept worse jumps until τ → 1, where we return to a classic Metropolis–Hastings algorithm. This modification to the original Metropolis–Hastings algorithm verifies the condition of ergodicity: regardless of the starting point the Markov chain will converge to the same stationary distribution.

We also tested the methodological soundness of our algorithm by applying it to data generated by known inputs. We created mock observations for the intensity, spectrum, and redshift constraint which we then used as constraints for our MCMC scheme; the results of the parameter space exploration are presented in Fig. B.1. We can see the code is able to satisfactorily recover the input parameter values: Schechter luminosity function {L_* = 10⁵³ erg s⁻¹, p = 1.5}, redshift distribution following the shape of the CSFRD; intrinsic E_p − L correlation scenario {E_p0 = 10^2.60 keV, σ_Ep = 0.30, α_A = 0.50}.

Fig. B.1. Corner plot from the MCMC exploration of fake observations generated from known inputs. The known true values are shown in green.

Appendix C: Fluence cut–sample size plane

The sample size and fluence cut of SHOALS do not allow us to discriminate between non-evolving (k_evol = 0) and strongly evolving (k_evol = 2) scenarios for the luminosity function of LGRBs. In this section we detail our methodology for determining the probability of discriminating between these two scenarios given a fluence cut and a sample size.

First, we cut both scenarios at a given fluence cut ℱ_cut. We then create a sub-sample from one scenario of size N_sub and compute the K–S test between the redshift distributions of this sub-sample and the sample of the other scenario. This K–S test is repeated for a number of bootstrap samples N_bs; the result is a distribution of p-values and D-statistics shown in the bottom panels of Fig. C.1. We define the probability of discriminating between the two scenarios at the 95% confidence level as the fraction of realisations with a p-value below 0.05; in the example of Fig. C.1, this probability would be 12%. We calculate this fraction for different values of N_sub and ℱ_cut and draw the contours to obtain the fluence cut–sample size plane shown in Fig. 8.

Fig. C.1. Test to distinguish between the two extreme scenarios in terms of cosmic evolution of the luminosity function. Top panel: redshift cumulative distribution function of the intrinsic Ep − L correlation scenarios cut at ℱcut = 10−6 ph cm−2 for kevol = 0 in dark blue and kevol = 2 in light blue. The kevol = 2 scenario is sub-sampled at Nsub = 117; the 95% confidence bound is calculated from bootstrapping. Bottom panel: distribution of p-values (right) and D-statistics (left) from the K–S tests performed on the bootstrapped samples. The back dashed vertical line shows a p-value of 0.05, below which it is possible to distinguish between the two distributions at the 95% confidence level.

Appendix D: Comparison with other models

For each model below, we generated the synthetic population with our own Monte Carlo code using N_GRB = 10⁶ LGRBs, assuming the same cosmology and the same functional forms for the intrinsic distributions (rate, luminosity function, etc.) as reported by the authors, and adopting the best fit parameters we found in the papers. We then confront the properties of the resulting population to the observational constraints used in the present study.

D.1. Daigne et al. (2006)

We used the best fit values presented in Table 1 of Daigne et al. (2006) for the six scenarios they explored which cover three different hypotheses concerning the LGRB efficiency (constant, mildly increasing and strongly increasing, labelled SFR1, SFR2, and SFR3, respectively) and two scenarios for the peak energy distribution (intrinsically correlated with the peak luminosity, labelled A for ‘Amati-like’, and independent log-normal, labelled LN). We used a Band function for the LGRB spectra with their empirical distributions for α and β. We used the same cosmology: Ω_m = 0.3, Ω_Λ = 0.7, and H₀ = 65 km s⁻¹ Mpc⁻¹. The observational constraints presented in Sect. 3 and the predictions from the model parameters are shown in Fig. D.1. It should be noted that our intensity constraint (leftmost panels in Fig. D.1) is the same as the one they used to adjust their population parameters; however, the other constraints are new. This comparison illustrates how using a redshift constraint, which was not possible in 2006, provides a clear advantage to better constrain the underlying intrinsic LGRB population.

Fig. D.1. Population generated using the model parameters of Daigne et al. (2006) compared to the observational constraints presented in Sect. 3. In the rightmost panels, the dashed lines correspond to the redshift distribution of the entire intrinsic population. Top panels: intrinsic Ep − L correlation scenarios and bottom panels: independent log-normal Ep scenarios.

D.2. Wanderman & Piran (2010)

We used the best fit values reported in Table 1 of Wanderman & Piran (2010) for the 1/2 bins case. We used a broken power law for the luminosity function, using their value L_min = 10⁵⁰ erg s⁻¹ and a fixed Band spectrum with E_p = 511 keV, α = −1 and β = −2.25. We used the same cosmology: Ω_m = 0.27, Ω_Λ = 0.73, and H₀ = 70 km s⁻¹ Mpc⁻¹. The observational constraints presented in Sect. 3 and the predictions from the model parameters are shown in Fig. D.2. It should be noted their population is adjusted on a sample of ∼120 Swift bursts; however, it reproduces very nicely the intensity constraint of observed CGRO/BATSE LGRBs if the normalisation is adjusted. The spectrum constraint is poorly reproduced because they used a fixed spectrum for all LGRBs, on the other hand the redshift constraint is surprisingly well reproduced, which suggests that their approximate criteria for the probability of redshift measurement as a function of gamma-ray peak flux is reasonable.

Fig. D.2. Population generated using the model parameters of Wanderman & Piran (2010) compared to the observational constraints presented in Sect. 3; the intensity constraint was renormalised.

D.3. Salvaterra et al. (2012)

We used the best fit values presented in Table 2 of Salvaterra et al. (2012) for the luminosity evolution scenario, the density evolution scenario and the metallicity threshold scenario. In all three scenarios we used the values for a broken power law luminosity function with L_min = 10⁴⁹ erg s⁻¹ and a Band spectrum with α = −1, β = −2.25 and with E_p following the intrinsic correlation of Eq. (9) with α_A = 0.49, E_p0 = 337 keV, σ_Ep = 028. We used the same cosmology: Ω_m = 0.3, Ω_Λ = 0.7, and H₀ = 70 km s⁻¹ Mpc⁻¹. It should be noted that the luminosity evolution scenario of Salvaterra et al. (2012) is slightly different in its definition than what we used in our model. In their implementation, it is L_cut, the break in the power law, that evolves with redshift and L_min stays fixed, whereas in our implementation defined in Eq. (8), the whole luminosity function evolves with redshift, including L_min and L_max. We therefore modified our code accordingly in order to properly reproduce their luminosity function and used the value of k_evol = 2.1 for the luminosity evolution scenario; for the density evolution scenario and for the metallicity threshold scenario (see the description and the definition of parameters in Salvaterra et al. 2012) we used δ_n = 1.7 and Z_th = 0.1, respectively. The observational constraints presented in Sect. 3 and the predictions from the model parameters are shown in Fig. D.3. In general, this population adjusted on the original Swift/BAT6 sample of 58 LGRBs and the log N − log P of Stern et al. (2001) reproduces well our observational constraints, despite some discrepancies at high and low peak flux. The simulated peak energy distribution is also lower on average than the one observed by Fermi/GBM. This illustrates the strength of simultaneously using constraints from CGRO/BATSE, Fermi/GBM and Swift/BAT that cover all major aspects of a population. We extend this approach in the present paper, adding a spectral constraint based on Fermi/GBM and using Swift data accumulated over a much much longer duration.

Fig. D.3. Model parameters of Salvaterra et al. (2012) applied to the observational constraints presented in Sect. 3.

D.4. Pescalli et al. (2016)

Pescalli et al. (2016) derived step functions for the LGRB formation rate and luminosity function, using a different approach than the previous studies cited above, based on the C⁻ method of Lynden-Bell (1971). In this work, contrary to Salvaterra et al. (2012), the authors constraints rely only on the luminosity and redshift distributions (and their join distribution: the luminosity-redshift plane) of the extended BAT6 sample of bright Swift LGRBs with $N_{15-150\mathrm{keV}}^{\mathrm{pk}}\ge 2.6\mathrm{ph}{\mathrm{s}}^{-1}{\mathrm{cm}}^{-2}$. We use the values they quote in Sect. 7 for a broken power law fit to their luminosity step function with L_min = 10⁴⁹ erg s⁻¹ and k_evol = 2.5. They do not quote any analytical form for the redshift distribution but they conclude it follows the star formation rate; we therefore used the form of Madau & Dickinson (2014). We used a Band spectrum with α = −1, β = −2.25 and with E_p following the intrinsic correlation of Eq. (9) with the values they derive in Table 1: α_A = 0.54, E_p0 = 309 keV, σ_Ep = 028. We used the same cosmology: Ω_m = 0.3, Ω_Λ = 0.7, and H₀ = 70 km s⁻¹ Mpc⁻¹. The observational constraints presented in Sect. 3 and the predictions from the model parameters are shown in Fig. D.4. Interestingly, we see in our intensity constraint that this population reproduces well the high peak flux regime (after renormalisation) but departs from the log N − log P towards lower peak fluxes. This comes as a surprise since the luminosity functions of Salvaterra et al. (2012) and Pescalli et al. (2016) appear consistent in the right panel of Fig. 1 of Pescalli et al. (2016). This is probably due to the difference in the low-luminosity slope between the two cases, which affects the number of faint bursts and thus the faint end of the log N − log P. The discrepancies in the log N − log P underline the advantage of constraining a population model on multiple complementary observational constraints that reach down to faint peak fluxes.

Fig. D.4. Model parameters of Pescalli et al. (2016) applied to the observational constraints presented in Sect. 3; the intensity constraint was renormalised.

All Tables

All Figures

Fig. 1. Distribution of Cvar and T90 obs for the GBM catalogue of Bhat et al. (2016) for burst with $N_{50-300\mathrm{keV}}^{\mathrm{pk}}\ge 0.9\mathrm{ph}{\mathrm{s}}^{-1}{\mathrm{cm}}^{-2}$. The data are binned in a 2D histogram whose shading represents the density of points; contours are also shown to guide the eye. The black dashed line corresponds to the linear fit of Eq. (11). The black filled lines in the side histograms are the 1D Gaussian kernel density estimations of the data.

In the text

Fig. 2. Observational constraints used in our population model. Panel a: intensity constraint; log N − log P diagram built from the offline search of CGRO/BATSE LGRBs of Stern et al. (2001), corrected for efficiency of detection at low fluxes. Panel b: spectral constraint; Ep number density distribution from the GBM spectral catalogue (Gruber et al. 2014) for long GRBs with $N_{50-300\mathrm{keV}}^{\mathrm{pk}}>0.9\mathrm{ph}{\mathrm{s}}^{-1}{\mathrm{cm}}^{-2}$. The actual histogram used in the fitting procedure is shown as black crosses while a finer-grained histogram is shown in light grey to guide the eye. The black line is a Gaussian kernel density estimation of the data. Panel c: redshift constraint; redshift number density distribution of the Swift/eBAT6 sample. The black line is a Gaussian kernel density estimation of the data.

In the text

Fig. 3. Efficiency-corrected log N − log P diagram for CGRO/BATSE from Stern et al. (2001) in black. The log N − log P of GBM is shown in green, adjusted to the one of CGRO/BATSE by multiplying by a constant whose value is obtained by minimising the χ2 between the two histograms over the grey shaded region. The black dashed vertical line represents the Npk cut to avoid biases due to faint flux incompleteness for our spectral constraint.

In the text

Fig. 4. Ep − L plane for the Swift/eBAT6 sample. The individual data points are colour-coded by redshift; the filled red line represents the Ep − Liso relation found by Pescalli et al. (2016) for this sample and the dashed red lines represents the scatter. A 2D Gaussian kernel density estimation of the data is shown as grey shaded contours. The black dashed, dot-dashed, and dotted lines represent the detection threshold for $N_{15-150\mathrm{keV}}^{\mathrm{pk}}=2.6\mathrm{ph}{\mathrm{s}}^{-1}{\mathrm{cm}}^{-2}$ and a fixed Band spectrum (α = 0.6, β = 2.5) at redshifts 0.3, 2, and 5, respectively. The side histograms represent the binned data and the colour of each bin represents the median redshift in that bin following the colour-code of the central panel; the black curve is the 1D Gaussian kernel density estimation.

In the text

Fig. 5. Best fits to the observational constraints from models with intrinsic Ep − L correlations and a luminosity function that does not evolve with redshift (kevol = 0). The red curve corresponds to the case where the LGRB rate follows the shape of CSFRD and the blue curve corresponds to the case where the LGRB efficiency is free to evolve with redshift. The black curves and grey histograms are the observational constraints described in Fig. 2. This figure is discussed in Sect. 4.

In the text

Fig. 6. Swift/eBAT6 Ep − L plane for models with no redshift evolution of the luminosity function (a, c) and with strong redshift evolution of the luminosity function (b, d). Top panels (a, b): independent log-normal Ep scenarios; bottom panels (c, d): intrinsic Ep − L correlation scenarios. The predictions for the entire intrinsic LGRB population is shown as grey contours and the predictions for the mock eBAT6 sample are shown in coloured filled contours. The observed eBAT6 sample is shown as coloured scatter points, the 2D Gaussian kernel density estimation shown in Fig. 4 is omitted for clarity. A full description of the observed data points, the side histograms and the dashed, dot-dashed and dotted curves of the figure is presented in Fig. 4.

In the text

	Fig. 7. Cumulative redshift distribution of the observed SHOALS sample in black; the 95% confidence bound calculated from bootstrapping is shown in grey and the upper limits are shown as arrows at the bottom of the plot. The predictions for the intrinsic Ep − L correlation scenarios with kevol = 0 and kevol = 2 are shown in blue and dashed light blue, respectively.
In the text

	Fig. 8. Fluence cut–sample size plane for the intrinsic Ep − L correlation scenario. The colour shading corresponds to the probability of discriminating between the kevol = 0 and the kevol = 2 scenarios at the 2σ (95%) confidence level. The star corresponds to the current SHOALS sample of Perley et al. (2016b).
In the text

Fig. 9. Luminosity cumulative distribution function in four different redshift bins for the intrinsic Ep − L correlation scenario. The full lines correspond to kevol = 2 and the dashed lines to kevol = 0. Leftmost panel: entire intrinsic LGRB population, centre left and centre right panels: bursts with $N_{15-150\mathrm{keV}}^{\mathrm{pk}}\ge 0.1\mathrm{ph}{\mathrm{s}}^{-1}{\mathrm{cm}}^{-2}$ and 2.6 ph s−1 cm−2 respectively. Rightmost panel: same as the centre right panel but adding the observed eBAT6 sample and where each curve was shifted by a constant C for clarity. The shaded area represents the 95% confidence region for the CDF of the observed eBAT6 sample, calculated using bootstrapping (without taking into account limits).

In the text

Fig. 10. Cosmic evolution of the LGRB rate and efficiency. Top panel: comoving LGRB rate density as a function of redshift for the best fitting models of the intrinsic Ep − L correlation scenarios. Bottom panel: LGRB efficiency η(z) as as function of redshift obtained with Eq. (6). The estimation of the fraction of star formation below Zth from Eq. (5) of Langer & Norman (2006) is shown, rescaled for comparison, in dashed. The shaded area represents the 95% confidence interval.

In the text

Fig. 11. All-sky rate of LGRBs at z ≥ 6, as a function of limiting peak flux for different energy bands in the intrinsic Ep − L correlation scenario. Each panel corresponds to a luminosity evolution scenario (from no evolution at kevol = 0 to strong evolution at kevol = 2). To obtain the actual observed rate for a given mission, one should take into account the field of view of the instrument and its duty cycle. The shaded area represents the 95% confidence interval.

In the text

Fig. A.1. Quality of Band fits to GRB spectra from GBM catalogue (Bhat et al. 2016). Left: spectral slopes α (top) and β (bottom) as a function of peak flux for the GBM spectral catalogue (Bhat et al. 2016). The entire catalogue of long GRBs above 0.9 ph s−1 cm−2 is in light-grey, the bursts complying with the ‘good’ criteria (i.e., small errors on the parameters) are shown in blue, the bursts for which Band is the best-fit spectral model are shown in orange. Right: reduced χ2 distribution for Band models from the GBM spectral catalogue (Bhat et al. 2016).

In the text

	Fig. B.1. Corner plot from the MCMC exploration of fake observations generated from known inputs. The known true values are shown in green.
In the text

Fig. C.1. Test to distinguish between the two extreme scenarios in terms of cosmic evolution of the luminosity function. Top panel: redshift cumulative distribution function of the intrinsic Ep − L correlation scenarios cut at ℱcut = 10−6 ph cm−2 for kevol = 0 in dark blue and kevol = 2 in light blue. The kevol = 2 scenario is sub-sampled at Nsub = 117; the 95% confidence bound is calculated from bootstrapping. Bottom panel: distribution of p-values (right) and D-statistics (left) from the K–S tests performed on the bootstrapped samples. The back dashed vertical line shows a p-value of 0.05, below which it is possible to distinguish between the two distributions at the 95% confidence level.

In the text

	Fig. D.1. Population generated using the model parameters of Daigne et al. (2006) compared to the observational constraints presented in Sect. 3. In the rightmost panels, the dashed lines correspond to the redshift distribution of the entire intrinsic population. Top panels: intrinsic Ep − L correlation scenarios and bottom panels: independent log-normal Ep scenarios.
In the text

	Fig. D.2. Population generated using the model parameters of Wanderman & Piran (2010) compared to the observational constraints presented in Sect. 3; the intensity constraint was renormalised.
In the text

	Fig. D.3. Model parameters of Salvaterra et al. (2012) applied to the observational constraints presented in Sect. 3.
In the text

	Fig. D.4. Model parameters of Pescalli et al. (2016) applied to the observational constraints presented in Sect. 3; the intensity constraint was renormalised.
In the text