Issue 
A&A
Volume 649, May 2021



Article Number  A166  
Number of page(s)  22  
Section  Extragalactic astronomy  
DOI  https://doi.org/10.1051/00046361/202039929  
Published online  02 June 2021 
Constraining the intrinsic population of long gammaray bursts: Implications for spectral correlations, cosmic evolution, and their use as tracers of star formation
^{1}
Sorbonne Université, CNRS, UMR7095, Institut d’Astrophysique de Paris, 75014 Paris, France
email: palmerio@iap.fr
^{2}
GEPI, Observatoire de Paris, PSL University, CNRS, 5 Place Jules Janssen, 92190 Meudon, France
Received:
17
November
2020
Accepted:
9
March
2021
Aims. Long gammaray bursts (LGRBs) have been shown to be powerful probes of the Universe, in particular for studying the star formation rate up to very high redshift (z ∼ 9). Since LGRBs are produced by only a small fraction of massive stars, it is paramount to have a good understanding of their underlying intrinsic population in order to use them as cosmological probes without introducing any unwanted bias. The goal of this work is to constrain and characterise this intrinsic population.
Methods. We developed a Monte Carlo model where each burst is described by its redshift and its properties at the peak of the light curve. We derived the best fit parameters by comparing our synthetic populations to carefully selected observational constraints based on the CGRO/BATSE, Fermi/GBM and Swift/BAT samples with appropriate flux thresholds. We explored different scenarios in terms of the cosmic evolution of the luminosity function and/or of the redshift distribution as well as including or not the presence of intrinsic spectralenergetics (E_{p} − L) correlations.
Results. We find that the existence of an intrinsic E_{p} − L correlation is preferred but with a shallower slope than observed (α_{A} ∼ 0.3) and a larger scatter (∼0.4 dex). We find a strong degeneracy between the cosmic evolution of the luminosity and of the LGRB rate, and show that a sample both larger and deeper than SHOALS by a factor of three is needed to lift this degeneracy.
Conclusions. The observed E_{p} − L correlation cannot be explained only by selection effects although these do play a role in shaping the observed relation. The degeneracy between the cosmic evolution of the luminosity function and of the redshift distribution of LGRBs should be included in the uncertainties of star formation rate estimates; these amount to a factor of 10 at z = 6 and up to a factor of 50 at z = 9.
Key words: methods: statistical / gammaray burst: general / galaxies: star formation
© J. T. Palmerio & F. Daigne 2021
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1. Introduction
Gammaray 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). Gammaray bursts exhibit two main radiative phases: (i) the prompt emission, which is commonly detected in the hard Xrays 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 Xrays, optical and radio (see e.g., Gehrels et al. 2009; Gehrels & Razzaque 2013). Gammaray 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, lowmass starforming 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 lowredshift 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 shortlived 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 redshiftcomplete observational samples of LGRBs as they require a rapid response from ground followup 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^{−6} 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 sourceframe peak energy distribution: lognormal 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 followup 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 luminosityredshift 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 wellmeasured spectral parameters, and a redshift constraint based on the extended BAT6 sample. Finally, we discuss the implications of our results in terms of intrinsic spectralenergetics 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^{1}: Ω_{m} = 0.27, Ω_{Λ} = 0.73, and H_{0} = 71 km s^{−1} Mpc^{−1}.
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 ‘sourceframe’ will be reserved to refer to properties as measured without the effect of cosmological redshift dilation (e.g., sourceframe E_{p}, versus observed E_{p}, noted E_{p obs}; when not specified the default is sourceframe).
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 sourceframe 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^{6} 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 timeintegrated 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 isotropicequivalent bolometric peak luminosity in units of erg s^{−1}, defined as , where L_{E} is the sourceframe power density in units of [erg s^{−1} keV^{−1}], (ii) E_{p}, the energy at which the E L_{E} spectrum peaks in the sourceframe, in units of keV, (iii) α, the lowenergy slope of the photon spectrum L_{E}/E, and (iv) β, the highenergy 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 timeintegrated properties that depend on the light curve of the LGRB are: (i) T_{90}, 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 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^{−1} cm^{−2}]) with:
where N_{Eobs} is the observed photon flux density (in units of [ph s^{−1} cm^{−2} keV^{−1}]) 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^{−1} cm^{−2}]) with:
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.
Summary of the properties and calculated quantities used to describe the LGRBs in our population model.
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; AmaralRogers 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:
with f(z) a broken exponential:
where z_{m} is the redshift of the break, a and b are the low and highredshift slopes, respectively, and 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:
where ṅ_{cc} is the comoving rate density of core collapses (in units of [yr^{−1} Mpc^{−3}]). We then assume that ṅ_{cc} is related to the CSFRD by:
where is the CSFRD in units of [M_{⊙} yr^{−1} Mpc^{−3}] and is the mean mass deduced from the stellar IMF:
Here, p_{cc}(z) is the fraction of formed stars ending with a core collapse, given by:
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:
In the following, we assume p_{cc} and 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:
with , 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
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):
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:
where A is a normalisation factor given by . 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 highluminosity 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^{2}. 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^{49} erg s^{−1}, which corresponds to the lowestluminosity burst in the Swift/eBAT6 sample (see Sect. 5.2 for a discussion on the lowluminosity population in our model); this value is similar to the ones assumed by other studies (typically taken between 10^{48} and 10^{50} erg s^{−1}). 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)^{kevol} as in Salvaterra et al. (2012). The redshiftevolving luminosity function is given by:
where the 1/(1 + z)^{kevol} prefactor 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 deevolved 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 signaltonoise ratio in the GBM catalogue (see left panels of Fig. A.1) and whose shape has only three parameters: the lowenergy slope α, the highenergy 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 isotropicequivalent 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 lognormal 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:
where E_{p0} is drawn from a lognormal distribution with scatter σ_{Ep}, and L_{0} is a constant fixed at 1.6 × 10^{52} erg s^{−1}. 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 lognormal E_{p} scenario assumes E_{p} is independent of L and follows a lognormal 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 lognormal 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^{−1} cm^{−2} to avoid faint flux incompleteness, see Sect. 3.2). This allows us to compute the probability density function of the intrinsic duration T_{90} by first fitting the observed distribution of T_{90 obs} from the GBM catalogue, cutting out LGRBs below 0.9 ph s^{−1} cm^{−2}, with a lognormal distribution, which yields a mean μ = 1.45 and standard deviation σ = 0.47. Then we assume the intrinsic distribution of T_{90} is lognormal with the same spread and correct for cosmological time dilation as , where 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 to the peak luminosity L. It can be estimated from the GBM catalogue with:
where 𝒩 is the photon fluence (in units of [ph cm^{−2}]) and N^{pk} is the peak photon flux (in units of [ph s^{−1} cm^{−2}]) defined in Eq. (1). The C_{var} − T_{90 obs} plane of the GBM catalogue with 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:
Fig. 1. Distribution of C_{var} and T_{90 obs} for the GBM catalogue of Bhat et al. (2016) for burst with . 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 lognormal distribution yielding μ = 0.04 and σ = 0.22. Thus to obtain C_{var}, we: (i) draw z and T_{90} as mentioned above and calculate T_{90 obs}, (ii) draw log from a lognormal 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_{90} 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_{90} and C_{var}; it is used as a crosscheck (see Sect. 3.5.2) rather than a constraint for this reason.
Summary of the various mock samples created for our population model.
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^{3}, 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; E_{p} number density distribution from the GBM spectral catalogue (Gruber et al. 2014) for long GRBs with . The actual histogram used in the fitting procedure is shown as black crosses while a finergrained 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 isotropicequivalent 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 offline search of all CGRO/BATSE records, they conclude that there is no turnover down to a peak flux of about ∼0.1 ph s^{−1} cm^{−2}. Using the catalogue^{4} 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^{5} (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 χ^{2} 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^{−1} cm^{−2}, 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. Efficiencycorrected 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 N^{pk} 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^{6} 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 socalled ‘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 Xshooter (Vernet et al. 2011) has mostly remedied this. It is thus crucial to use a wellcontrolled 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 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 , 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 (Goldstein et al. 2013), we can use this corrected log N − log P to obtain the total observed allsky LGRB rate:
for BATSE LGRBs (i.e., brighter than 0.067 ph s^{−1} cm^{−2} in the 50–300 keV range). We can use to normalise our population with:
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 . Defining the intrinsic rate of all LGRBs above L_{min} as:
where N_{GRB} is the number of LGRBs in our simulated population, we can obtain with:
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:
The observed rate of LGRBs for a given sample is then simply the total rate times the average exposure factor, f_{exp}:
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 isotropicequivalent luminosity of the bursts and their peak energy. This is relevant, in particular when trying to distinguish between scenarios with intrinsic correlated or independent lognormal 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. E_{p} − L plane for the Swift/eBAT6 sample. The individual data points are colourcoded by redshift; the filled red line represents the E_{p} − L_{iso} 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, dotdashed, and dotted lines represent the detection threshold for 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 colourcode 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_{90} 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 lognormal 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 nonevolving 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 overpredicts the number of bright and lowredshift 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 lognormal 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 E_{p} − L correlations and a luminosity function that does not evolve with redshift (k_{evol} = 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 lognormal 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.
Best fit parameter values in the case of a lognormal E_{p} scenario.
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 lognormal 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 lognormal 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 lognormal 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 highredshift slope b varies strongly with k_{evol} and is similar in both E_{p} scenarios and finally the lowredshift 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 increases with k_{evol} in the independent lognormal 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 lognormal 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 luminosityevolving and efficiencyevolving 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 lognormal 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 lognormal E_{p} scenarios with mild or no evolution, we find a dearth of highL 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 lognormal 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 E_{p} − 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 lognormal E_{p} scenarios; bottom panels (c, d): intrinsic E_{p} − 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, dotdashed 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, dotdashed, 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
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 lowluminosity bursts, which are also mostly lowE_{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 subclasses of GRBs, such as lowluminosity GRBs (LLGRBs, Liang et al. 2007; Virgili et al. 2009; Stanway et al. 2014), Xray Rich GRBs or Xray Flashes (XRRGRBs 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 LLGRBs and XRRGRBs/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 XRRGRBs 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 XRRGRBs 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 lowluminosity tail of the classical LGRB population but rather a fullfledged 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 lowerenergy 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), HiZGUNDAM (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 gammaray 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 (LloydRonning 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^{−2}] in the same energy interval [E_{min, obs};E_{max, obs}]) with:
This assumes that the timeintegrated prompt spectrum is the same as the peakbrightness prompt spectrum which is a strong additional assumption. However, using bursts from the observed Fermi/GBM^{bright} sample (with ) we find that the peak energy of these two spectra are tightly correlated: . 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^{−2}] between E_{min, obs} and E_{max, obs}) which can be consequently calculated from:
where F^{pk} is the peak energy flux (in units of [erg s^{−1} cm^{−2}] in the same energy interval) given by Eq. (3). The observed SHOALS sample is defined as having ℱ_{15−150 keV} > 10^{−6} erg cm^{−2} 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 (pvalues 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 E_{p} − L correlation scenarios with k_{evol} = 0 and k_{evol} = 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^{−6} erg cm^{−2}), 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 KStests; 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 E_{p} − L correlation scenario. The colour shading corresponds to the probability of discriminating between the k_{evol} = 0 and the k_{evol} = 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 lowerluminosity 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 and 2.6 ph s^{−1} cm^{−2}, 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 E_{p} − L correlation scenario. The full lines correspond to k_{evol} = 2 and the dashed lines to k_{evol} = 0. Leftmost panel: entire intrinsic LGRB population, centre left and centre right panels: bursts with 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 highredshift 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 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 starformation 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 E_{p} − 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 Z_{th} 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 for the intrinsic E_{p} − L correlation scenarios and for the independent lognormal 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^{49} erg s^{−1}; the following equation can be used to extrapolate to lower values of L_{min} for a Schechter or powerlaw function assuming that the additional LGRBs are not detected in the observed samples:
where p is the slope of the luminosity function and erg s^{−1} 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 lowluminosities 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 η_{0} = (3−6) × 10^{−6}, 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; LloydRonning 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 restframe 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 nonevolving 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 highredshift LGRBs, a primary goal for using the full potential of LGRBs as a cosmological tool. In Fig. 11 we show the expected allsky 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 highredshift LGRBs is higher for softer energy bands, and depends in particular on the lowerenergy 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 codedmask telescope ECLAIRs is designed to trigger at these lower energies with a lowenergy threshold of 4 keV (Godet et al. 2014). This holds even more for other future missions, such as THESEUS^{7}, 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 lowerenergy threshold (potentially as low as 0.5 keV), or the HiZGUNDAM (Yonetoku et al. 2014), with its 1–10 keV wide field Xray 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^{−2} ph s^{−1} cm^{−2}; this implies that the whole of the population is seen and going to lower peak fluxes will not increase the proportion of highredshift LGRBs. This could be an interesting test to distinguish between different luminosity evolution scenarios. However, the possibility of getting unbiased, redshiftcomplete samples of LGRBs down to such low peak fluxes remains to this day unlikely.
Fig. 11. Allsky rate of LGRBs at z ≥ 6, as a function of limiting peak flux for different energy bands in the intrinsic E_{p} − L correlation scenario. Each panel corresponds to a luminosity evolution scenario (from no evolution at k_{evol} = 0 to strong evolution at k_{evol} = 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 lognormal distribution independently of other properties (independent lognormal 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 spectralenergetics 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., XRRGRBs, XRFs), which are expected to be weak following the spectralenergetics correlation; this may imply a different slope of the luminosity function at low luminosities or a separate population alltogether, 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 nonevolving 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 nonevolving 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 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 lowenergy 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 lowenergy threshold with a good sensitivity could help to partially lift the degeneracy between luminosity and rate evolution.
These are the same cosmological parameters as in Daigne et al. (2006), Wanderman & Piran (2010), Salvaterra et al. (2012), and Pescalli et al. (2016), allowing a direct comparison of our work with these previous studies. We checked that adopting updated parameters, Ω_{m} = 0.31, Ω_{Λ} = 0.69, and H_{0} = 68 km s^{−1} Mpc^{−1} (Planck Collaboration XIII 2016), has a negligible impact on our results.
To be rigorous, what we present in panel a of Fig. 2 is a log R − log N diagram, where N is the peak photon flux N^{pk} and R is the rate of GRBs per Δlog N^{pk} bin; for simplicity we refer to this as a log N − log P.
Acknowledgments
This work has been partially supported by the French Space Agency (CNES). JTP would like to thank Sébastien Carassou, Tom Charnock and JensKristian Krogager for invaluable discussion regarding the statistical framework used in this paper. The corner plots used in this paper for representing multidimensional datasets were performed with the PYTHON corner.py code (ForemanMackey 2016). Parts of the results in this work make use of the colormaps in the CMasher package (van der Velden 2020).
References
 AmaralRogers, A., Willingale, R., & O’Brien, P. T. 2016, MNRAS, 464, 2000 [Google Scholar]
 Amati, L. 2006, MNRAS, 372, 233 [Google Scholar]
 Amati, L., Frontera, F., Tavani, M., et al. 2002, A&A, 390, 81 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Amati, L., O’Brien, P., Götz, D., et al. 2018, AdSpR, 62, 191 [Google Scholar]
 Atteia, J. L., Heussaff, V., Dezalay, J. P., et al. 2017, ApJ, 837, 119 [Google Scholar]
 Band, D. L., & Preece, R. D. 2005, ApJ, 627, 319 [Google Scholar]
 Band, D., Matteson, J., Ford, L., et al. 1993, ApJ, 413, 281 [NASA ADS] [CrossRef] [Google Scholar]
 Barlow, R., & Beeston, C. 1993, Comput. Phys. Commun., 77, 219 [Google Scholar]
 Barraud, C., Daigne, F., Mochkovitch, R., & Atteia, J. L. 2005, A&A, 440, 809 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Behroozi, P. S., Wechsler, R. H., & Conroy, C. 2013, ApJ, 770, 57 [Google Scholar]
 Beloborodov, A. M., & Mészáros, P. 2017, Space Sci. Rev., 207, 87 [Google Scholar]
 Bhat, P. N., Meegan, C. A., von Kienlin, A., et al. 2016, ApJS, 223, 28 [Google Scholar]
 Blanchard, P. K., Berger, E., & Fong, W.F. 2016, ApJ, 817, 144 [Google Scholar]
 Bloom, J. S., Kulkarni, S. R., Djorgovski, S. G., et al. 1999, Nature, 401, 453 [Google Scholar]
 Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2015, ApJ, 803, 34 [Google Scholar]
 Bouwens, R. J., Oesch, P. A., Labbé, I., et al. 2016, ApJ, 830, 67 [Google Scholar]
 Bromm, V., & Loeb, A. 2006, ApJ, 642, 382 [Google Scholar]
 Butler, N. R., & Kocevski, D. 2007, ApJ, 668, 400 [Google Scholar]
 Butler, N. R., Kocevski, D., Bloom, J. S., & Curtis, J. L. 2007, ApJ, 671, 656 [Google Scholar]
 Carassou, S., de Lapparent, V., Bertin, E., & Le Borgne, D. 2017, A&A, 605, A9 [EDP Sciences] [Google Scholar]
 Chrimes, A. A., Stanway, E. R., & Eldridge, J. J. 2020, MNRAS, 491, 3479 [NASA ADS] [CrossRef] [Google Scholar]
 Cucchiara, A., Levan, A. J., Fox, D. B., et al. 2011, ApJ, 736, 7 [Google Scholar]
 Daigne, F., & Mochkovitch, R. 1998, MNRAS, 296, 275 [NASA ADS] [CrossRef] [Google Scholar]
 Daigne, F., Rossi, E. M., & Mochkovitch, R. 2006, MNRAS, 372, 1034 [Google Scholar]
 de Mink, S. E., Cantiello, M., Langer, N., et al. 2009, A&A, 497, 243 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Drovandi, C. C., Pettitt, A. N., & Lee, A. 2015, Stat. Sci., 30, 72 [Google Scholar]
 ForemanMackey, D. 2016, J. Open Source Softw., 1, 24 [Google Scholar]
 Frontera, F., Amati, L., Guidorzi, C., Landi, R., & in’t Zand, J. 2012, ApJ, 754, 138 [Google Scholar]
 Fruchter, A. S., Levan, A. J., Strolger, L., et al. 2006, Nature, 441, 463 [Google Scholar]
 Fryer, C. L., & Heger, A. 2005, ApJ, 623, 302 [Google Scholar]
 Fynbo, J. P. U., Starling, R. L. C., Ledoux, C., et al. 2006, A&A, 451, L47 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gehrels, N., & Razzaque, S. 2013, Front. Phys., 8, 661 [Google Scholar]
 Gehrels, N., RamirezRuiz, E., & Fox, D. 2009, ARA&A, 47, 567 [Google Scholar]
 Ghirlanda, G., Nava, L., Ghisellini, G., Firmani, C., & Cabrera, J. I. 2008, MNRAS, 387, 319 [Google Scholar]
 Ghirlanda, G., Ghisellini, G., Nava, L., et al. 2012, MNRAS, 422, 2553 [Google Scholar]
 Ghirlanda, G., Ghisellini, G., Salvaterra, R., et al. 2013, MNRAS, 428, 1410 [Google Scholar]
 Ghirlanda, G., Salvaterra, R., Ghisellini, G., et al. 2015, MNRAS, 448, 2514 [Google Scholar]
 Giannios, D., & Spruit, H. C. 2005, A&A, 430, 1 [EDP Sciences] [Google Scholar]
 Godet, O., Nasser, G., Atteia, J., et al. 2014, in Space Telescopes and Instrumentation 2014: Ultraviolet to Gamma Ray, SPIE Conf. Ser., 9144, 914424 [Google Scholar]
 Goldstein, A., Preece, R. D., Mallozzi, R. S., et al. 2013, ApJS, 208, 21 [Google Scholar]
 Graham, J. F., & Fruchter, A. S. 2017, ApJ, 834, 170 [Google Scholar]
 Greiner, J., Krühler, T., Klose, S., et al. 2011, A&A, 526, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gruber, D., Goldstein, A., von Ahlefeld, V. W., et al. 2014, ApJS, 211, 12 [Google Scholar]
 Hao, J.M., Cao, L., Lu, Y.J., et al. 2020, ApJS, 248, 21 [Google Scholar]
 Harrison, F. A., Bloom, J. S., Frail, D. A., et al. 1999, ApJ, 523, L121 [Google Scholar]
 Hartoog, O. E., Malesani, D., Fynbo, J. P. U., et al. 2015, A&A, 580, A139 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Heussaff, V. 2015, PhD Theses, Université Paul Sabatier  Toulouse III, France [Google Scholar]
 Heussaff, V., Atteia, J.L., & Zolnierowski, Y. 2013, A&A, 557, A100 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hjorth, J., Sollerman, J., Møller, P., et al. 2003, Nature, 423, 847 [Google Scholar]
 Hjorth, J., Malesani, D., Jakobsson, P., et al. 2012, ApJ, 756, 187 [Google Scholar]
 Izzo, L., Thöne, C. C., Schulze, S., et al. 2017, MNRAS, 472, 4480 [Google Scholar]
 Izzo, L., Auchettl, K., Hjorth, J., et al. 2020, A&A, 639, L11 [EDP Sciences] [Google Scholar]
 Jakobsson, P., Levan, A., Fynbo, J. P. U., et al. 2006, A&A, 447, 897 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Japelj, J., Vergani, S. D., Salvaterra, R., et al. 2016, A&A, 590, A129 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Jenke, P. A., Linares, M., Connaughton, V., et al. 2016, ApJ, 826, 228 [Google Scholar]
 Kaneko, Y., Magdalena González, M., Preece, R. D., Dingus, B. L., & Briggs, M. S. 2008, ApJ, 677, 1168 [Google Scholar]
 Katsukura, D., Sakamoto, T., Tashiro, M. S., & Terada, Y. 2020, ApJ, 889, 110 [Google Scholar]
 Kirkpatrick, S., Gelatt, C. D., & Vecchi, M. P. 1983, Science, 220, 671 [Google Scholar]
 Kistler, M. D., Yüksel, H., Beacom, J. F., & Stanek, K. Z. 2008, ApJ, 673, L119 [Google Scholar]
 Kobayashi, S., Piran, T., & Sari, R. 1997, ApJ, 490, 92 [Google Scholar]
 Kommers, J. M., Lewin, W. H. G., Kouveliotou, C., et al. 2000, ApJ, 533, 696 [Google Scholar]
 Kouveliotou, C., Meegan, C. A., Fishman, G. J., et al. 1993, ApJ, 413, L101 [Google Scholar]
 Krühler, T., Malesani, D., Fynbo, J. P. U., et al. 2015, A&A, 581, A125 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Krühler, T., Kuncarayakti, H., Schady, P., et al. 2017, A&A, 602, A85 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kumar, P., & Zhang, B. 2015, Phys. Rep., 561, 1 [Google Scholar]
 La Barbera, F., Ferreras, I., Vazdekis, A., et al. 2013, MNRAS, 433, 3017 [NASA ADS] [CrossRef] [Google Scholar]
 Lamb, D. Q., & Reichart, D. E. 2000, ApJ, 536, 1 [Google Scholar]
 Lan, G.X., Zeng, H.D., Wei, J.J., & Wu, X.F. 2019, MNRAS, 488, 4607 [Google Scholar]
 Langer, N., & Norman, C. A. 2006, ApJ, 638, L63 [NASA ADS] [CrossRef] [Google Scholar]
 Le Floc’h, E., Duc, P.A., Mirabel, I. F., et al. 2003, A&A, 400, 499 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Le Floc’h, E., Charmandaris, V., Forrest, W. J., et al. 2006, ApJ, 642, 636 [NASA ADS] [CrossRef] [Google Scholar]
 Liang, E.W., Zhang, B.B., & Zhang, B. 2007, ApJ, 670, 565 [Google Scholar]
 LloydRonning, N., Hurtado, V. U., Aykutalp, A., Johnson, J., & Ceccobello, C. 2020, MNRAS, 494, 4371 [Google Scholar]
 Lu, R.J., Wei, J.J., Liang, E.W., et al. 2012, ApJ, 756, 112 [Google Scholar]
 Lyman, J. D., Levan, A. J., Tanvir, N. R., et al. 2017, MNRAS, 467, 1795 [NASA ADS] [Google Scholar]
 LyndenBell, D. 1971, MNRAS, 155, 95 [Google Scholar]
 Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415 [NASA ADS] [CrossRef] [Google Scholar]
 Mazets, E. P., Golenetskii, S. V., Il’Inskii, V. N., et al. 1981, Astrophys. Space Sci., 80, 3 [Google Scholar]
 McKinney, J. C., & Uzdensky, D. A. 2012, MNRAS, 419, 573 [Google Scholar]
 Melandri, A., Sbarufatti, B., D’Avanzo, P., et al. 2012, MNRAS, 421, 1265 [Google Scholar]
 Melandri, A., Malesani, D. B., Izzo, L., et al. 2019, MNRAS, 490, 5366 [Google Scholar]
 Mészáros, P., & Rees, M. J. 2000, ApJ, 530, 292 [Google Scholar]
 Mochkovitch, R., & Nava, L. 2015, A&A, 577, A31 [EDP Sciences] [Google Scholar]
 Nakar, E., & Piran, T. 2005, MNRAS, 360, L73 [Google Scholar]
 Nava, L., Ghirlanda, G., Ghisellini, G., & Firmani, C. 2008, MNRAS, 391, 639 [Google Scholar]
 Nava, L., Salvaterra, R., Ghirlanda, G., et al. 2012, MNRAS, 421, 1256 [Google Scholar]
 Oesch, P. A., Bouwens, R. J., Illingworth, G. D., et al. 2014, ApJ, 786, 108 [Google Scholar]
 Oesch, P. A., Bouwens, R. J., Illingworth, G. D., et al. 2015, ApJ, 808, 104 [Google Scholar]
 Paczynski, B. 1998, ApJ, 494, L45 [Google Scholar]
 Palmerio, J. T., Vergani, S. D., Salvaterra, R., et al. 2019, A&A, 623, A26 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Perley, D. A., Levan, A. J., Tanvir, N. R., et al. 2013, ApJ, 778, 128 [Google Scholar]
 Perley, D. A., Krühler, T., Schulze, S., et al. 2016a, ApJ, 817, 7 [NASA ADS] [CrossRef] [Google Scholar]
 Perley, D. A., Tanvir, N. R., Hjorth, J., et al. 2016b, ApJ, 817, 8 [Google Scholar]
 Pescalli, A., Ghirlanda, G., Salvaterra, R., et al. 2016, A&A, 587, A40 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Piran, T. 2005, Rev. Mod. Phys., 76, 1143 [Google Scholar]
 Planck Collaboration XIII. 2016, A&A, 594, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Podsiadlowski, P., Ivanova, N., Justham, S., & Rappaport, S. 2010, MNRAS, 406, 840 [NASA ADS] [Google Scholar]
 Porciani, C., & Madau, P. 2001, ApJ, 548, 522 [Google Scholar]
 Preece, R. D., Briggs, M. S., Mallozzi, R. S., et al. 2000, ApJS, 126, 19 [Google Scholar]
 Rees, M. J., & Meszaros, P. 1994, ApJ, 430, L93 [Google Scholar]
 Robertson, B. E., & Ellis, R. S. 2012, ApJ, 744, 95 [Google Scholar]
 Sakamoto, T., Lamb, D. Q., Kawai, N., et al. 2005, ApJ, 629, 311 [Google Scholar]
 Sakamoto, T., Hullinger, D., Sato, G., et al. 2008, ApJ, 679, 570 [Google Scholar]
 Salafia, O. S., Barbieri, C., Ascenzi, S., & Toffano, M. 2020, A&A, 636, A105 [CrossRef] [EDP Sciences] [Google Scholar]
 Salpeter, E. E. 1955, ApJ, 121, 161 [Google Scholar]
 Salvaterra, R., Della Valle, M., Campana, S., et al. 2009, Nature, 461, 1258 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Salvaterra, R., Campana, S., Vergani, S. D., et al. 2012, ApJ, 749, 68 [Google Scholar]
 Savaglio, S., Glazebrook, K., & Le Borgne, D. 2009, ApJ, 691, 182 [Google Scholar]
 Schechter, P. 1976, ApJ, 203, 297 [Google Scholar]
 Selsing, J., Malesani, D., Goldoni, P., et al. 2019, A&A, 623, A92 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Shahmoradi, A., & Nemiroff, R. J. 2011, MNRAS, 411, 1843 [Google Scholar]
 Springel, V., & Hernquist, L. 2003, MNRAS, 339, 289 [Google Scholar]
 Spruit, H. C., Daigne, F., & Drenkhahn, G. 2001, A&A, 369, 694 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Stanway, E. R., Levan, A. J., Tanvir, N., et al. 2014, MNRAS, 446, 3911 [Google Scholar]
 Stern, B. E., Tikhomirova, Y., Kompaneets, D., Svensson, R., & Poutanen, J. 2001, ApJ, 563, 80 [Google Scholar]
 Svensson, K. M., Levan, A. J., Tanvir, N. R., Fruchter, A. S., & Strolger, L.G. 2010, MNRAS, 405, 57 [NASA ADS] [Google Scholar]
 Tanvir, N. R., Fox, D. B., Levan, A. J., et al. 2009, Nature, 461, 1254 [Google Scholar]
 Tanvir, N. R., Fynbo, J. P. U., de Ugarte Postigo, A., et al. 2019, MNRAS, 483, 5380 [Google Scholar]
 Usov, V. V. 1994, MNRAS, 267, 1035 [Google Scholar]
 van der Velden, E. 2020, J. Open Source Softw., 5, 2004 [Google Scholar]
 Vangioni, E., Olive, K. A., Prestegard, T., et al. 2015, MNRAS, 447, 2575 [Google Scholar]
 Vergani, S. D., Salvaterra, R., Japelj, J., et al. 2015, A&A, 581, A102 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Vergani, S. D., Palmerio, J., Salvaterra, R., et al. 2017, A&A, 599, A120 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Vernet, J., Dekker, H., D’Odorico, S., et al. 2011, A&A, 536, A105 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Vielfaure, J. B., Vergani, S. D., Japelj, J., et al. 2020, A&A, 641, A30 [EDP Sciences] [Google Scholar]
 Virgili, F. J., Liang, E.W., & Zhang, B. 2009, MNRAS, 392, 91 [Google Scholar]
 von Kienlin, A., Meegan, C. A., Paciesas, W. S., et al. 2014, ApJS, 211, 13 [Google Scholar]
 Wanderman, D., & Piran, T. 2010, MNRAS, 406, 1944 [NASA ADS] [Google Scholar]
 Wei, J., Cordier, B., Antier, S., et al. 2016, ArXiv eprints [arXiv:1610.06892] [Google Scholar]
 White, N. E. 2020, The Gamow Explorer: A GammaRay Burst Mission to Study the High Redshift Universe [Google Scholar]
 Woosley, S. E. 1993, ApJ, 405, 273 [Google Scholar]
 Yonetoku, D., Murakami, T., Nakamura, T., et al. 2004, ApJ, 609, 935 [Google Scholar]
 Yonetoku, D., Murakami, T., Tsutsui, R., et al. 2010, PASJ, 62, 1495 [NASA ADS] [Google Scholar]
 Yonetoku, D., Mihara, T., Sawano, T., et al. 2014, in Space Telescopes and Instrumentation 2014: Ultraviolet to Gamma Ray, eds. T. Takahashi, J. W. A. den Herder, & M. Bautz, SPIE Conf. Ser., 9144, 91442S [Google Scholar]
 Zhang, B., & Yan, H. 2011, ApJ, 726, 90 [Google Scholar]
Appendix A: Band function
The Band function (Band et al. 1993) is given by:
where α is the lowenergy slope, β the highenergy slope, E_{p} the peak of the E L_{E} spectrum, and A is a normalisation that satisfies . 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^{−1} cm^{−2}, illustrates that even when Band is not the best fitting model it still provides a good fit to the data (i.e., the reduced χ^{2} 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 lightgrey, 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 bestfit 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:
The likelihood for the entire histogram becomes:
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 loglikelihood is thus:
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^{−3}. 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 loglikelihood. 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 τ = τ_{0} that decreases at a userdefined 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^{53} erg s^{−1}, 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. 
Summary of the flat prior bounds used for the parameters of our population model.
Appendix C: Fluence cut–sample size plane
The sample size and fluence cut of SHOALS do not allow us to discriminate between nonevolving (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 subsample from one scenario of size N_{sub} and compute the K–S test between the redshift distributions of this subsample 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 pvalues and Dstatistics 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 pvalue 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 E_{p} − L correlation scenarios cut at ℱ_{cut} = 10^{−6} ph cm^{−2} for k_{evol} = 0 in dark blue and k_{evol} = 2 in light blue. The k_{evol} = 2 scenario is subsampled at N_{sub} = 117; the 95% confidence bound is calculated from bootstrapping. Bottom panel: distribution of pvalues (right) and Dstatistics (left) from the K–S tests performed on the bootstrapped samples. The back dashed vertical line shows a pvalue 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^{6} 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 ‘Amatilike’, and independent lognormal, 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_{0} = 65 km s^{−1} Mpc^{−1}. 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 E_{p} − L correlation scenarios and bottom panels: independent lognormal E_{p} 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^{50} erg s^{−1} 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_{0} = 70 km s^{−1} Mpc^{−1}. 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 gammaray 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^{49} erg s^{−1} 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_{0} = 70 km s^{−1} Mpc^{−1}. 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 LyndenBell (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 luminosityredshift plane) of the extended BAT6 sample of bright Swift LGRBs with . We use the values they quote in Sect. 7 for a broken power law fit to their luminosity step function with L_{min} = 10^{49} erg s^{−1} 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_{0} = 70 km s^{−1} Mpc^{−1}. 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 lowluminosity 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
Summary of the properties and calculated quantities used to describe the LGRBs in our population model.
Summary of the flat prior bounds used for the parameters of our population model.
All Figures
Fig. 1. Distribution of C_{var} and T_{90 obs} for the GBM catalogue of Bhat et al. (2016) for burst with . 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; E_{p} number density distribution from the GBM spectral catalogue (Gruber et al. 2014) for long GRBs with . The actual histogram used in the fitting procedure is shown as black crosses while a finergrained 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. Efficiencycorrected 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 N^{pk} cut to avoid biases due to faint flux incompleteness for our spectral constraint. 

In the text 
Fig. 4. E_{p} − L plane for the Swift/eBAT6 sample. The individual data points are colourcoded by redshift; the filled red line represents the E_{p} − L_{iso} 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, dotdashed, and dotted lines represent the detection threshold for 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 colourcode 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 E_{p} − L correlations and a luminosity function that does not evolve with redshift (k_{evol} = 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 E_{p} − 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 lognormal E_{p} scenarios; bottom panels (c, d): intrinsic E_{p} − 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, dotdashed 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 E_{p} − L correlation scenarios with k_{evol} = 0 and k_{evol} = 2 are shown in blue and dashed light blue, respectively. 

In the text 
Fig. 8. Fluence cut–sample size plane for the intrinsic E_{p} − L correlation scenario. The colour shading corresponds to the probability of discriminating between the k_{evol} = 0 and the k_{evol} = 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 E_{p} − L correlation scenario. The full lines correspond to k_{evol} = 2 and the dashed lines to k_{evol} = 0. Leftmost panel: entire intrinsic LGRB population, centre left and centre right panels: bursts with 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 E_{p} − 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 Z_{th} 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. Allsky rate of LGRBs at z ≥ 6, as a function of limiting peak flux for different energy bands in the intrinsic E_{p} − L correlation scenario. Each panel corresponds to a luminosity evolution scenario (from no evolution at k_{evol} = 0 to strong evolution at k_{evol} = 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 lightgrey, 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 bestfit 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 E_{p} − L correlation scenarios cut at ℱ_{cut} = 10^{−6} ph cm^{−2} for k_{evol} = 0 in dark blue and k_{evol} = 2 in light blue. The k_{evol} = 2 scenario is subsampled at N_{sub} = 117; the 95% confidence bound is calculated from bootstrapping. Bottom panel: distribution of pvalues (right) and Dstatistics (left) from the K–S tests performed on the bootstrapped samples. The back dashed vertical line shows a pvalue 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 E_{p} − L correlation scenarios and bottom panels: independent lognormal E_{p} 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 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.