© The Authors 2026

1. Introduction

Gamma-ray bursts (GRBs) are high energy gamma-ray transients of astrophysical origin. They are thought to be produced as a result of cataclysmic astrophysical events, and they are commonly classified according to their time duration and spectral hardness (Kouveliotou et al. 1993; von Kienlin et al. 2020): Short GRBs (SGRBs) are shorter in time duration and on average harder in energy spectrum compared to long GRBs (LGRBs). While a sizeable number of LGRBs have been solidly associated with extreme cases of core-collapse supernovae (CCSNe), SGRBs are believed to mainly be associated with compact binary coalescence (CBC) events, namely binary neutron star (BNS) or neutron star-black hole (NSBH) mergers (Berger 2014). This connection was cemented by the joint observation and unambiguous association of GRB 170817A and the gravitational wave signal GW170817 from a BNS merger (Abbott et al. 2017).

Studies of the SGRB population generally aim to study the intrinsic luminosity distribution (also known as the ‘luminosity function’) of those events, inferring its parameters from a sample of GRBs observed at cosmological distances by various instruments. While some studies describe the luminosity distribution as an empirically defined function (Wanderman & Piran 2015, hereafter W15; Ghirlanda et al. 2016, hereafter G16), others have tested the structured jet hypothesis by utilising more physically motivated models to describe the emitted luminosity (Tan & Yu 2020; Salafia et al. 2023, hereafter S23).

Bearing in mind the CBC origin scenario, the cosmological rate density evolution of the SGRB population depends on the time needed by the compact objects to form from massive stars and for the binary system to merge. Simple theoretical considerations for a BNS system whose orbital separation shrinks due to gravitational radiation suggest that time delays (τ_d) between the star formation and the merger should be distributed as a power-law dP/dτ_d ∝ τ_d^−α_τ with α_τ ∼ 1 (Piran 1992). This model has been widely used in SGRB population studies (e.g. Paul 2018, Tan & Yu 2020, W15). For example, W15, who studied a sample of SGRBs detected by the Fermi Gamma-ray Burst Monitor (Fermi/GBM), the Swift Burst Alert Telescope (Swift/BAT), and the Compton Gamma-Ray Observatory Burst And Transient Source Experiment (CGRO/BATSE), found a delay-time distribution (DTD) power-law index α_τ ∼ 1, keeping the minimum time delay fixed at τ_d^min = 20 Myr. Some works suggest, though, that the power-law index of the DTD might be steeper, leading to shorter coalescence times. D’Avanzo et al. (2014), for example, selected a sample of Swift/BAT using quality cuts that maximise the probability of a redshift measurement and found that its DTD is best represented by α_τ ∼ 1.5 and τ_d^min = 10 Myr. More recent studies of the association between SGRB events with their host galaxies (Zevin et al. 2022) also inferred a steeper power-law DTD, with α_τ ∼ 1.5 − 2.2, but with larger minimum time delays, τ_d^min ∼ 105 − 250 Myr. An analysis of the chemical enrichment of r-process heavy elements in the Milky Way, on the other hand, under the assumption that they are mostly generated in BNS mergers, finds that fast mergers are required, with α_τ ≳ 2.0 and τ_min ≲ 40 Myr (Chen et al. 2025).

Another widely adopted model of the SGRB DTD is a log-normal distribution $\mathrm{d}P/\mathrm{d}{\tau }_{\mathrm{d}}=exp[-1/2{(ln{\tau }_{\mathrm{d}}-ln{\mu }_{\tau })}^2/2{\sigma }_{\tau }^2)]/{\sigma }_{\tau }{\tau }_{\mathrm{d}}\sqrt{2\pi }$, for which W15 finds μ_τ ∼ 3 − 4 Gyr and σ_τ ≲ 0.2. Similar values are reported by Luo et al. (2022), where a simulated SGRB population with a structured jet model fitting with the observed Fermi/GBM photon flux distribution is shown to be compatible with a log-normal DTD, disfavouring power-law and Gaussian DTD models.

Some works find different evolution branches for the BNS population, which can correspond to different BNS formation channels. Studies of BNS systems in our Galaxy show that the BNS population might be composed of two distinct sub-populations: a ‘fast’ population, with either a power-law DTD with α_τ ∼ 2 (Maoz & Nakar 2025) or a log-normal one with μ_τ ∼ 300 Myr and σ_τ ∼ 1 (Beniamini & Piran 2019), and a ‘slow’ population, with τ_d typically above ∼1 Gyr and with α_τ ∼ 1. In particular, Maoz & Nakar (2025) estimates that within Galactic BNS systems that merge within a Hubble time, the fast component is about ten to 100 times more abundant than the slow one. Beniamini & Piran (2019) reached a similar but slightly weaker conclusion that almost half of the Milky Way BNS population must belong to the fast-merging channel. Multi-band optical and near-infrared observations of SGRB host galaxies aimed at measuring their stellar masses and population ages (Leibler & Berger 2010) have shown a similar duality, where long (τ_d ∼ 3 Gyr) and short time delays (τ_d ∼ 0.2 Gyr) are respectively associated with early and late-type galaxies.

In this paper we perform a study of the SGRB population, assuming two different models of the DTD and two distinct models of the luminosity function, to test the robustness of our conclusions. The parameters are inferred within a hierarchical Bayesian framework (Mandel et al. 2019, hereafter M19).

In Sect. 2 we describe our methodology and sample selection. Results are shown in Sect. 3 and discussed in Sect. 4, where we also demonstrate a few possible sources of biases, and we analyse their impact on inferring the parameters of our population.

2. Methods

2.1. Bayesian hierarchical inference method

Our approach to characterise the SGRB population closely follows S23, with some updates that we describe below. The parameters of the population model (formally called ‘hyper-parameters’) constitute the elements of the hyper-parameter vector λ′_pop. Each SGRB is described by a vector of source parameters, λ_src, i, which are estimated based on the measured data, d_i (here i = 1, ..., N_obs and N_obs is the number of events in the sample). Following Eqs. 7 and 8 in M19, the posterior probability distribution function (PDF) can be written as

$$\begin{array}{c}\hfill P(\mathit{\lambda }{\prime }_{\mathrm{pop}}|\{{\mathit{d}}_i\})=\frac{\pi (\mathit{\lambda }{\prime }_{\mathrm{pop}})P(\{{\mathit{d}}_i\}|\mathit{\lambda }{\prime }_{\mathrm{pop}})}{P(\{{\mathit{d}}_i\})}=\frac{\pi (\mathit{\lambda }{\prime }_{\mathrm{pop}})}{P(\{{\mathit{d}}_i\})}{\displaystyle \prod _{i=1}^{N_{\mathrm{obs}}}}\frac{{\mathcal{N}}_i({\mathit{d}}_i|\mathit{\lambda }{\prime }_{\mathrm{pop}})}{\mathcal{D}(\mathit{\lambda }{\prime }_{\mathrm{pop}})}.\end{array}$$(1)

