Constraining the intrinsic population of Long Gamma-Ray Bursts: implications for spectral correlations, cosmic evolution and their use as tracers of star formation

Long Gamma-Ray Bursts (LGRBs) have been shown to be powerful probes of the Universe, in particular to study the star formation rate up to very high redshift ($z \sim 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. We developed a Monte Carlo model where each burst is described by its redshift and its properties at the peak of the lightcurve. 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 cosmic evolution of the luminosity function and/or of the redshift distribution as well as including or not the presence of intrinsic spectral-energetics (Ep-L) correlations. We find that the existence of an intrinsic Ep-L correlation is preferred but with a shallower slope than observed($\alpha_A \sim 0.3$) and a larger scatter ($\sim 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. The observed We conclude that Ep-L correlation cannot be explained only by selection effects although these do play a role in shaping the observed relation. The degeneracy between 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$.


Introduction
Gamma-Ray Bursts (GRBs) are the most powerful electromagnetic phenomena, associated to the relativistic ejection following the birth of a stellar mass compact object (black hole or magnetar, see e.g. Piran 2005;Kumar & Zhang 2015). GRBs exhibit two main radiative phases: (i) the prompt emission which is commonly detected in the hard X-rays and soft γ-rays (1 keV -10 MeV) and usually lasts a few hundreds of milliseconds to a few hundreds of seconds; (ii) the afterglow emission, which is initially bright, decays rapidly, and peaks successively in X-rays, optical and radio (see e.g. Gehrels et al. 2009;Gehrels & Razzaque 2013). GRBs 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, Salvaterra et al. 2009;Tanvir et al. 2009 and using photometry only at z ∼ 9, 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 * E-mail: palmerio@iap.fr to cosmological expansion (Lamb & Reichart 2000); 1 day after the prompt emission on Earth corresponds to 6 hours in the source frame at z = 3 and 2 hours at z = 10, essentially catching the afterglow earlier in its light curve and thus brighter as redshift increases.
Among the various classes of GRBs, the case of long GRBs (LGRBs, which have a duration of the prompt emission longer than ∼ 2 s, Mazets et al. 1981;Kouveliotou et al. 1993) is the most promising for the study of the distant Universe. These are the most frequent GRBs and both theoretical progenitor models ("collapsar" model, Woosley 1993;Paczynski 1998) and observations have firmly associated them with the collapse of certain massive stars. Indeed, they generally occur in faint, blue, low-mass star-forming galaxies (Le Floc'h et al. 2003;Savaglio et al. 2009;Palmerio et al. 2019), and in the UV-brightest regions of their hosts (Fruchter et al. 2006;Svensson et al. 2010;Lyman et al. 2017). In addition, the majority of low-redshift LGRBs are found in association with core-collapse supernovae Hjorth et al. 2003;Krühler et al. 2017;Izzo et al. 2017;Melandri et al. 2019;Izzo et al. 2020). Due to the short-lived nature of their massive star progenitors, LGRBs are expected to occur up to very high redshift, possibly in association with the first generation of stars (Bromm & Loeb 2006), and thus can be used as lighthouses to study galaxies (e.g. Le Floc'h et al. 2006;Perley et al. 2013;Vergani et al. 2017) and the intergalactic medium at high redshift (e.g. Fynbo et al. 2006;Hartoog et al. 2015;Selsing et al. 2019;Tanvir et al. 2019;Vielfaure et al. 2020).
From their association with certain massive stars, LGRBs as a population are considered potential tracers of star formation (e.g. Lamb & Reichart 2000;Porciani & Madau 2001;Kistler et al. 2008;Robertson & Ellis 2012;Vangioni et al. 2015). However, the precise link between star formation and LGRBs remains shrouded in uncertainty for two main reasons. The first reason is that this link depends on many factors which 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 follow-up and significant telescope time to get meaningful statistics. In practice we often have to chose between either biased or 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 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 App. 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 source-frame peak energy distribution: log-Normal or correlated to the peak luminosity. Using the observed redshift distribution from early Swift results of Jakobsson et al. (2006) as a cross-check, they found evidence for an increas-ing LGRB production efficiency with redshift. With the advent of Swift and thanks to dedicated follow-up campaigns, samples of GRBs with a significant fraction of redshift measurements were obtained, allowing for two new approaches to emerge. The first approach is to model the redshift recovery efficiency as in Wanderman & Piran (2010). These authors directly inverted the observed luminosity-redshift distribution of a sample of 120 Swift LGRBs, assuming a fixed Band spectrum as in Porciani & Madau (2001), to obtain the intrinsic luminosity function and redshift distribution. They found some evidence that the LGRB rate does not follow the CSFRD, in agreement with Daigne et al. (2006). Using a different approach, Salvaterra et al. (2012) designed a complete and unbiased sample of 58 bright Swift LGRBs named BAT6 (later extended to 99 LGRBs by Pescalli et al. 2016), essentially paying the cost of sample size in order to achieve high redshift completeness and avoid the uncertainties caused by modelling the redshift recovery efficiency. They found evidence for strong redshift evolution of the luminosity function although they noted the degeneracy with the redshift evolution of the LGRB efficiency and suggested both might be at play.
Building on these works, we present in this paper our own population model which combines the latest, largest and most diverse datasets of LGRBs available to date with currently relevant scenarios in order to investigate some of the questions raised above. We fit our population model parameters using an intensity constraint based on a large sample of CGRO/BATSE LGRBs including faint events, a spectral constraint based on a sample of bright Fermi/GBM LGRBs with well measured spectral parameters, and a redshift constraint based on the extended BAT6 sample. Finally, we discuss the implications of our results in terms of intrinsic spectral-energetics correlations, of soft bursts, of cosmic evolution of the LGRB luminosity and/or rate, of the use of LGRBs as tracers of star formation at high redshift and of 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 while 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 are provided in the Appendix. All errors are reported at the 1σ confidence level unless stated otherwise. We use a standard cosmology: Ω m = 0.27, Ω Λ = 0.73, and H 0 = 71 km s −1 Mpc −1 .

Modelling the intrinsic LGRB population
We should specify to lift any ambiguity that throughout the rest of this paper the term intrinsic is used to qualify the entire, true LGRB population, without any selection criteria (e.g. intrinsic redshift distribution versus redshift distribution of the Swift sample) while the term source-frame will be reserved to refer to properties as measured without the effect of cosmological redshift dilation (e.g. source-frame Table 1: Summary of the properties and calculated quantities used to describe the LGRBs in our population model. In each table, the Pk type refers to quantities defined at peak brightness of the LGRB while the Ti signifies the quantity is time-integrated (i.e. it depends on the lightcurve of the LGRB). The first table lists the properties describing a single LGRB, always defined in the source frame. The top three properties' distributions are adjusted by comparing the population to the observational constraints described in Sect. 3. They are completed by the next two properties, with fixed distribution (see Sect. 2.4) to entirely describe the emission at the peak brightness. The last two properties are only used for additional cross-checks described in Sect. 3.5. Their distributions are computed (see Sect. 2.5 The peak energy of the E 2 N Eobs spectrum in the observer frame, E pobs = E p /(1 + z). N pk ph s −1 cm −2 Pk The peak photon flux, calculated from z, L, E p , α, β for any E min,obs ; E max,obs (see Eq. 1) F pk erg s −1 cm −2 Pk The peak energy flux, calculated from z, L, E p , α, β for any E min,obs ; E max,obs (see Eq. 3) T 90obs s Ti The duration over which 5 to 95 % of photons are received in the observer frame, The photon fluence, calculated from z, L, E p , α, β, T 90 and C var for any E min,obs ; E max,obs (see Eq. 15) F erg cm −2 Ti The energy fluence, calculated from z, L, E p , α, β, T 90 and C var for any E min,obs ; E max,obs (see Eq. 16) E iso erg Ti The isotropic-equivalent bolometric energy, E iso = C var T 90 L E p , versus observed E p , noted E pobs ; when not specified the default is source-frame).
Each burst in our model is described by the following properties: a redshift, a peak luminosity and a spectrum. For the majority of the work presented in this paper, these properties alone are enough to compute all the observed quantities we are interested in, and in particular the peak flux. However, for some specific studies described in Sect. 3.5, other observed quantities such as duration or fluence are necessary. For this purpose, we add two properties to fully describe an LGRB: a source-frame duration and a variability coefficient which 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.

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 which 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: -L, the isotropic-equivalent bolometric peak luminosity in units of erg s −1 , defined as L = ∞ 0 L E dE, where L E is the source-frame power density in units of [erg s −1 keV −1 ]. -E p , the energy at which the E L E spectrum peaks in the source-frame, in units of keV. -α, the low-energy slope of the photon spectrum L E /E. -β, the high-energy slope of the photon spectrum L E /E. -In addition to these four quantities, the intrinsic shape of the photon spectrum L E /E is assumed to be given by a Band function (Band et al. 1993).
The peak brightness of an LGRB depends on the timescale with which its lightcurve 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 impact significantly the calculated values of the peak flux; a rough estimate puts this correction between 1, 7, 13, 18% at redshifts 0.1, 1, 3, 6 respectively (Heussaff 2015).
The two time-integrated properties which depend on the lightcurve of the LGRB are: -T 90 , the duration over which 5 to 95% of photons are emitted in the source frame. -C var , the variability coefficient defined as the mean lu-minosityL 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 E obs 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 Tab. 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.

Intrinsic redshift distribution
The intrinsic redshift distribution is perhaps one of the most important ingredients of an LGRB population model because any constraint on this distribution has critical implications for the identification of LGRB progenitors and the use of LGRBs as tracers of star formation. However, it remains poorly constrained due to the small fraction of GRBs with measured redshifts (∼ 30 % for Swift GRBs) leading to significant uncertainties regarding its shape; for these reasons many authors have used various functional forms to represent it (see e.g. Porciani & Madau 2001;Daigne et al. 2006;Wanderman & Piran 2010;Salvaterra et al. 2012;Amaral-Rogers et al. 2016).
In this work we chose to use a simple functional form, with only 3 free parameters, which could adequately represent the cosmic star formation rate density (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, and a and b are the low-and high-redshift slopes respectively andṅ 0 GRB is a normalisation given by our model (see Sect. 3.4).
We can relateṅ GRB to the cosmic star formation rate density by introducing the LGRB efficiency η(z) defined aṡ n GRB (z) = η(z)ṅ cc (z) 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 cosmic star formation rate density by: whereρ * (z) is the cosmic star formation rate density in units of [M yr −1 Mpc −3 ] andm is the mean mass deduced from the stellar initial mass function (IMF): p cc (z) is the fraction of formed stars ending with a corecollapse, 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 andm to be constant with redshift, meaning any discrepancy between the redshift distribution of LGRBs and the cosmic star formation rate density 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 cosmic star formation rate density is then taken to have a similar shape as the broken exponential defined above: withρ 0 * = 2.8 × 10 −2 M yr −1 Mpc −3 , a * = 1.1, b * = −0.57 and z m, * = 1.9; these values are obtained by fitting the broken exponential form of Eq. 4 to the slightly more complex functional form of Springel & Hernquist (2003) with the parameter values of model 3 from Vangioni et al. (2015) which are deduced of the observed cosmic formation rate extracted by (Behroozi et al. 2013) from galaxy luminosity functions measurements including high redshift measurements at z 8 − 10 (Oesch et al. 2014;Oesch et al. 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 cosmic star formation rate density values a * , b * and z m, * , we can obtain the evolving LGRB efficiency η(z) by rearranging Eq. 5:

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 normalization factor given by 1 = ∞ 0 φ(L) dL. In practice this means we are actually working with a probability density but we will continue refering to it as a luminosity function out of convenience. This functional form has the advantage of having only 3 parameters: the minimum luminosity L min , the break luminosity L * and the slope p. It is worth noting that this form is very similar to a simple power law with the difference that the break at the high-luminosity end is smooth rather than sudden which is slightly more realistic (see the discussion in Atteia et al. 2017). We also tested luminosity models with a broken power law but found little success in deriving meaningful constraints on all parameters simultaneously 1 . 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 lowest luminosity burst in the Swift/eBAT6 sample (see Sect. 5.2 for a discussion on the low-luminosity population in our model); this value is similar to the ones assumed by other studies (typically taken between 10 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) k evol as in Salvaterra et al. (2012). The redshift-evolving luminosity function is given by: where the 1/(1+z) k evol pre-factor comes from the condition that the probability density remain normalised to unity. In this case, the values of the parameters quoted are always given for the de-evolved luminosity function, at z = 0.

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 S/N in the GBM catalogue (see left pannels of Fig. A.1) and whose shape has only 3 parameters: the low-energy slope α, the high-energy slope β, and the peak energy E p . The Band function is presented in more detail in App. A.

Intrinsic peak energy distribution
In the past, various authors have found relations between the spectral and energetic properties of GRBs and in particular between the peak energy E p of the GRB spectrum and the isotropic-equivalent energy E iso (Amati et al. 2002;Amati 2006;Lu et al. 2012) or luminosity L iso (Yonetoku et al. 2004;Yonetoku et al. 2010;Frontera et al. 2012). Some authors have suggested that these correlations are caused by strong selection effects (Nakar & Piran 2005;Band & Preece 2005;Shahmoradi & Nemiroff 2011;Heussaff et al. 2013), while others have argued that selection effects do not suffice to explain the observed correlation Nava et al. 2008;Ghirlanda et al. 2012). For our models we tested two different scenarios regarding the E p distribution of LGRBs: a scenario with an intrinsic correlation between E p and L iso and a scenario where the intrinsic E p distribution follows a log-Normal distribution independently of L iso .
Several population models (e.g. Salvaterra et al. 2012;Pescalli et al. 2016;Lan et al. 2019) assume that the observed correlation between the peak luminosity of a burst and its peak energy is intrinsic. The correlation has 3 free parameters and can be written as: where E p0 is drawn from a log-Normal distribution with scatter σ E p , 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, σ E p = 0.28 (Pescalli et al. 2016). In the intrinsic E p − L correlation scenario studied in the present paper, we decided to relax this condition and leave the parameters of this relation free to vary to allow for an intrinsic E p − L correlation to be affected by additional selection effects to shape the observed correlation.
The independent log-Normal E p scenario assumes E p is independent of L and follows a log-Normal distribution with 2 free parameters: the mean E p0 [keV] and the spread σ E p [dex]. Note that the correlated scenario becomes the independent log-Normal scenario when α A goes to 0. The motivation behind this scenario is to see whether we can obtain an observed E p − L correlation without imposing any intrinsic correlation.

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).

Intrinsic duration and variability coefficient distributions
The intrinsic duration and variability coefficients are used solely for cross-checks 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 in-  Bhat et al. (2016) for burst with N pk 50−300 keV ≥ 0.9 ph s −1 cm −2 . The data are binned in a 2D histogram whose shading represents the density of points; contours are also shown to guide the eye. The black dashed line corresponds to the linear fit of Eq. 11. The black filled lines in the side histograms are the 1D Gaussian kernel density estimations of the data.
completeness, 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 90obs from the GBM catalogue, cutting out LGRBs below 0.9 ph s −1 cm −2 , with a log-Normal distribution which yields a mean µ = 1.45 and standard deviation σ = 0.47. Then we assume the intrinsic distribution of T 90 is log-Normal with the same spread and correct for cosmological time dilation as µ − log(1 +z GBM ), wherez GBM is the median redshift of our mock Fermi/GBM bright sample.
The variability coefficient C var is defined as the ratio of the mean luminosityL to the peak luminosity L. It can be estimated from the GBM catalogue with: where N 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 90obs plane of the GBM catalogue with N pk 50−300 keV ≥ 0.9 ph s −1 cm −2 is shown in Fig. 1 in the case of a Band spectral fit using the flux and fluence values in the 50-300 keV band. We fit the correlation between C var and T 90obs with a linear regression and find: We crosschecked these values on different catalogues (BATSE 5B, Goldstein et al. 2013) and different energy bands (GBM 10-1000 keV, Bhat et al. 2016) and find that the slope is within 10% around −0.4. We fit the decorrelated C var with a log-Normal distribution yielding µ = 0.04 and σ = 0.22. Thus to obtain C var we: (i) draw z and T 90 as

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 Tab. 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 Tab. 2. The first three mock samples of Tab. 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 cross-check (see Sect. 3.5.2) rather than a constraint for this reason.

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 which is validated given the good quality fits described in Sect. 4.

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 2 , represented in panel (a) of Fig. 2.
This diagram is a good way to estimate the peak isotropic-equivalent luminosity function of GRBs, however there is a difficulty residing in the fact that while peak fluxes are proportional to the luminosity, they also depend on redshift. This means a burst with a high peak flux could be low luminosity at low redshift, or high luminosity at high redshift. Fruitful studies (Kommers et al. 2000;Stern et al. 2001) have focused on the turnover at low peak flux, trying to determine if it is real (i.e. due to a minimum luminosity of GRBs) or if it is caused by the lower efficiency of detectors at these fluxes. In the case of Stern et al. (2001), who went down to the lowest peak flux limit by conducting an off-line search of all CGRO/BATSE records, they conclude that there is no turnover down to a peak flux of about ∼ 0.1 ph s −1 cm −2 . Using the catalogue 3 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 10 objects in each bin for proper Gaussian statistics to be applicable.

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 pobs 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 3rd Fermi/GBM catalogue 4 (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 normalized 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 min- imization on the bright bins (shown within the grey shaded area in Fig. 3) to find the best value. This normalization 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 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 pobs distribution that is unbiased from faint flux incompleteness, at the price of sample size. The resulting E pobs distribution is shown in panel (b) of Fig. 2.

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 5 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 effect and biases which are difficult to model (see however Wanderman & Piran 2010). For instance, the ability to measure a redshift for GRBs relies fundamentally on the capacity to locate it, which biases this distribution against so-called dark bursts (Greiner et al. 2011;Melandri et al. 2012) which exhibit highly extinguished optical afterglows. Another selection effect, called the redshift desert, is due to the fact that most emission and common absorption lines are shifted outside the window of optical spectrographs around z ∼ 2, although the advent of newer spectrographs such as X-Shooter (Vernet et al. 2011) have 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 there-5 http://www.mpe.mpg.de/~jcg/grbgen.html  Stern et al. (2001) in black. The log Nlog P of GBM is shown in green, adjusted to the one of CGRO/BATSE by multiplying by a constant whose value is obtained by minimizing 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.
fore used the redshift distribution from the extended BAT6 sample.
The Swift/BAT6 sample ) is a complete sample of Swift LGRBs with a selection based on the peak flux N pk 15−150 keV > 2.6 ph s −1 cm −2 and favorable 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 which do not depend on the redshift of the LGRB: rticle number, page 8 of 24 J. T. Palmerio & F. Daigne: Constraining the intrinsic population of LGRBs -The burst must be well localized by Swift/XRT and the information was distributed quickly -There is low galactic foreground extinction (A V < 0.5) -The burst declination is between -70 • and +70 • -The burst's angular distance to the sun is greater that 55 • -There are no nearby bright stars This results in 58 LGRBs for the original BAT6, later extended to 99 LGRBs by Pescalli et al. (2016). This extended BAT6 sample (Swift/eBAT6) has 82 bursts with redshifts, yielding an 83% redshift completeness, and its redshift distribution is shown in panel (c) of Fig. 2. It should be noted that while the eBAT6 bursts are statistically representative of the population of Swift bursts with N pk 15−150 keV > 2.6 ph s −1 cm −2 , this population constitutes only the brightest 25% of Swift bursts.

Event rates for our mock samples
One of the advantages of using the CGRO/BATSE log Nlog P of Stern et al. (2001) is that they derived a detection efficiency for all bursts within the field of view. Applying this efficiency correction while also taking into account the mean observed solid angle of the sky Ω BATSE = 8.17 sr and the live time of the search T BATSE live = 6.54 yr (Goldstein et al. 2013), we can use this corrected log N-log P to obtain the total observed all-sky LGRB rate: LGRB yr −1 in 4 π for BATSE LGRBs (i.e. brighter than 0.067 ph s −1 cm −2 in the 50-300 keV range). We can use R tot BATSE 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 R tot BATSE . 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ṅ 0 GRB 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.

Additional cross-checks
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 but see Sect. 5.1 and 5.3 for more details. Redshift (z) Fig. 4: E p − L plane for the Swift/eBAT6 sample. The individual data points are colour-coded 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, dot-dashed and dotted lines represent the detection threshold for N pk 15−150 keV = 2.6 ph s −1 cm −2 and a fixed Band spectrum (α = 0.6, β = 2.5) at redshifts 0.3, 2, and 5 respectively. The side histograms represent the binned data and the colour of each bin represents the median redshift in that bin following the colour-code of the central panel; the black curve is the 1D Gaussian kernel density estimation.
The Swift/eBAT6 sample is used as a redshift constraint, but there is more information to be extracted from this complete sample. More specifically, we use the E p − L plane since it contains information about the observed correlation between the isotropic-equivalent luminosity of the bursts and their peak energy. This is relevant, in particular when trying to distinguish between scenarios with intrinsic correlated or independent log-Normal E p distribu-tions (see Sect. 3.2). Figure 4 shows the E p − L plane for the Swift/eBAT6 sample (data is 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, σ E p = 0.28, their fit is represented by the blue 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).

The SHOALS redshift distribution
The Swift Gamma-Ray Burst 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 which 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 cross-check a posteriori to make sure the best fit populations also adequately represent SHOALS.

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 App. B. Due to computational considerations we tried to reduce the number of free parameters as much as possible in a justifiable way. For instance, in the independent log-Normal E p scenario, we noticed the value of the standard deviation σ E p 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 4 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 Tab. 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 cosmic star formation rate density), 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 low-redshift LGRBs. This is an important result connected to the physics of LGRB progenitors and was already found by Daigne et al. (2006) after two years of Swift observations. It was then confirmed by following studies (Wanderman & Piran 2010;Salvaterra et al. 2012) and we confirm it again on larger, more precise samples and find that it is true both in the independent log-Normal E p scenario and the scenario with intrinsic E p − L correlations. Our other models can then be used to better understand the underlying evolution (efficiency and/or luminosity) responsible for this discrepancy.
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 Tab. 3 for the independent log-Normal E p scenario and in Tab. 4 for the intrinsic E p − L correlation scenario; the quality of the fits in each case is comparable.
-The scenarios where the redshift distribution is fixed proportional to the cosmic star formation rate density find strong evolution of the luminosity function (k evol = 1.6 ± 0.1 for the independent log-Normal E p and k evol = 2.1 ± 0.1 for the intrinsic E p − L scenario). The value of k evol = 2.1 ± 0.1 is in agreement with previous studies which find k evol = 2.1 ± 0.6, ) and k evol = 2.5, (Pescalli et al. 2016) using similar E p scenarios. -Results obtained with an LGRB rate strictly proportional to the cosmic star formation rate are consistent with results obtained by fixing k evol and leaving the LGRB rate free to vary. In the intrinsic E p − L case, this corresponds to a strong cosmic evolution of the luminosity function (k evol = 2) and in the independent log-Normal E p case, this corresponds to an intermediary case between k evol = 1 and k evol = 2. For this reason in the following, the discussion will focus on the scenarios where the LGRB rate is left free to vary. -The intrinsic E p − L correlation scenarios find an intrinsic slope of the correlation significantly shallower than observed (α A 0.3 versus ∼ 0.5 observed) and a larger dispersion (σ E p 0.45 versus ∼ 0.28), close to the best value derived in the independent log-Normal E p scenario. The value of these best parameters are independent of the value of the scenario for the cosmic evolution of the luminosity and/or the rate.
-The behaviour of the parameters of the luminosity function and the redshift distribution is similar in both scenarios; as k evol increases, L * and the slopes of the redshift distribution, a and b, decrease. Comparing the two scenarios for the peak energy E p , we find very similar results regarding the luminosity function for the same value of k evol : p differs by less than 0.1, L * is typically bigger by a factor ∼ 2 in the intrinsic E p − L correlation scenario. We also find similar results regarding the redshift distribution for the same value of k evol : the break redshift z m 2.1 − 2.2 independently of k evol and of the E p scenario, the high redshift slope b varies strongly with k evol and is similar in both E p scenarios and finally the low redshift slope a also varies strongly with k evol but is typically steeper by +0.2 in the intrinsic E p − L correlation scenario. -Finally, the results are comparable for both scenarios regarding the peak energy distribution, which is expected −0.08 † The parameter b in these two cases has a multi-peaked marginalized posterior distribution; the median and 1 σ errors reported here are not necessarily representative of the best fitting value but are quoted for simplicity.  : 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 cosmic star formation rate density 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.
given that σ E p is the same and that the correlation slope α A is shallow. -Interestingly, the local LGRB rateṅ 0 GRB increases with k evol in the independent log-Normal E p scenarios while it stays roughly constant in the intrinsic E p − L correlation scenarios.
We conclude that we cannot distinguish between an independent log-Normal E p scenario and an intrinsic E p − L correlation scenario solely using the constraints presented in Sect. 3; we discuss this in more detail in Sect. 5.1 using additional cross-checks. 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 efficiency-evolving scenarios. This degeneracy has been an issue in LGRB population models based on smaller, complete samples of bright LGRBs from Swift (see e.g. Salvaterra et al. 2012;Pescalli et al. 2016). Despite using fainter and larger samples, we find that lifting this degeneracy remains difficult; we discuss this in more detail using additional cross-checks in Sect. 5.3.

Is there an intrinsic spectral-energetic correlation?
The observed E p − L plane of the extended BAT6 sample is a powerful tool to compare the predictions of the intrinsic E p − L correlation and the independent log-Normal E p scenarios. There has been significant debate concerning the causes of this observed correlation, whether due to observational selection effects or not (e.g. Nakar 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 rticle number, page 12 of 24 plane is shown for our mock eBAT6 sample in Fig. 6 for our best fit models with the two extreme values k evol = 0 and 2 (no/strong luminosity evolution) and the two scenarios regarding the peak energy E p . Looking at this plane for the independent log-Normal E p scenarios with mild or no evolution, we find a dearth of high-L bursts compared to the observations (panel (a) of Fig. 6); in the case of strong evolution, this discrepancy is less pronounced (panel (b) of Fig. 6). On the other hand, the luminosity distributions of the intrinsic E p − L correlation scenarios agree well with the observations for all values of k evol (panels (c) and (d) of Fig. 6). For this reason this scenario is preferred and the following discussions will be focused on this scenario, although we checked that the general trend was the same for the independent log-Normal E p scenario. This suggests that the observed E p − L correlation cannot be explained only from selection effects, in agreement with Ghirlanda et al. (2008); Nava et al. (2008); Ghirlanda et al. (2012).
Nonetheless, looking at the best fit values in Tab. 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 dashed, dot-dashed and dotted 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 E p (500 − 600 keV) L 1.6 10 52 erg/s 0.3 (14)

The population of soft bursts
It is remarkable that in all scenarios, the intrinsic distribution of LGRBs (shown as grey contours in Fig. 6) is dominated by largely undetected low-luminosity bursts, which are also mostly low-E p bursts in the preferred scenario with a moderate intrinsic E p − L correlation, in agreement with previous studies (Daigne et al. 2006). This could be linked to observed subclasses of GRBs such as low-luminosity GRBs (LL-GRBs, Liang et al. 2007;Virgili et al. 2009;Stanway et al. 2014), X-ray Rich GRBs or X-ray Flashes (XRR-GRBs and XRFs respectively, Barraud et al. 2005;Sakamoto et al. 2005Sakamoto et al. , 2008. With the spectral range and sensitivity of current instruments, their observational samples are still quite small (Virgili et al. 2009), making the precise description of their properties challenging. For XRFs, the γ-ray spectra are not very well characterised as the E p is often below the threshold of detectors and furthermore, the afterglows of XRR/XRFs are rarely observed (Sakamoto et al. 2008).
If the intrinsic spectrum-luminosity correlation is indeed valid, one would expect a link between LL-GRBs and XRR-GRBs/XRFs. However we note that the observed frequency of these events is larger than the predictions of our model. For instance, Sakamoto et al. (2005) proposed a classification based on the softness S defined as the ratio of energy fluences in the 2-30 keV and 30-400 keV bands. Among bursts detected both by the HETE2/WXM and FREGATE instruments, they find 22% of classical GRBs with S < 0.3, 42% of XRR-GRBs with 0.3 < S < 1 and 36% of XRFs with S < 1. For the corresponding mock sample (see Tab. 2) in our model, the fraction of XRFs is much lower (∼ 4%) and even XRR-GRBs are underrepresented (∼ 30%). This suggests that either the luminosity function of LGRBs has a different slope at low luminosity which is unconstrained by our model or that XRFs are not just the low-luminosity tail of the classical LGRB population but rather a full-fledged distinct population of their own. These two scenarios have different implications in terms of emission mechanisms and progenitor physics and distinguishing between them would require modelling the population of XRFs in details. Such a dedicated study remains difficult to this day due to the small size of samples of such soft events often without any afterglow or redshift measurements (but see Jenke et al. 2016;Katsukura et al. 2020). We can hope that the situation will improve in the future thanks to a new generation of satellites with a lower energy threshold and good localisation capabilities such as SVOM/ECLAIRs (4 -120 keV, Godet et al. 2014;Wei et al. 2016) and THESEUS (2 keV -20 MeV, Amati et al. 2018).

Cosmic evolution of
LGRBs: rate or luminosity?
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 cross-check 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 90obs following the prescription presented in Sect. 2.5 and, using the peak photon flux N pk following Eq. 1, we get the photon fluence N (in units of [ph cm −2 ] in the same energy interval  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). E min,obs ; E max,obs ) with: This assumes that the time-integrated prompt spectrum is the same as the peak-brightness prompt spectrum which is a strong additional assumption. However, using bursts from the observed Fermi/GBM bright sample (with N pk 50−300 keV ≥ 0.9 ph s −1 cm −2 ) we find that the peak energy of these two spectra are tightly correlated: log(E Pk p /E T i p ) = 0.04 ± 0.2. Additionally, we compared the peak flux -fluence planes of our mock samples with observations and found good agreement suggesting that our method is sound.
The SHOALS sample selection is based on the energy fluence F (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 F 15−150 keV > 10 −6 erg cm −2 and favorable observing conditions (see Perley et al. 2016a for more details). We construct our simulated SHOALS sample using this energy fluence threshold and look at the predictions regarding the redshift distribution. We restrict ourselves to scenarios with k evol = 0 (no evolution) and k evol = 2 (strong evolution) as they are the ones that will exhibit the most differences, the resulting distributions are shown in Fig. 7. All populations that fit the observational constraints described in Sect. 3 predict a mock SHOALS redshift distribution that cannot be distinguished from the observed SHOALS at the 95% confidence level as is attested by performing K-S tests between the cumulative distributions (p-values range from 0.7 to 0.95).
This means that even with the sample size and fluence cut of SHOALS (117 objects, F 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 KS-tests; the results are shown in Fig. 8 while the detailed methodology is shown in App. 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.
Another approach to lifting the degeneracy is to look directly at the luminosity distribution in different redshift bins. Fig. 9 shows the luminosity CDF in four redshift bins for the intrinsic E p − L correlation scenario at different peak flux cuts. The scenario with k evol = 0 is shown in dashed while k evol = 2 is shown in full. As expected, if one had access to the entire population (leftmost panel), it would be easy to distinguish between no luminosity evolution and strong luminosity evolution since the CDF would be the same in all redshift bins for the case k evol = 0. The difficulty resides in the fact that as we increase the peak flux cut (as shown in middle panel of Fig. 9), we remove mostly lower luminosity bursts and shift the distributions towards higher luminosities, essentially creating the same effect as luminosity evolution. This is illustrated in the center left and center right panels of Fig. 9 where the same plots are shown but for bursts with N pk 15−150 keV ≥ 0.1 ph s −1 cm −2 and 2.6 ph s −1 cm −2 respectively. The rightmost panel of Fig. 9 shows the same curves as in the center 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 favor 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.

The
LGRB rate and its connection to the cosmic star formation rate 5.4.1. The LGRB production efficiency by massive stars The top panel of Fig. 10 shows the LGRB comoving rate density as a function of redshift for the intrinsic LGRB population in the intrinsic E p − L correlation scenario for different values of k evol . The case of k evol = 2 follows closely the star formation rate density, while the high redshift slope gets shallower as k evol decreases. Models with lower values of k evol have higher global rates of LGRBs. Using Eq. 6, we can derive the LGRB efficiency as a function of redshift, shown in the bottom panel of Fig. 10. Assuming p cc andm 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;Japelj et al. 2016;Perley et al. 2016a;Vergani et al. 2017;Palmerio et al. 2019) which find that metallicity is the driving factor behind the LGRB production efficiency. The fraction of star-formation happening below a given metallicity (from Eq. 5 of Langer & Norman 2006) is shown in the bottom panel of Fig. 10 as dashed lines, arbitrarily scaled. We can see some similarity between our derived efficiency and these curves, although the behaviour at high redshift is different. This might be an indication of an additional break at z ≥ 6 in the redshift distribution of LGRBs, although with the amount of data available at these redshifts to date, it would be hard to constrain. Another possibility is that the effect of metallicity is more complicated than a simple threshold. Indeed, metallicity can also play an indirect role by influencing the stellar IMF (e.g. La Barbera et al. 2013) or the fraction of binary progenitors (Chrimes et al. 2020).

The Local LGRB rate
The local LGRB rate density for LGRBs pointing towards us we derived with our population model isṅ 0 GRB = 0.7 − 0.8 yr −1 Gpc −3 for the intrinsic E p − L correlation scenarios andṅ 0 GRB = 1.0 − 1.5 yr −1 Gpc −3 for the independent log-Normal E p scenarios (see Tab. 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 power-law function assuming that the additional LGRBs are not detected in the observed samples: where p is the slope of the luminosity function and L ref min = 5 × 10 49 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 App. D). We should note however that the slope at low-luminosities may be different (see discussion in Sect. 5.2) and therefore any extrapolation to lower L min should keep this in mind.
Using Eq. 5 we find η 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;Lloyd-Ronning et al. 2020), we get 0.8−1.6 LGRBs pointing in any direction per thousand core-collapse at z = 0 and up to a factor 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.

Measuring the CSFRD with LGRBs ?
Due to their association with massive stars, LGRBs have been used as probes of the cosmic star formation rate density (Kistler et al. 2008;Robertson & Ellis 2012;Hao et al. 2020), in particular up to high redshift where estimations from rest-frame UV measurements are plagued by dust uncertainties (Bouwens et al. , 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.

Strategy for high redshift detection
Using our intrinsic LGRB populations, we can investigate which type of strategies are most effective at observing high log L [erg/s] + C N pk 15−150 keV ≥ 2.6 ph/s/cm 2 C = 0 C = 1 C = 2 C = 3 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. The leftmost panel shows the entire intrinsic LGRB population, while the center left and center right panels show bursts with N pk 15−150 keV ≥ 0.1 ph s −1 cm −2 and 2.6 ph s −1 cm −2 respectively. The rightmost panel is the same as the center 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).
redshift LGRBs, a primary goal for using the full potential of LGRBs as a cosmological tool. In Fig. 11 we show the expected all-sky yearly rate of LGRBs at z ≥ 6 as a function of limiting peak flux for different energy bands. We see that, at a given peak flux limit, the rate of high redshift LGRBs is higher for softer energy bands, and depends in particular on the lower energy threshold of the band, in line with the results of Ghirlanda et al. (2015). This bodes well for the SVOM mission to be launched in 2022 (Wei et al. 2016), whose coded-mask telescope ECLAIRs is designed to trigger at these lower energies with a low-energy threshold of 4 keV (Godet et al. 2014), and even more for the THE-SEUS 6 mission in discussion for the next decade, whose XGIS instrument will detect GRBs in the 2 keV -20 MeV energy band with a larger effective area (Amati et al. 2018). 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 high redshift LGRBs. This could be an interesting test to distinguish between different luminosity evolution scenarios. However, the possibility of getting unbiased, redshift-complete samples of LGRBs down to such low peak fluxes remains to this day unlikely.

Conclusion
We created a model for the intrinsic population of LGRBs which we constrained by comparing to carefully selected observational constraints, specifically tailored to constrain all major aspects of the population. We explored two main scenarios concerning the distribution of the peak energy E p : one where the peak energy follows a log-Normal distribution independently of other properties (independent log-Normal E p ) and one where the peak energy is correlated to the luminosity (intrinsic E p − L correlation). In addition, in each case we explored four scenarios for the evolution of the luminosity function: no evolution (k evol = 0), mild 6 https://www.isdc.unige.ch/theseus/ 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 cross-checks 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. XRR-GRBs, XRFs) which are expected to be weak following the spectral-energetics correlation; this may imply a different slope of the luminosity function at low luminosities or a separate population all-together, either case is difficult to constrain with current samples. -We confirm the degeneracy between cosmic evolution of the LGRB rate and of the luminosity function. Under a strong cosmic evolution of the luminosity function (k evol = 2) we find that the LGRB rate follows the cosmic star formation rate; however under a nonevolving luminosity function we find a increasingly higher LGRB production efficiency from stars at high redshift. -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 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.
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ṅ 0 GRB ∼ 0.8 yr −1 Gpc −3 in the preferred intrinsic E p − L scenario which translates to ∼ 1 LGRB pointing in any direction per thousand core-collapse at z = 0, assuming a mean jet opening angle of θ jet 5 • . -We confirm that the best strategy for detecting high redshift bursts is to use detectors with a soft low-energy threshold of ∼ a few keV such as the ECLAIRs telescope on-board the SVOM satellite to be launched in 2022 or the XGIS instrument of the future THESEUS mission. We find that the combination of such a low energy threshold with a good sensitivity could help to partially lift the degeneracy between luminosity and rate evolution. LGRBs/yr in 4π at z LGRBs/yr in 4π at z 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.
where α is the low energy slope, β the high energy slope, E p the peak of the E L E spectrum, E c = α−β α+2 E p and A is a normalisation that satisfies 1 = ∞ 0 f (E) dE. 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).

Appendix B.1: Likelihood
We decided to follow a methodology already used for galaxy catalogs 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, o i and s i are the observed and simulated number count in bin i respectively. The log-likelihood is thus: where the factor ln(o i !) is neglected since it is a constant and our goal is to maximize ln L. 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. multi-plying 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 3.1) has the strongest weight of our constraints due to and β (bottom) as a function of peak flux for the GBM spectral catalogue (Bhat et al. 2016). The entire catalogue of long GRBs above 0.9 ph s −1 cm −2 is in light-grey, the bursts complying with the "good" criteria (i.e. small errors on the parameters) are shown in blue, the bursts for which Band is the best-fit spectral model are shown in orange. Right: reduced χ 2 distribution for Band models from the GBM spectral catalogue (Bhat et al. 2016).
its large sample size (N ∼ 7000 once the efficiency correction is taken into account). In order to strengthen the impact of the redshift constraint (weak due to its small sample size, N = 82), we added a weight of 10 to its log-likelihood. We tested different weights and this one was chosen as a balance between having a notable impact and being unrealistically constraining with respect to the other constraints.