Here π(λ_pop′) is the ‘hyper-prior’, that is, the joint prior on the hyper-parameters, and P({d_i}) is the Bayesian evidence of the data (in practice, a normalisation constant). The likelihood of the data given the hyper-parameters, P({d_i}|λ′_pop), is written as a product of single-source terms, each of which is expressed as a ratio of a numerator over a denominator:

$$\begin{array}{c}\hfill {\mathcal{N}}_i({\mathit{d}}_i|\mathit{\lambda }{\prime }_{\mathrm{pop}})={\displaystyle \int P}({\mathit{d}}_i|{\mathit{\lambda }}_{\mathrm{src}})P_{\mathrm{pop}}({\mathit{\lambda }}_{\mathrm{src}}|\mathit{\lambda }{\prime }_{\mathrm{pop}})\mathrm{d}{\mathit{\lambda }}_{\mathrm{src}},\end{array}$$(2)

and

$$\begin{array}{c}\hfill \mathcal{D}(\mathit{\lambda }{\prime }_{\mathrm{pop}})={\displaystyle \int P_{\det }}({\mathit{\lambda }}_{\mathrm{src}})P_{\mathrm{pop}}({\mathit{\lambda }}_{\mathrm{src}}|\mathit{\lambda }{\prime }_{\mathrm{pop}})\mathrm{d}{\mathit{\lambda }}_{\mathrm{src}}.\end{array}$$(3)

Here, P_pop(λ_src|λ′_pop) is the ‘population probability’ that specifies the probability density of the source parameters for any choice of the hyper-parameters, and P_det(λ_src) is the selection function. It represents the probability that an event with source parameters λ_src, sampled from the population, is selected and hence included in the sample. The selection function therefore contains all the information on selection effects that stem from both the detection process and from any additional selection cuts that define the sample. Further, P(d_i|λ_src) is the likelihood of measuring the data for a single event given the source parameters. This framework has been implemented in the publicly available grbpop code¹, which we used as a starting point to implement different SGRB population models.

2.2. Including the observed number of events in the inference

To take into account the information on the number of events in the sample and hence include the SGRB local rate density (R₀) as one of the parameters of the population, we updated the definition of our posterior probability by embedding the Poissonian probability of the number of observed events. Following Eq. (11) in M19, given the total rate of astrophysical events in a population,

$$\begin{array}{c}\hfill R({\mathit{\lambda }}_{\mathrm{pop}})={\displaystyle {\int }_0^{\infty }}\frac{\dot{\rho }(z,{\mathit{\lambda }}_{\mathrm{pop}})}{1+z}\frac{\mathrm{d}V}{\mathrm{d}z}\mathrm{d}z,\end{array}$$(4)

with λ_pop ≡ (λ′_pop, R₀), the number of expected detections for a sample of a given detector is

$$\begin{array}{c}\hfill N_{\det }({\mathit{\lambda }}_{\mathrm{pop}})={\eta }_{\mathrm{DC}}·T·R({\mathit{\lambda }}_{\mathrm{pop}})·\mathcal{D}(\mathit{\lambda }{\prime }_{\mathrm{pop}}),\end{array}$$(5)

where η_DC is the product of the detector duty cycle times the average fraction of the sky accessible to it,²R(λ_pop) is the total cosmic rate of events, T is the duration of the observation period, and 𝒟(λ′_pop) is the term defined in Eq. (3), which can be shown to be equivalent to the fraction of events in the population that pass the detection and selection cuts. The probability of observing N_obs events given N_det(λ_pop) is then a Poissonian probability, and the posterior PDF becomes

$$\begin{array}{c}\hfill P({\mathit{\lambda }}_{\mathrm{pop}}|\{{\mathit{d}}_i\})=\frac{\pi (\mathit{\lambda }{\prime }_{\mathrm{pop}})}{P(\{{\mathit{d}}_i\})}{\displaystyle \prod _{i=1}^{N_{\mathrm{obs}}}}\left[\frac{{\mathcal{N}}_i({\mathit{d}}_i|\mathit{\lambda }{\prime }_{\mathrm{pop}})}{\mathcal{D}(\mathit{\lambda }{\prime }_{\mathrm{pop}})}\right]e^{-N_{\det }}{(N_{\det })}^{N_{\mathrm{obs}}},\end{array}$$(6)

where the $\frac{1}{N_{\mathrm{obs}}!}$ term is omitted due to the distinguishability of the events in the sample (M19). We defined the isotropic equivalent of the peak luminosity, L(p_[E0, E1], E_p, z), following Eq. (17) in S23 and modelling the SGRB photon spectrum as a power-law with an exponential cut-off Ghirlanda et al. (2004), with the low-energy photon index set to the median value of α based on the spectral analysis of the Fermi/GBM SGRB catalogue (S23).

2.3. Population models

2.3.1. Rate density evolution

Guided by the expectation that short GRBs are generated mainly in compact binary mergers, such as BNS or NSBH systems, the redshift probability distribution was modelled considering a fixed cosmic star formation history (CSFH) and by convolving it with the distribution of time delays (τ_d) between the binary system formation and the eventual merger. We considered two DTD models: a power law with a minimum time delay cut-off, namely