Appendix 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 user-defined rate. As the effective temperature τ decreases, the chain is less and less likely to accept worse jumps until τ → 1, where we return to a classic Metropolis-Hastings algorithm. This modification to the original Metropolis-Hastings algorithm verifies the condition of ergodicity: regardless of the starting point the Markov Chain will converge to the same stationary distribution.
We also tested the methodological soundness of our algorithm by applying it to data generated by known inputs. We created mock observations for the intensity, spectrum and redshift constraint which we then used as constraints for our MCMC scheme; the results of the parameter space exploration are presented in Fig. B.1. We can see the code is able to satisfactorily recover the input parameter values: Schechter luminosity function {L * = 10 53 erg s −1 , p = 1.5}, redshift distribution following the shape of the cosmic star formation rate density; intrinsic E p − L correlation scenario {E p0 = 10 2.60 keV, σ E p = 0.30, α A = 0.50}.

Appendix C: Fluence cut -sample size plane
The sample size and fluence cut of SHOALS do not allow us to discriminate between non-evolving (k evol = 0) and strongly evolving (k evol = 2) scenarios for the luminosity function of LGRBs. In this section we detail our methodology for determine 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 F 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 p-values and D-statistics shown in the bottom panels of Fig. C.1. We define the probability of discriminating between the two scenarios at the 95% confidence level as the fraction of realisations with a p-value below 0.05; in the example of Fig. C.1, this probability would be 12%. We calculate this fraction for different values of N sub and F cut and draw the contours to obtain the fluence cut -sample size plane shown in Fig. 8.

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.

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, SFR3 respectively) and two scenarios for the peak energy distribution (intrinsically correlated with the peak luminosity, labelled A for "Amati-like", and independent log-normal, labelled LN). We used a Band function for the LGRB spectra with their empirical distributions for α and β. We used the same cosmology: Ω m = 0.3, Ω Λ = 0.7, and H 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.  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. The top panels correspond to the intrinsic E p − L correlation scenarios and the bottom panels to the independent log-Normal E p scenarios.

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 gamma-ray peak flux is reasonable.

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, σ E p = 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 on a much much longer duration.

Pescalli et al. 2016
Pescalli et al. (2016) derived step functions for the LGRB formation rate and luminosity function, using a different approach than the previous studies cited above, based on the C − method of Lynden-Bell (1971). In this work, contrary to Salvaterra et al. (2012), the authors constraints rely only on the luminosity and redshift distributions (and their join distribution: the luminosity-redshift plane) of the extended BAT6 sample of bright Swift LGRBs with N pk 15−150 keV ≥ 2.6 ph s −1 cm −2 . We use the values they quote in Section 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, σ E p = 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 Figure 1 of Pescalli et al. (2016). This is probably due to the difference in the low-luminosity slope between the two cases, which affects the number of faint bursts and thus the faint end of the log N-log P. The discrepancies in the log N-log P underline the advantage of constraining a population model on multiple complementary observational constraints that reach down to faint peak fluxes.