$$\begin{array}{c}\hfill P({\tau }_{\mathrm{d}}|\mathit{\lambda }{\prime }_{\mathrm{pop}})\propto \{\begin{array}{cc}{\displaystyle 0,}\hfill & {\tau }_{\mathrm{d}}<{\tau }_{\mathrm{d}}^{\min };\hfill \\ {\displaystyle {\tau }_{\mathrm{d}}^{-{\alpha }_{\tau }},}\hfill & {\tau }_{\mathrm{d}}\ge {\tau }_{\mathrm{d}}^{\min },\hfill \end{array}\end{array}$$(7)

where the free parameters are τ_d^min and the index α_τ, and a log-normal DTD, that is

$$\begin{array}{c}\hfill P({\tau }_{\mathrm{d}}|\mathit{\lambda }{\prime }_{\mathrm{pop}})=exp[-\frac{1}{2}{\left(\frac{ln{\tau }_{\mathrm{d}}-ln{\mu }_{\tau }}{{\sigma }_{\tau }}\right)}^2]{({\tau }_{\mathrm{d}}\sqrt{2\pi {\sigma }_{\tau }^2})}^{-1},\end{array}$$(8)

where the parameters to constrain are the median value μ_τ and the dispersion σ_τ. The rate density distribution is therefore

$$\begin{array}{c}\hfill {\displaystyle \dot{\rho }(z,\mathit{\lambda }{\prime }_{\mathrm{pop}})\propto {\displaystyle {\int }_z^{\infty }}\psi (z\prime )P(t_{\mathrm{LB}}(z\prime )-t_{\mathrm{LB}}(z)|\mathit{\lambda }{\prime }_{\mathrm{pop}})\frac{\mathrm{d}t}{\mathrm{d}z\prime }\mathrm{d}z\prime ,}\end{array}$$(9)

where ψ(z) is the CSFH and t_LB(z) is the cosmological lookback time at redshift z. We chose to adopt the CSFH model from Madau & Fragos (2017),

$$\begin{array}{c}\hfill \psi (z)\propto \frac{{(1+z)}^{a_{\psi }}}{1+{\left(\frac{1+z}{1+z_{\psi }}\right)}^{b_{\psi }}},\end{array}$$(10)

with a_ψ = 2.6, b_ψ = 6.2, and z_ψ = 2.2. The redshift probability distribution of SGRBs is then computed as

$$\begin{array}{c}\hfill {\displaystyle P(z|\mathit{\lambda }{\prime }_{\mathrm{pop}})\propto \frac{\dot{\rho }(z,\mathit{\lambda }{\prime }_{\mathrm{pop}})}{1+z}\frac{\mathrm{d}V}{\mathrm{d}z},}\end{array}$$(11)

where dV/dz is the differential comoving volume at redshift z.

2.3.2. Luminosity function models: ELF and QUSJ

To describe the luminosity function (i.e. the luminosity PDF) of our population, we chose two different models. The first one, which we refer to as the ‘empirical luminosity function’ model (ELF hereafter), is based on a doubly broken power-law luminosity function, and it extends the W15 luminosity function model to lower luminosities. Following³ Abbott et al. (2022), we express the luminosity function in this model as

$$\begin{array}{c}\hfill {\displaystyle P(L|\mathit{\lambda }{\prime }_{\mathrm{pop}})\propto \{\begin{array}{cc}0\hfill & L<L_0;\hfill \\ {\displaystyle {\left(\frac{L_{\ast \ast }}{L_{\ast }}\right)}^{-{\alpha }_{\mathrm{BPL}}-1}{\left(\frac{L}{L_{\ast \ast }}\right)}^{-{\gamma }_{\mathrm{BPL}}-1},}\hfill & L_0\le L\le L_{\ast \ast };\hfill \\ {\displaystyle {\left(\frac{L}{L_{\ast }}\right)}^{-{\alpha }_{\mathrm{BPL}}-1},}\hfill & L_{\ast \ast }<L\le L_{\ast };\hfill \\ {\displaystyle {\left(\frac{L}{L_{\ast }}\right)}^{-{\beta }_{\mathrm{BPL}}-1},}\hfill & L>L_{\ast }.\hfill \end{array}}\end{array}$$(12)

Here L₀, L_**, and L_* are the minimum luminosity, low-luminosity break, and high-luminosity break, respectively. The power-law indices of the three branches are expressed through the γ_BPL, α_BPL, and β_BPL parameters to facilitate comparison with Abbott et al. (2022). As a way to incorporate some information from GRB 170817A in this model, we chose our prior on L₀ to ensure that L₀ < L_17A, where $L_{17\mathrm{A}}=1.2_{-0.4}^{+0.5}\times {10}^{47}\mathrm{erg}{\mathrm{s}}^{-1}$ is the peak luminosity of GRB 170817A (S23).

To fully specify the population properties, the PDF of the peak photon energy E_p of the GRB spectral energy distribution across the population must also be specified. Following S23, we modelled this as a log-normal distribution with a median value that depends on L:

$$\begin{array}{c}\hfill P(E_{\mathrm{p}}|L,\mathit{\lambda }{\prime }_{\mathrm{pop}})=\frac{exp[-\frac{1}{2}{\left(\frac{ln(E_{\mathrm{p}})-ln({\stackrel{\sim }{E}}_{\mathrm{p}}(L))}{{\sigma }_{\mathrm{c}}}\right)}^2]}{E_{\mathrm{p}}\sqrt{2\pi {\sigma }_{\mathrm{c}}^2}},\end{array}$$(13)

where ${\stackrel{\sim }{E}}_{\mathrm{p}}(L)=E_{\mathrm{p}}^{★}{(L/L_{\ast })}^y$ allows for a Yonetoku-like correlation (Yonetoku et al. 2004) when y ≠ 0. The source parameters L and E_p are therefore parametrized with a common probability density: P(L, E_p|λ′_pop) = P(E_p|L, λ′_pop)P(L|λ′_pop).

The second model we considered is the quasi-universal structured jet population model (QUSJ hereafter) from S23, which also features a correlated joint distribution of L and E_p, induced by the jet structure. The latter is defined by two dimensionless functions of the ratio of the viewing angle, θ_v, to the ‘jet core’ angle θ_c, written as ℓ(θ_v/θ_c) and η(θ_v/θ_c), that describe how the luminosity and E_p depend on the jet viewing angle (see S23 for more details). In the analysis based on the QUSJ model, we include the information on the luminosity, peak photon energy, and viewing angle of GRB 170817A in the form of priors on the population parameters, following S23.

In both cases the redshift probability distribution, P(z|λ′_pop), is assumed to be independent of L and E_p. The distribution of source parameters across the population is therefore

$$\begin{array}{c}\hfill P_{\mathrm{pop}}({\mathit{\lambda }}_{\mathrm{src}}|\mathit{\lambda }{\prime }_{\mathrm{pop}})=P(L,E_{\mathrm{p}}|\mathit{\lambda }{\prime }_{\mathrm{pop}})P(z|\mathit{\lambda }{\prime }_{\mathrm{pop}}).\end{array}$$(14)

2.3.3. W15 population model

To gain insight on the difference between our results and those of previous studies, we set out to reproduce the W15 constraints using hierarchical Bayesian inference but with the same sample and selection effects modelling used in that work. For this purpose, the luminosity probability distribution P(L|λ′_pop) we employed is the same as in Eq. (12), but we fixed L₀ = L_** = 5 × 10⁴⁹ erg s⁻¹. The parameters to be constrained for P(L|λ′_pop) therefore are only α_BPL, β_BPL, and L_*.

The redshift probability distribution, P(z|λ′_pop), was modelled as in Sect. 2.3.1. The only difference is for the power-law DTD model, where the minimum time delay was fixed to τ_d^min = 20 Myr. The SGRB photon spectrum dṄ/dE(E,_p) was modelled as a Band function Band et al. (1993), whose parameters were fixed to E_p = 800 keV, α_Band = −0.5, and β_Band = −2.25, as in the original paper (motivated by the analysis of Nava et al. 2011).

Since the spectral energy peak in this model is fixed, the source parameters are only the GRB luminosity and redshift, that is, λ_src = (L, z). In Appendix A we give technical details on how this difference is handled in the inference.

2.4. Samples and detection efficiencies

When building our reference SGRB samples, we needed to carefully model the selection effects that moulded them. We therefore considered two samples for the analysis, built in the same way as in the ‘flux-limited sample analysis’ of S23.

The likelihood numerators, 𝒩_i(d_i|λ′_pop), for each event in the samples and denominators, 𝒟(λ′_pop), were then computed in the same way as in S23 for both the ELF and the QUSJ case since the only term that changes between the two models is P_pop(λ_src|λ′_pop).

2.4.1. Observer-frame sample

The first sample consists of short GRBs detected by Fermi/GBM with publicly available spectral analysis results. We selected events with T₉₀ < 2 s from the GBM’s first 10 years of observation, namely from the 12 July 2008 to 11 July 2018. We applied a completeness cut on the 64 ms peak photon flux in the [50 − 300] keV energy band, selecting only the events with p_[50, 300] > p_lim, GBM = 3.50 cm⁻² s⁻¹, which is the value above which we can consider our sample as complete in flux (S23). Additional cuts on the observed spectral peak energy (E_p, obs) were applied, considering only events with best-fit values 50 keV ≤ E_p, obs ≤ 10 MeV, which is the spectral range where the effective area of the Fermi/GBM detectors is optimal. After all cuts, the sample contained 210 short GRB events. Hereon, we refer to this as the ‘observer-frame sample’.

The detection efficiency for this sample was modelled as

$$\begin{array}{cc}& P_{\det }(L,E_{\mathrm{p}},z)=\mathrm{\Theta }(p_{[50,300]}(L,E_{\mathrm{p}},z)-p_{\lim ,\mathrm{GBM}})\times \hfill \\ \hfill & \mathrm{\Theta }(E_{\mathrm{p}}/(1+z)-50\mathrm{keV})\mathrm{\Theta }(10\mathrm{MeV}-E_{\mathrm{p}}/(1+z)),\hfill \end{array}$$(15)

where Θ(x) is the Heaviside step function.

Given the number of events contained in the observer-frame sample, we considered its event detection rate to estimate the astrophysical local rate density (R₀) of the SGRB population. The N_det for the Poissonian term in Eq. (6) was therefore computed using the 𝒟(λ′_pop) for this sample, an observation period of T = 10 yr, and a multiplicative factor of η_DC = 0.59 accounting for both the accessible field of view and the duty cycle of Fermi/GBM (Burns et al. 2016).

2.4.2. Rest-frame sample

The second sample consist of short GRBs detected by Fermi/GBM and Swift/BAT that pass a set of stringent cuts that maximise the redshift completeness without distorting the redshift distribution. In practice, this is the sub-sample of the S-BAT4ext catalogue (Ferro et al. 2023) of events that were also detected by Fermi/GBM. The selection includes a Swift/BAT completeness flux threshold in the [15 − 150] keV energy band, p_[15, 150] > p_lim, BAT = 3.50 cm⁻² s⁻¹, in addition to the same Fermi/GBM completeness cuts as in the observer-frame sample. The sample consists of 18 short GRBs events, 16 of which have a measured redshift. For those events we used the results of the spectral analysis performed in S23, which yields, for each event i, a posterior on the source parameters P(L, E_p, z | d_i), given a prior π(L, E_p, z)∝(1 + z)⁻¹L⁻¹. We refer to this sample as the ‘rest-frame sample’. We did not use the number of events in this sample to inform our estimate of R₀ because of the difficulty in correctly modelling the duty cycle factor η_DC incorporating all the cuts that define the sample.

Because we required the detection by Fermi/GBM, the rest-frame sample detection efficiency includes the same cuts as the observer-frame sample. In addition, it features a factor

$$\begin{array}{c}\hfill \mathrm{\Theta }(p_{[15,150]}(L,E_{\mathrm{p}},z)-p_{\lim ,\mathrm{BAT}})\end{array}$$(16)

to account for the Swift/BAT completeness cut. All other S-BAT4 cuts are independent on the source properties and hence do not enter the definition of this efficiency.

We note that GRB 170817A is not included in either of the two samples because its peak flux is below the Fermi/GBM completeness threshold. In the case of the quasi-universal jet population model, our results still include information on this event in the form of a prior on the hyper-parameters, which we built following the procedure used by S23 in their flux-limited sample analysis.

2.4.3. W15 samples

The observer-frame samples considered for the W15 population analysis consist of short GRBs observed by Fermi/GBM, CGRO/BATSE, and Swift/BAT. For the three detectors, the detection efficiency was modelled as a hard threshold on the 64 ms peak photon flux in their respective most sensitive detector energy band, as in Eq. (16). The photon flux thresholds considered for CGRO/BATSE and Fermi/GBM in the [E₀,  E₁]=[50, 300] keV detector energy band for this analysis were, respectively, p_lim, BATSE = 1.5 cm⁻² s⁻¹ and p_lim, GBM = 2.37 cm⁻² s⁻¹, while for Swift/BAT it was p_lim, BAT = 2.5 cm⁻² s⁻¹ in the [E₀^d,  E₁^d]=[15,  150] keV band.

The BATSE and Fermi sample contain GRBs with a measured p_[E0, E1] above the respective detector threshold and T₉₀ < 2 s. While for BATSE these selection effects have been applied to the entirety of its catalogue, for a total of 414 events, the Fermi/GBM events were considered up to the 10 April 2013, as in W15, yielding 146 bursts.

The rest-frame sample was built following W15. The sample contains events detected by Swift/BAT with a measured redshift and L, with T₉₀ < 2 s, up to GRB 131004A. A cut on the probability of having a non-collapsar origin higher than 60%, estimated following Bromberg et al. (2013), was then applied, selecting a sub-sample of 12 events.

Since the spectral peak energy in this analysis is fixed to a punctual value in the source frame, the computation of the likelihood was a bit different from that in the ELF and QUSJ cases. Full details on the likelihood evaluation and the obtained results are described in Appendix A. In the following, we refer to the results of this analysis as W15*.

2.5. Priors and posterior evaluation

We chose broad priors on most parameters, considering either a uniform or a uniform-in-log distribution within a defined range for each parameter, with the following exceptions: (i) For θ_c and θ_w, we considered an isotropic prior (i.e. uniform in the subtended solid angle) π(θ_c/w)∝sin(θ_c/w). (ii) For the minimum luminosity of the ELF model, we chose the upper bound of the prior based on the luminosity of GRB 170817A, as stated previously. (iii) For the QUSJ model, we followed the approach of S23 to condition the prior on the observed properties of GRB 170817A (see their section 2.5.3). The bounds that we report in Table 1 refer to the priors before that conditioning. The population parameters were considered to be independent from each other, with the exceptions of (L₀, L_**, L_*) for the ELF model and (θ_c, θ_w) for the QUSJ model. In fact, for the triplet of luminosity breaks, we imposed L₀ < L_** < 10⁵¹ erg s⁻¹ < L_*, while for the break angles in the quasi-universal jet structure, we set θ_c < θ_w. The choice to put a boundary between L_** and L_* at L = 10⁵¹ erg s⁻¹ was motivated by wanting to avoid that those parameters have local peaks in the posteriors in the same parameter region. Common parameters in the two population models share common priors. The full set of priors for each parameter with their bounds are shown in Table 1.

The posterior PDF was estimated through dynamical nested sampling using the open source python package dynesty (Higson et al. 2019); dynesty. The number of initial live points was set to 100 times the dimensionality of the posterior PDF, for a total of 1200 and 1400 initial live points, respectively, for the ELF and the QUSJ models. As stopping criteria, we chose dlogz = 10⁻⁵ and 30 000 effective samples. Batches of live points were periodically and automatically added to reach those criteria, with each batch containing one fifth of the initial live points used for the sampling. The obtained posterior PDFs for all the cases considered are portrayed as corner plots in Appendix B, with contours representing the 90% and 50% confidence levels, while the estimated median values of the parameters with their 90% credible intervals are shown in Table 2.

3. Results

A summary of the constraints obtained on the population parameters for the ELF and QUSJ models, with both the power-law and log-normal DTDs, is given in Table 2. The results for the W15* analysis are reported in Appendix A, where we demonstrate a good agreement with the results in the original paper.

3.1. Time delay and redshift distribution

We illustrate the constraints obtained with our analyses on the shape of the rate density evolution with redshift in Fig. 1 (left-hand panels), where the $\dot{\rho }(z)/R_0$ curves (i.e. the SGRB rate density normalised to its value at z = 0) are compared with the normalized CSFH. The results are quite consistent between the ELF and QUSJ models, where the curves tend to closely follow the CSFH, while the W15* analysis yields a very different constraint. The DTD parameter constraints (right-hand panels) show that this is the result of much shorter time delays in our models with respect to those obtained in the W15* analysis. To make this more evident and to enable a comparison between different DTD models, we also show in the figure, in addition to the posterior probability on the DTD parameters, the distribution of the average time delay:

$$\begin{array}{c}\hfill ⟨{\tau }_{\mathrm{d}}⟩(\mathit{\lambda }{\prime }_{\mathrm{pop}})={\displaystyle \int {\tau }_{\mathrm{d}}}P({\tau }_{\mathrm{d}}|\mathit{\lambda }{\prime }_{\mathrm{pop}})\mathrm{d}{\tau }_{\mathrm{d}}.\end{array}$$(17)

Fig. 1. Rate density distributions as functions of redshift (left panels) and corner plots of the posterior PDFs for the DTD parameters and ⟨τd⟩ (right panels). Top and bottom panels show the results obtained with a power-law and a log-normal DTD, respectively. Rate density distributions are normalised to 1 at z = 0, and the CSFH from Madau & Fragos (2017) is shown in purple for comparison. Red and blue curves represent the distributions obtained considering, respectively, the ELF and the QUSJ models. W15* results are displayed in yellow.

In the power-law DTD scenario, both luminosity function models show a steep power-law index α_τ > 1.5, while the minimum time delays are below τ_d < 370 Myr (90% credible limits). The average time delays are generally below ∼1 Gyr, with a median and 90% credible range of $⟨{\tau }_{\mathrm{d}}⟩={77}_{-65}^{+619}\mathrm{Myr}$ (ELF) or $⟨{\tau }_{\mathrm{d}}⟩={105}_{-93}^{+665}\mathrm{Myr}$ (QUSJ). The time delays associated with the W15* DTD parameters, on the other hand, were found to be significantly larger, yielding $⟨{\tau }_{\mathrm{d}}⟩=4.2_{-2.8}^{+2.0}\mathrm{Gyr}$.

Similar results were obtained in the log-normal DTD case, where the median value of the log-normal distribution is μ_τ < 700 Myr, while σ_c is only loosely constrained to be somewhat below the upper end of its prior range, which is equivalent to 1.5 dex. For the ELF and QUSJ models, the average time delays are respectively $⟨{\tau }_{\mathrm{d}}⟩\sim {92}_{-79}^{+877}\mathrm{Myr}$ and $⟨{\tau }_{\mathrm{d}}⟩\sim {116}_{-102}^{+761}\mathrm{Myr}$, while the W15* analysis yields $⟨{\tau }_{\mathrm{d}}⟩=3.1_{-0.8}^{+0.7}\mathrm{Gyr}$.

3.2. Luminosity function

The red (resp. blue) filled areas in Figure 2 show the 90% credible bands of the luminosity function obtained assuming the ELF (resp. QUSJ) model for the power-law DTD (left-hand panel) and the log-normal DTD (right-hand panel). In each panel, we also show the corresponding result of the W15* analysis. The luminosity of GRB 170817A (S23) is shown with a vertical grey band for comparison.

Fig. 2. Luminosity probability distributions obtained considering either a power-law or a log-normal DTD (left and right panel, respectively). Results obtained with the broken power-law model are shown in red, and the ones obtained with the structured jet model are in blue. The shaded areas represent the 90% credible intervals. Luminosity functions from W15 are also shown in yellow. The measured luminosity for GRB 170817A is shown in grey as a reference point.

The two different DTD models considered have a limited impact on the shape of the luminosity function. Above L_iso ∼ 5 × 10⁵² erg s⁻¹, both of our luminosity function models show good agreement with each other and with the W15* results. Between L_iso ∼ 10⁴⁹ erg s⁻¹ and L_iso ∼ 10⁵² erg s⁻¹, the trend of the ELF tends to follow the one from W15*, while the structured jet model curves have a slightly shallower slope. This is a result of the choice of parametrization and a consequence of the larger amount of information on GRB 170817A included in the QUSJ model. Below L_iso ∼ 10⁴⁹ erg s⁻¹, the ELF model shows very large rate uncertainties compared to the structured jet scenario. This is again a consequence of the information from GRB 170817A included in that model, and in part it also reflects the stronger predictivity of the QUSJ model at low luminosity, where the minimum luminosity is set by the maximum possible viewing angle and by the slope of the apparent structure.

3.3. Local rate density

As shown in Fig. 3, the overall local rate density inferred using the ELF has a larger uncertainty compared to the result using the QUSJ model, regardless of the chosen DTD model. This is a direct consequence of the reduced predictivity of the ELF model on the minimum SGRB luminosity compared to the QUSJ model. It is instructive to compare these rate density constraints with the most recent local BNS merger rate density constraint inferred from gravitational wave observations, based on the fourth gravitational wave transient catalogue (GWTC-4; Abac et al. 2025, R₀ ∈ [7.6,  250] Gpc⁻³ yr⁻¹ as shown by the vertical grey band in the figure). For the ELF model with the power-law DTD and the QUSJ model with the log-normal DTD, most of the posterior support is within the GWTC-4 BNS rate density constraint. For the QUSJ model with the power-law DTD, the support is concentrated at somewhat higher local rate densities, but a substantial fraction of the probability is still within the gravitational wave constraint. The ELF model with the power-law DTD rate posterior peaks within the GWTC-4 BNS rate density boundaries, reaching almost a plateau up to the upper prior PDF bound. This is a result of the poor constraints on the L_** parameter in this case, which is correlated with log R₀ (see Fig. B.3 in Appendix B). In general, these results support a scenario where BNS mergers constitute the progenitors of the vast majority of SGRBs, in line with the conclusions of S23.

Fig. 3. Comparisons between local rate density probability distributions. Distributions obtained with the broken power-law and the structured jet luminosity models are shown in red and blue, respectively, while results corresponding to power-law and log-normal DTD are respectively displayed with dashed and dotted lines. The BNS rate density inferred from gravitational wave observations (Abac et al. 2025) is shown in grey (90% credible intervals).

4. Discussion

The estimated time delays based on all four combinations of the models considered range between ∼11 − 18 Myr and ∼764 − 854 Myr. The lower end of these distributions can be explained by fast BNS merger channels, which may happen in binary star systems forming with a short orbital separation and evolving quickly through a common envelope and a mass transfer phase (Tauris et al. 2017), where the merger time can further decrease thanks to neutron star natal kicks (Andrews & Zezas 2019; Beniamini & Piran 2024). The parameter ranges found for the power-law DTD are compatible with those estimated from the Milky Way r-process element abundances (Chen et al. 2025), suggesting that BNS mergers are responsible for a significant fraction of those elements in our Galaxy.

Nonetheless, our findings for both DTD models contrast with time delay estimations in other short GRB population studies, such as W15, Luo et al. (2022) and Tan & Yu (2020). The differences could in part be traced back to the different SGRB samples with a measured redshift considered. For example, the largest measured redshift for a GRB of the W15 sample is z = 1.131, while our sample contains four GRBs with a measured redshift beyond that value, up to z = 2.28 ± 0.14.

However, the different treatment of the selection effects could also play a role in the estimate of the SGRB time delays. On one hand, the chosen peak-flux thresholds for Fermi/GBM in W15 and Tan & Yu (2020) and for Swift/BAT in W15 are below the respective completeness flux thresholds. Moreover, many previous studies relied on rest-frame samples (i.e. samples of GRBs with measured redshifts) constructed without a careful assessment of their selection effects and how these could have distorted their redshift distribution. To test for the impact of incorrectly modelling the selection effects on inferring the time delay parameters, we performed the additional analysis runs described below.

To demonstrate the biases deriving from using a flux-incomplete sample, we simulated a sample of short GRBs with an ELF and a power-law DTD using the following parameters: α_BPL = 2, β_BPL = 3, γ_BPL = 1, L₀ = 10⁴⁶ erg s⁻¹, L_** = 10⁴⁸ erg s⁻¹, L_* = 10⁵² erg s⁻¹, τ_d^min = 0.1 Myr, α_τ = 1, E_p^★ = 1 MeV, y = 0.3, and σ_c = 1. We considered Fermi/GBM as the only detector, and we modelled its detection efficiency following S23. Figure 4 shows the assumed detection efficiency as a function of p_[50, 300] (black line; in the plot we assume E_p, obs = 100 keV, but the dependence on the latter parameter is weak; see S23). This allowed us to obtain a sample of 271 simulated SGRBs that were ‘detected’. We divided this sample into an observer-frame and a rest-frame sample of 241 and 30 events, respectively. For the inference runs, we considered a hard flux threshold detection efficiency model P_det(L, E_p, z) = Θ(p_[50, 300](L,E_p,z)−p_lim) for different values of p_lim, namely 3.5, 3.0, 2.5, and 1.5 cm⁻² s⁻¹ (the assumed detection efficiency models are shown by the thick coloured lines in Fig. 4). For each p_lim value considered, we ran the inference method only on the sub-sample of simulated SGRBs brighter than that threshold.

Fig. 4. Comparison between the detection efficiency in our simulated sample and that used in the inference. The solid black line shows the detection efficiency model pdet, GBM from S23 as a function of the SGRB photon flux, assuming Ep, obs = 100 keV. The thick coloured lines show the hard-threshold detection efficiency models assumed in our inference of the simulated sample described in Sect. 4, with different colours indicating different assumed photon flux thresholds, as in Figure 5. For each of these analyses, only the SGRBs above the assumed threshold were included in the inference.

Figure 5 shows the posterior distribution of the average time delay for each of these four inference runs. When p_lim = p_lim, GBM = 3.5 cm⁻² s⁻¹, the sample is ‘complete’ in flux (in practice, the detection efficiency is ≳80% above the chosen threshold), and the posterior distribution peak coincides with the true value of the simulated population ⟨τ_d^true⟩ = 2.78 Gyr. When the assumed p_lim is less than the value of p_lim, GBM, the inference method systematically overestimates the inferred time delays. This simple simulation shows that the long time delays inferred in previous SGRB population studies might be, at least in part, the result of a bias due to an incorrect modelling of the selection effects.

Fig. 5. Posterior distribution of the average time delays inferred from a simulated SGRB population while assuming different photon peak flux cuts. The dashed black line shows the ‘true’ average time delay value corresponding to the power-law DTD parameters from which the SGRB events have been sampled. The blue curve shows the result obtained with the same flux threshold cut as in our analysis, which ensures flux completeness of the sample. The yellow, green, and pink curves were obtained with lower threshold values, as reported in the legend.

5. Conclusions

We have performed a hierarchical Bayesian study of the short GRB population in order to determine the time delays between those events and the CSFH. To minimise observational biases, we selected both a flux-complete and a high redshift-completeness sample of SGRBs. We considered two different models for the delay time distribution, namely, a power-law with a minimum time delay and a log-normal distribution, and we modelled the luminosity probability distribution of our population either as a doubly broken power law or starting from a quasi-universal jet structure. We chose the CSFH from Madau & Fragos (2017).

The inferred parameters for a power-law delay time distribution model are compatible with the estimates from the r-process elements in our Galaxy (Chen et al. 2025), and they correspond to time delays ranging from ∼12 Myr to ∼686 − 770 Myr. The steep power-index (α_τ > 1.5) deviates from the ‘conventional’ value α_τ = 1 adopted for a binary system that evolves while losing energy solely through gravitational radiation. This suggests that compact binary systems can merge through faster channels, for example forming with a low initial orbital separation and undergoing through a common envelope and mass transfer phase (Tauris et al. 2017). In some cases, the time required by the binary to merge might further reduce when neutron star natal kicks occur (Andrews & Zezas 2019; Beniamini & Piran 2024). The same fast merger channel scenarios are compatible with our findings for a log-normal delay time distribution, where the expected time delays go from ∼13 − 14 Myr to ∼877 − 969 Myr.

For each of the delay time distributions considered, we found that the choice of luminosity function model has a limited impact on the time delay results, supporting the conclusion that our finding does not depend specifically on our parametrization of the luminosity function.

Finally, we tested the impact of ill-modelled selection effects on the inferred time delays of a SGRB population by repeating the study for a simulated sample and using different values for its flux cut. We found that the inferred time delays become shifted towards higher values when using a flux threshold lower than the completeness one. We therefore stress the importance of building samples whose selection effects are understood and correctly modelled when performing statistical inference for parameter estimation, in order to minimise the observational biases when studying a population of astrophysical events.

Acknowledgments

MP acknowledges the fundings from the Fonds de la Recherche Scientifique – FNRS, Belgium, under grant No. 4.4501. OS acknowledges support from the Italian National Institute for Astrophysics (INAF) through ‘Finanziamento per la ricerca fondamentale’ grant number 1.05.23.04.04, and also funding by the European Union-Next Generation EU, PRIN 2022 RFF M4C21.1 (202298J7KT – PEACE).

References

Appendix A: Likelihoods and results for W15*

A.1. Likelihood evaluation

The likelihood numerator described in Eq. 2 for the Fermi/GBM and CGRO/BATSE samples are computed in a similar way to what done in S23 but considering the photon spectral peak energy in the source frame fixed to E_p^* = 800 keV. We can, in fact, write the posterior probability on the source parameters such as P(λ_src|d_i) = P(d_i|λ_src)π(L)π(z)/P(d_i). Then, assuming that P(λ_src|d_i)∝π(z) and neglecting the uncertainty on the measured p_[E0, E1], we can write

$$\begin{array}{c}\hfill P({\mathit{d}}_i|{\mathit{\lambda }}_{\mathrm{src}})\propto \frac{\delta (L-{\widehat{L}}_i(z))}{\pi (L)}P({\mathit{d}}_i),\end{array}$$(A.1)

where ${\widehat{L}}_i(z)=L(p_i,E_{\mathrm{p}}^{\ast },z)$. We can compute the prior on L by assuming a uniform prior in p_[E0, E1] and applying a coordinate transform

$$\begin{array}{c}\hfill \pi (L)=\frac{\partial p_i}{\partial L}\pi (p_i)\propto \frac{p_i}{L}\end{array}$$(A.2)

so that the numerator for the likelihood terms in the product for each event of the sample we obtain is

$$\begin{array}{c}\hfill {\mathcal{N}}_i^{\mathrm{BATSE}-\mathrm{GBM}}({\mathit{d}}_i|{\mathit{\lambda }}_{\mathrm{pop}}^{\prime })=P({\mathit{d}}_i){\displaystyle {\int }_0^{\infty }}\frac{{\widehat{L}}_i(z)}{p_i}{P}_{\mathrm{pop}}({\widehat{L}}_i(z),z|{\mathit{\lambda }}_{\mathrm{pop}}^{\prime })\mathrm{d}z.\end{array}$$(A.3)

The likelihood numerator for the Swift/BAT sample is built in an analogous way. By neglecting the uncertainty on the measured luminosity and redshift, we have

$$\begin{array}{c}\hfill P({\mathit{d}}_i|{\mathit{\lambda }}_{\mathrm{src}})\propto \frac{\delta (L-L_i)\delta (z-z_i)}{\pi (L,z)}P({\mathit{d}}_i).\end{array}$$(A.4)

If we assume a prior uniform in ln(L) and z, we have that π(L, z)∝L⁻¹, and therefore the likelihood numerator for Swift/BAT events is

$$\begin{array}{c}\hfill {\mathcal{N}}_i^{\mathrm{BAT}}({\mathit{d}}_i|{\mathit{\lambda }}_{\mathrm{pop}}^{\prime })=P({\mathit{d}}_i)L_i{P}_{\mathrm{pop}}(L_i,z_i|{\mathit{\lambda }}_{\mathrm{pop}}^{\prime }).\end{array}$$(A.5)

For the events in the three samples, the likelihood denominator can be computed as in Eq. 3:

$$\begin{array}{c}\hfill \mathcal{D}({\mathit{\lambda }}_{\mathrm{pop}}^{\prime })={\displaystyle \int \mathrm{\Theta }}(p_{[E_0,E_1]}({\mathit{\lambda }}_{\mathrm{src}},E_{\mathrm{p}}^{\ast })-p_{\lim ,k})P_{\mathrm{pop}}({\mathit{\lambda }}_{\mathrm{src}}|{\mathit{\lambda }}_{\mathrm{pop}}^{\prime })\mathrm{d}{\mathit{\lambda }}_{\mathrm{src}}.\end{array}$$(A.6)

The astrophysical SGRB local rate density R₀ is estimated using the rate of observed SGRBs in the Fermi/GBM sample, with its corresponding 𝒟(λ′_pop) and considering an effective observing time of T ⋅ η_DC = 3.65 yr (W15).

A.2. Priors and posterior evaluation

The prior PDFs for the population parameters λ′_pop have been chosen independent between each other, i.e. π(λ′_pop) = π(α_BPL)π(β_BPL)π(L_*)π(α_t)π(R₀) and π(λ′_pop) = π(α_BPL)π(β_BPL)π(L_*)π(μ_t)π(σ_t)π(R₀) respectively for the power-law and the log-normal DTD models. All the priors have been chosen either uniform or uniform in logarithm, within the bounds that we can deduce from W15 plots. The set of priors used is listed in table A.1.

As in Sect. 2.5, the evaluation of the posterior PDF is performed through dynamical nested sampling. The number of initial live points is set to 200 times the dimensionality of the posterior PDF, for a total of 1000 and 1200 initial live points respectively for the power-law and the log-normal DTD models. As stopping criteria, we choose dlogz = 10⁻⁵ and 30000 effective samples and batches of live points were periodically and automatically added to reach those criteria, with each batch containing one fifth of the initial live points.

A.3. Results

The obtained posterior PDF corner plots are shown in Fig. B.1 and Table A.2 shows the values estimated for the parameters of the population. Although the CSFH we used to reproduce those results (Madau & Fragos 2017) is not featured in the W15 study, its trend is very similar to the one from Planck Collaboration XXX (2014), which is one of the models used for the analysis in W15 and referred to as SFR2. We therefore take into account this set of results from W15 to compare ours. The uncertainty taken into account on the median values is 1σ as in W15, in order to make a direct comparison between the two works.

The results are quite similar to those of W15, except for a shorter median time delay and a larger variance for the log-normal DTD and the local rate density whose estimate is one order of magnitude higher than the one in W15. We attribute the difference in the estimate of R₀ to the different methods used, as in W15 they build their estimate on the fraction of non-collapsar events based on the criteria by Bromberg et al. (2013).

If we perform the same analysis with a wider range of parameters for our prior PDFs, we obtain considerably different results. In fact, if we consider a uniform prior in the interval [ − 3,  5] for α_BPL and a uniform-in-log prior for L_* between [10⁵⁰,  10⁵⁴] erg s⁻¹, the posterior starts to exhibit multiple local maxima as shown in Fig. B.1.

By looking at the luminosity functions obtained with these extended prior bounds (Fig. A.1), we can see that both in the power-law and the log-normal case they peak around L_iso ∼ 2 × 10⁵⁰ erg s⁻¹, following a single power-law trend after.

Fig. A.1. Luminosity functions obtained with the sample and model from W15, obtained considering a power-law and a log-normal DTD (respectively left and right panel). The yellow dashed curves correspond to the original parameter range considered for the study, while the green continuous curves represent the results obtained with a wider parameter range. Shaded areas depict the 90% percentiles.

The redshift distributions obtained, displayed in Fig. A.2, present different behaviours for the two delay time distribution models. The power-law model with fixed minimum time delay brings to a considerable uncertainty in the peak of the distribution, ranging from z ∼ 0 to z ∼ 2 when considering a 90% credible interval. Overall the curves look a bit shifted towards lower redshift when compared to the ones obtained with the tighter parameter range. On the other hand, the curves obtained with the log-normal model and a wider parameter range have a similar behaviour to what has been found in W15, but with a larger uncertainty on the distribution. The larger credible intervals are related to the wider log σ_τ values found with the choice of a less restrictive prior for the other parameters (see Fig. B.1).

Fig. A.2. Redshift distribution of SGRBs for the W15* population normalised to 1 at z = 0. The power-law and log-normal DTD scenarios are shown, respectively, in the left and right panel. The yellow dashed curves depict the original parameter range considered for the study, while the green continuous curves show the results obtained with a wider parameter range. Shaded areas correspond to the 90% percentiles.

Appendix B: Posterior PDF corner plots

Fig. B.1. Posterior PDF contour plots for the W15* results along with the marginalised posterior PDFs on the sides for each parameter, in the power-law DTD (left panel) and log-normal DTD (right panel) cases. Red and blue contours respectively represent the results obtained using the "original" and extended priors on the parameters. Darker and lighter shades of the contours depict, respectively, the 50% and 90% confidence levels. The error bars on the marginalised posteriors mark the median of the parameters along with their 1 σ uncertainty.

Fig. B.2. Corner plot of the posterior PDF for the ELF + power-law DTD model, with the relative marginalised posteriors for each parameter. Contours depict the 50% and 90% confidence levels. Vertical dashed lines on the marginalised posteriors represent the median values with the 90% credible intervals.

Fig. B.3. Corner plot of the posterior PDF for the ELF + log-normal DTD model, with the relative marginalised posteriors for each parameter. Contours depict the 50% and 90% confidence levels. Vertical dashed lines on the marginalised posteriors represent the median values with the 90% credible intervals.

Fig. B.4. Corner plot of the posterior PDF for the QUSJ + power-law DTD model, with the relative marginalised posteriors for each parameter. Contours depict the 50% and 90% confidence levels. Vertical dashed lines on the marginalised posteriors represent the median values with the 90% credible intervals.

Fig. B.5. Corner plot of the posterior PDF for the QUSJ + log-normal DTD model, with the relative marginalised posteriors for each parameter. Contours depict the 50% and 90% confidence levels. Vertical dashed lines on the marginalised posteriors represent the median values with the 90% credible intervals.

All Tables

All Figures

Fig. 1. Rate density distributions as functions of redshift (left panels) and corner plots of the posterior PDFs for the DTD parameters and ⟨τd⟩ (right panels). Top and bottom panels show the results obtained with a power-law and a log-normal DTD, respectively. Rate density distributions are normalised to 1 at z = 0, and the CSFH from Madau & Fragos (2017) is shown in purple for comparison. Red and blue curves represent the distributions obtained considering, respectively, the ELF and the QUSJ models. W15* results are displayed in yellow.

In the text

Fig. 2. Luminosity probability distributions obtained considering either a power-law or a log-normal DTD (left and right panel, respectively). Results obtained with the broken power-law model are shown in red, and the ones obtained with the structured jet model are in blue. The shaded areas represent the 90% credible intervals. Luminosity functions from W15 are also shown in yellow. The measured luminosity for GRB 170817A is shown in grey as a reference point.

In the text

Fig. 3. Comparisons between local rate density probability distributions. Distributions obtained with the broken power-law and the structured jet luminosity models are shown in red and blue, respectively, while results corresponding to power-law and log-normal DTD are respectively displayed with dashed and dotted lines. The BNS rate density inferred from gravitational wave observations (Abac et al. 2025) is shown in grey (90% credible intervals).

In the text

Fig. 4. Comparison between the detection efficiency in our simulated sample and that used in the inference. The solid black line shows the detection efficiency model pdet, GBM from S23 as a function of the SGRB photon flux, assuming Ep, obs = 100 keV. The thick coloured lines show the hard-threshold detection efficiency models assumed in our inference of the simulated sample described in Sect. 4, with different colours indicating different assumed photon flux thresholds, as in Figure 5. For each of these analyses, only the SGRBs above the assumed threshold were included in the inference.

In the text

Fig. 5. Posterior distribution of the average time delays inferred from a simulated SGRB population while assuming different photon peak flux cuts. The dashed black line shows the ‘true’ average time delay value corresponding to the power-law DTD parameters from which the SGRB events have been sampled. The blue curve shows the result obtained with the same flux threshold cut as in our analysis, which ensures flux completeness of the sample. The yellow, green, and pink curves were obtained with lower threshold values, as reported in the legend.

In the text

	Fig. A.1. Luminosity functions obtained with the sample and model from W15, obtained considering a power-law and a log-normal DTD (respectively left and right panel). The yellow dashed curves correspond to the original parameter range considered for the study, while the green continuous curves represent the results obtained with a wider parameter range. Shaded areas depict the 90% percentiles.
In the text

Fig. A.2. Redshift distribution of SGRBs for the W15* population normalised to 1 at z = 0. The power-law and log-normal DTD scenarios are shown, respectively, in the left and right panel. The yellow dashed curves depict the original parameter range considered for the study, while the green continuous curves show the results obtained with a wider parameter range. Shaded areas correspond to the 90% percentiles.

In the text

Fig. B.1. Posterior PDF contour plots for the W15* results along with the marginalised posterior PDFs on the sides for each parameter, in the power-law DTD (left panel) and log-normal DTD (right panel) cases. Red and blue contours respectively represent the results obtained using the "original" and extended priors on the parameters. Darker and lighter shades of the contours depict, respectively, the 50% and 90% confidence levels. The error bars on the marginalised posteriors mark the median of the parameters along with their 1 σ uncertainty.

In the text

	Fig. B.2. Corner plot of the posterior PDF for the ELF + power-law DTD model, with the relative marginalised posteriors for each parameter. Contours depict the 50% and 90% confidence levels. Vertical dashed lines on the marginalised posteriors represent the median values with the 90% credible intervals.
In the text

	Fig. B.3. Corner plot of the posterior PDF for the ELF + log-normal DTD model, with the relative marginalised posteriors for each parameter. Contours depict the 50% and 90% confidence levels. Vertical dashed lines on the marginalised posteriors represent the median values with the 90% credible intervals.
In the text

	Fig. B.4. Corner plot of the posterior PDF for the QUSJ + power-law DTD model, with the relative marginalised posteriors for each parameter. Contours depict the 50% and 90% confidence levels. Vertical dashed lines on the marginalised posteriors represent the median values with the 90% credible intervals.
In the text

	Fig. B.5. Corner plot of the posterior PDF for the QUSJ + log-normal DTD model, with the relative marginalised posteriors for each parameter. Contours depict the 50% and 90% confidence levels. Vertical dashed lines on the marginalised posteriors represent the median values with the 90% credible intervals.
In the text