© The Authors 2023

1. Introduction

Two main clusters of gamma-ray bursts (GRBs) have long been identified in the two-dimensional plane of duration versus hardness ratio¹ of the large sample collected by the Burst Alert and Transient Source Experiment (BATSE) on board the Compton Gamma-Ray Observatory (Kouveliotou et al. 1993). The bimodality in the BATSE GRBs is also apparent when considering the durations only (more precisely, the time, T₉₀, over which 5% to 95% of the background-subtracted counts are collected), with the histogram of the logarithms of the durations featuring two peaks separated by a valley at around T₉₀ = 2 s. This has since become the customary separation between the ‘long’ (LGRB) and the ‘short’ (SGRB) events in the GRB class, although the actual position of the valley varies somewhat across different detectors (e.g. Bromberg et al. 2013). The evidence accumulated during the following decades painted a picture of two different progenitor systems: starting with GRB 980425 (Galama et al. 1998), a number of LGRBs have been firmly associated with type Ib/c core-collapse supernovae (e.g. Bloom et al. 2002; Malesani et al. 2004; Mirabal et al. 2006; Kann et al. 2011; Cano et al. 2014; D’Elia et al. 2018; Hu et al. 2021), securing the scenario of a massive star progenitor (Woosley 1993); the progenitors of SGRBs remained elusive for a longer time, though all hints consistently pointed (e.g. Nakar 2007; Fong & Berger 2013; Berger 2014; D’Avanzo 2015) to a compact binary merger progenitor (Eichler et al. 1989; Mochkovitch et al. 1993). Such a progenitor has been confirmed by the astounding association (Abbott et al. 2017a,b) of the SGRB 170817A (Goldstein et al. 2017) – detected by the Gamma-ray Burst Monitor (GBM; Meegan et al. 2009) on board the Fermi spacecraft and by the International Gamma-Ray Astrophysics Laboratory (INTEGRAL; Savchenko et al. 2017) – with the first-ever binary neutron star (BNS) merger detected by humankind, GW170817 (Abbott et al. 2019), whose gravitational wave (GW) signal was captured on 17 August 2017 by the Advanced Laser Interferometer Gravitational wave Observatory (aLIGO; Aasi et al. 2015) and localised thanks to Advanced Virgo (Acernese et al. 2014).

The properties of SGRBs – especially their shorter duration and harder spectrum with respect to LGRBs (Ghirlanda et al. 2009; Calderone et al. 2015) – make them very hard to detect with current facilities: of the more than 280 GRBs revealed by Fermi/GBM every year, only about 40 are SGRBs. The softer sensitivity band of the Burst Alert Telescope (BAT) on board the Neil GehrelsSwift Observatory (Gehrels et al. 2004), together with its smaller field of view, allows it to identify and localise only a handful of SGRBs per year. Moreover, even when localised by Swift, the fraction of SGRBs that end up with a secure redshift determination is relatively low, due to a combination of a fainter X-ray ‘afterglow’ (Costa et al. 1997) – whose detection is a requirement for a precise localisation – and a larger typical offset from the host galaxy (Fong et al. 2013, 2022; Fong & Berger 2013; D’Avanzo et al. 2014; Berger 2014) with respect to LGRBs, which renders the host galaxy identification ambiguous in cases where multiple galaxies stand at similar offsets from the position of the afterglow. As a result, the intrinsic properties of the SGRB population are much more uncertain than for LGRBs. Indeed, a variety of attempts at constraining the properties of the SGRB population throughout the years, in some cases with very different methodologies and reference samples, has yielded varied results (e.g. Schmidt 2001; Guetta & Piran 2005, 2006; Virgili et al. 2011; Yonetoku et al. 2014; D’Avanzo et al. 2014; Wanderman & Piran 2015; Shahmoradi & Nemiroff 2015; Ghirlanda et al. 2016; Zhang & Wang 2018; Paul 2018; Tan & Yu 2020). Some of these results are in clear tension with each other: in this work, we consider Wanderman & Piran (2015; hereafter W15) and Ghirlanda et al. (2016; hereafter G16), which reach very different conclusions, as our benchmarks.

A prominent challenge in unveiling the intrinsic properties of SGRBs is the uncertainty about which processes shape the luminosity function, that is, the probability distribution from which the luminosity of each event in the population is sampled. Lipunov et al. (2001) and Rossi et al. (2002) were the first to realise that the highly relativistic nature of GRB jets would make their angular structure an important factor in determining the luminosity function, in addition to the intrinsic spread in luminosities. Indeed, if the energy density and the typical Lorentz factor of a GRB jet are functions of the angular separation from the jet axis, then the apparent energetics are viewing-angle-dependent by virtue of relativistic aberration effects (e.g. Salafia et al. 2015). Since jets are isotropically oriented in space, this naturally produces a large spread in the apparent luminosities, with a well-defined dependence on the angular structure (Pescalli et al. 2015). In the presence of a narrow distribution of intrinsic jet luminosities, the luminosity function is then mainly shaped by the angular structure (e.g. Salafia et al. 2020; Tan & Yu 2020).

An angular profile in the jet properties arises naturally in essentially any physically viable jet formation scenario (see Salafia & Ghirlanda 2022 for a recent review). Moreover, a non-trivial jet energy density angular profile is required (e.g. Mooley et al. 2018; Lamb et al. 2018, 2019; Ghirlanda et al. 2019; Takahashi & Ioka 2020, 2021; Beniamini et al. 2022) to explain the observed properties of the non-thermal afterglow of the SGRB associated with GW170817. Hence, the question of whether the observed distribution of SGRB properties can be traced back, at least in part, to the differing viewing angles is of particular relevance, as it would provide a route to a unification of these sources and a way to disentangle the intrinsic diversity in their properties (and hence those of their progenitor) from the apparent diversity due to extrinsic factors, particularly the viewing angle.

In the future, joint GW-GRB observations will provide direct information on the structure of SGRB jets thanks to the measurements of the inclination of the merging binary’s angular momentum (which is most likely a proxy of the jet viewing angle) that can be inferred from GW analysis (Williams et al. 2018; Hayes et al. 2020; Farah et al. 2020; Biscoveanu et al. 2020). Still, such information must be combined self-consistently with that encoded in the SGRB population observed in gamma-rays only. In this work, we describe our investigation of the population properties of SGRBs within a ‘quasi-universal jet’ scenario. We assumed that all SGRB jets share the same angular profile of luminosity as a function of the viewing angle, while we allowed for a spread in the on-axis luminosities (as in e.g. Tan & Yu 2020; Hayes et al. 2023). This produces a particular parametrisation of the SGRB population properties, which we derive and describe in detail in Sects. 2.2 and 2.3. In order to constrain the parameters of this model, we considered the sample of SGRBs detected by Fermi and carefully modelled the underlying selection effects (Sect. 2.4). This allowed us to fit the model to the data through a hierarchical Bayesian approach (Sect. 2.5). In Sect. 3 we describe in detail the results of the fit, and in Sect. 4 we discuss several implications.

Throughout this work, we assume a flat Friedmann-Lemaître-Robertson-Walker cosmology with Planck Collaboration XIII (2016) parameters, that is, H₀ = 67.74 km Mpc⁻¹ s⁻¹ and Ω_m, 0 = 0.3075.

2. Methodology

2.1. Apparent versus intrinsic structure

The dominant form of energy in GRB jets and the processes that dissipate such energy, leading to the observed ‘prompt’ gamma-ray emission are still a matter of debate (see Kumar & Zhang 2015, for a recent review). Still, regardless of the particular dissipation and emission process, the observed emission is affected by relativistic aberration effects in a way that depends on the Lorentz factor profile and the viewing angle (e.g. Woods & Loeb 1999; Salafia et al. 2015). For example, an observer looking, from a viewing angle θ_v (angle between the line of sight and the jet axis), at an axisymmetric jet with a bulk Lorentz factor profile Γ(θ) (where θ is the angle from the jet axis) that radiates an energy per unit solid angle dE/dΩ(θ), would measure an isotropic-equivalent gamma-ray energy (Salafia et al. 2015),

$$\begin{array}{c}\hfill E_{\mathrm{iso}}({\theta }_{\mathrm{v}})={\displaystyle {\int }_0^{\pi /2}{\int }_0^{2\pi }}\frac{{\delta }_{\mathrm{D}}^3(\theta ,\varphi ,{\theta }_{\mathrm{v}})}{\mathrm{\Gamma }(\theta )}\frac{\mathrm{d}E}{\mathrm{d}\mathrm{\Omega }}(\theta )\mathrm{d}\varphi sin\theta \mathrm{d}\theta ,\end{array}$$(1)

where ϕ is the azimuthal angle of a spherical coordinate system whose z-axis coincides with the jet axis, ${\delta }_{\mathrm{D}}(\theta ,\varphi ,{\theta }_{\mathrm{v}})={\mathrm{\Gamma }}^{-1}(\theta ){[1-\mathit{\beta }(\theta ,\varphi )·{\mathit{e}}_{\mathrm{LoS}}({\theta }_{\mathrm{v}})]}^{-1}$ is the Doppler factor, β(θ, ϕ) is the local jet bulk velocity vector – with a magnitude $\beta (\theta )=\sqrt{1-1/{\mathrm{\Gamma }}^2(\theta )}$ – and e_LoS is a unit vector pointing to the observer. Hence, when considering the emitted energy in gamma-rays, the ‘apparent structure’ E_iso(θ_v) depends on the intrinsic structure (Γ(θ),dE/dΩ(θ)) through a functional that is not invertible in general. The situation is even more nuanced when considering the luminosity, as the effective angular profile L_iso(θ_v) depends also on the degree of overlap between different pulses (Salafia et al. 2016), and hence on the intrinsic variability.

For these reasons, in the absence of strong theoretical constraints on the intrinsic jet structure and on the dissipation and emission processes, the most straightforward approach – which we adopt in this work – is that of parameterising directly the apparent luminosity structure L(θ_v) (we drop the ‘iso’ suffix from here on for simplicity), which also reduces the number of parameters. An assessment of the intrinsic jet structures that are compatible with a given apparent luminosity profile can be then carried out a posteriori.

2.2. Statistical model of an SGRB population with a quasi-universal jet

Within our framework, we describe each SGRB by four physical quantities, namely its viewing angle θ_v, its peak isotropic-equivalent luminosity L, the photon energy E_p at the peak of the spectral energy distribution (SED; i.e. the νF_ν spectrum) and its redshift z. These collectively represent what we refer to as the source parameter vector, ${\lambda }_{\text{src}}^{\star }=({\theta }_{\text{v}},L,E_{\text{p}},z)$. For most events, the viewing angle is unknown and we therefore consider a reduced source parameter vector λ_src = (L, E_p, z). We assumed the diversity in these parameters to be the combined result of intrinsic heterogeneity in the physical properties of jets within the population and extrinsic diversity induced by the differing viewing angles under which these jets are observed.

In order to represent the intrinsic heterogeneity in SGRB jets, we opted for parameterising the joint probability density P(L_c, E_p, c) of their on-axis (‘core’) peak luminosity L_c and peak SED photon energy E_p, c as follows. We assumed L_c to be distributed as a power law with index −A, with a lower exponential cutoff below $L_{\text{c}}^{\star }$, namely

$$\begin{array}{c}\hfill P(L_{\mathrm{c}}|A,L_{\mathrm{c}}^{\star })=\frac{A}{\stackrel{\sim }{\mathrm{\Gamma }}(1-1/A)L_{\mathrm{c}}^{\star }}exp[-{\left(\frac{L_{\mathrm{c}}^{\star }}{L_{\mathrm{c}}}\right)}^A]{\left(\frac{L_{\mathrm{c}}}{L_{\mathrm{c}}^{\star }}\right)}^{-A},\end{array}$$(2)

where $\stackrel{\sim }{\mathrm{\Gamma }}$ indicates a Gamma function, and the integrated probability is normalised to unity. This probability density is defined in such a way that it peaks at $L_{\text{c}}^{\star }$ regardless of the value of A. This choice of parametrisation follows from the fact that, in the quasi-universal jet scenario, the high end of the luminosity function reflects the luminosity distribution of jets seen close to on-axis (see Sect. 2.3), and the high-luminosity end has been found to be well described by a power law in most previous studies of the SGRB population (e.g. Schmidt 2001; Guetta & Piran 2006; Virgili et al. 2011; Yonetoku et al. 2014; D’Avanzo et al. 2014; Wanderman & Piran 2015; Ghirlanda et al. 2016; Zhang & Wang 2018; Paul 2018; Tan & Yu 2020). The reason for introducing a low-end cut-off follows from the quasi-universality assumption (in which case low luminosities are due to off-axis viewing angles rather than to intrinsically weak jets) and can be physically linked to a minimum jet luminosity required for the jet to break out from the progenitor vestige (see Salafia & Ghirlanda 2022 for a recent review).

The probability distribution on E_p, c, conditional on L_c, was assumed log-normal and centred at a L_c-dependent value ${\stackrel{\sim }{E}}_{\mathrm{p},\mathrm{c}}=E_{\mathrm{p},\mathrm{c}}^{\star }{(L_{\mathrm{c}}/L_{\mathrm{c}}^{\star })}^{\mathit{y}}$, where y sets the slope of the relation. Hence,

$$\begin{array}{c}\hfill P(E_{\mathrm{p},\mathrm{c}}|L_{\mathrm{c}},E_{\mathrm{p},\mathrm{c}}^{\star },y,{\sigma }_{\mathrm{c}})=\frac{exp[-\frac{1}{2}{\left(\frac{ln(E_{\mathrm{p},\mathrm{c}})-ln({\stackrel{\sim }{E}}_{\mathrm{p},\mathrm{c}})}{{\sigma }_{\mathrm{c}}}\right)}^2]}{E_{\mathrm{p},\mathrm{c}}\sqrt{2\pi {\sigma }_{\mathrm{c}}^2}},\end{array}$$(3)

where σ_c sets the dispersion of E_p, c around ${\stackrel{\sim }{E}}_{\mathrm{p},\mathrm{c}}$. This assumption allows for a ‘Yonetoku’ correlation² with slope y between the logarithms of the on-axis peak SED photon energy and the luminosity, which may be induced, for example, by the underlying emission process. The log-normal form of the scatter around the relation was chosen for its simplicity. The case with no correlation (hence with log-normally distributed values of E_p, c, un-correlated with L_c) is represented by y = 0 and it is therefore naturally included. The joint probability density of the core quantities is

$$\begin{array}{c}\hfill P(L_{\mathrm{c}},E_{\mathrm{p},\mathrm{c}}|{\mathit{\lambda }}_{\mathrm{pop}})=P(E_{\mathrm{p},\mathrm{c}}|L_{\mathrm{c}},{\mathit{\lambda }}_{\mathrm{pop}})P(L_{\mathrm{c}}|{\mathit{\lambda }}_{\mathrm{pop}}),\end{array}$$(4)

where the population parameter vector λ_pop contains $L_{\text{c}}^{\star }$, A, $E_{\text{p},\text{c}}^{\star }$, y, σ_c and all other parameters that fully specify the SGRB population model.

The next, key assumption of the model is that all jets share a universal ‘structure’ that specifies the dependence of L and E_p on the viewing angle θ_v. In practice, we assumed the viewing-angle-dependent luminosity and SED peak photon energy to be expressed as

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill L({\theta }_{\mathrm{v}},{\mathit{\lambda }}_{\mathrm{pop}})=& L_{\mathrm{c}}\ell ({\theta }_{\mathrm{v}},{\mathit{\lambda }}_{\mathrm{pop}})\text{and}\hfill \\ \hfill E_{\mathrm{p}}({\theta }_{\mathrm{v}},{\mathit{\lambda }}_{\mathrm{pop}})=& E_{\mathrm{p},\mathrm{c}}\eta ({\theta }_{\mathrm{v}},{\mathit{\lambda }}_{\mathrm{pop}}),\hfill \end{array}\end{array}$$(5)

where ℓ and η are functions of the viewing angle and of some parameters included in the λ_pop vector. These functions, which we assumed to be redshift-independent, define the universal apparent structure of the jet.

For a population of isotropically oriented jets, whose viewing angle probability distribution is P(θ_v) = sin θ_v, we have

$$\begin{array}{c}\hfill \begin{array}{cc}& P({\theta }_{\mathrm{v}},L,E_{\mathrm{p}}|{\mathit{\lambda }}_{\mathrm{pop}})=P(L,E_{\mathrm{p}}|{\theta }_{\mathrm{v}},{\mathit{\lambda }}_{\mathrm{pop}})P({\theta }_{\mathrm{v}})=\hfill \\ \hfill & \delta (L-L_{\mathrm{c}}\ell ({\theta }_{\mathrm{v}},{\mathit{\lambda }}_{\mathrm{pop}}))\delta (E_{\mathrm{p}}-E_{\mathrm{p},\mathrm{c}}\eta ({\theta }_{\mathrm{v}},{\mathit{\lambda }}_{\mathrm{pop}}))P(L_{\mathrm{c}},E_{\mathrm{p},\mathrm{c}}|{\mathit{\lambda }}_{\mathrm{pop}})sin{\theta }_{\mathrm{v}}.\hfill \end{array}\end{array}$$(6)

The induced joint luminosity and peak photon energy distribution, given the isotropic viewing angles, is then given by

$$\begin{array}{c}\hfill \begin{array}{cc}& P(L,E_{\mathrm{p}}|{\mathit{\lambda }}_{\mathrm{pop}})={\displaystyle {\int }_0^{\pi /2}}P(L,E_{\mathrm{p}}|{\theta }_{\mathrm{v}},{\mathit{\lambda }}_{\mathrm{pop}})sin{\theta }_{\mathrm{v}}\mathrm{d}{\theta }_{\mathrm{v}}\hfill \\ \hfill & ={\displaystyle {\int }_0^{\pi /2}}\frac{exp\{-{[yln(\ell ({\theta }_{\mathrm{v}})L_{\mathrm{c}}^{\star }/L)-ln(E_{\mathrm{p},\mathrm{c}}^{\star }\eta ({\theta }_{\mathrm{v}})/E_{\mathrm{p}})]}^2/2{\sigma }_{\mathrm{c}}^2\}}{E_{\mathrm{p}}\sqrt{2\pi {\sigma }_{\mathrm{c}}^2}}\hfill \\ \hfill & \times \frac{A}{\stackrel{\sim }{\mathrm{\Gamma }}(1-1/A)\ell ({\theta }_{\mathrm{v}})L_{\mathrm{c}}^{\star }}exp[-{\left(\frac{\ell ({\theta }_{\mathrm{v}})L_{\mathrm{c}}^{\star }}{L}\right)}^A]{\left(\frac{L}{\ell ({\theta }_{\mathrm{v}})L_{\mathrm{c}}^{\star }}\right)}^{-A}sin{\theta }_{\mathrm{v}}\mathrm{d}{\theta }_{\mathrm{v}}.\hfill \end{array}\end{array}$$(7)

In general, this must be solved numerically, but in the next section we analyse two cases where the intrinsic dispersion is negligible (i.e. σ_c → 0 and A → ∞) and ℓ and η take simple forms, so that an analytical integration is possible: this will help in demonstrating the main features of the P(L, E_p) distribution induced by such a quasi-universal structure scenario.

As a consequence of the assumption of redshift independence of the jet structure parameters, the probability distribution of the population source parameters is P_pop(L, E_p, z | λ_pop) = P(L, E_p | λ_pop)P(z | λ_pop), where the redshift probability distribution can be expressed as

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

Here dV/dz is the differential comoving volume and $\dot{\rho }(z,{\mathit{\lambda }}_{\mathrm{pop}})$ is the SGRB rate density at redshift z. We parametrise the latter as a smoothly broken power law, namely

$$\begin{array}{c}\hfill \dot{\rho }(z,{\mathit{\lambda }}_{\mathrm{pop}})=R_0\frac{{(1+z)}^a}{1+{\left(\frac{1+z}{1+z_p}\right)}^{b+a}},\end{array}$$(9)

where a, b, and z_p are free parameters and R₀ is the local rate density of SGRBs with any viewing angle. This parametrisation reflects the functional form of the fitting function to the cosmic star formation rate (CSFR) from Madau & Dickinson (2014). Allowing for variations in its parameters a, b, and z_p covers a variety of possible evolutionary histories that rise, peak and decay with increasing (1 + z), in a way that is agnostic of the progenitor nature, but can still accommodate, for example, typical rate evolutions obtained by convolving the cosmic star formation history with a merger delay time distribution under the assumption of a compact binary merger progenitor (as done in e.g. D’Avanzo et al. 2014; Wanderman & Piran 2015).

2.3. Apparent jet structure models and the implied luminosity-peak energy distribution

2.3.1. Luminosity function

A simple and widely adopted parametric form for the jet structure is a Gaussian one, namely

$$\begin{array}{c}\hfill \begin{array}{cc}& \ell =exp(-\frac{1}{2}{\left(\frac{{\theta }_{\mathrm{v}}}{{\theta }_{\mathrm{c}}}\right)}^2),\hfill \\ \hfill & \eta =exp(-\frac{1}{2}{\left(\frac{{\theta }_{\mathrm{v}}}{{\theta }_{\mathrm{c},E_{\mathrm{p}}}}\right)}^2).\hfill \end{array}\end{array}$$(10)

In the absence of a dispersion in the core quantities, which formally corresponds to the limit σ_c → 0 and A → ∞, and in the case where ℓ and η are monotonic, Eq. (7) reduces to a change of variables from θ_v to either L or E_p applied to the viewing angle probability P(θ_v) = sin θ_v, that is (Pescalli et al. 2015; Salafia & Ghirlanda 2022),

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill P(L,E_{\mathrm{p}})& ={[{\left(L_{\mathrm{c}}\frac{\partial \ell }{\partial {\theta }_{\mathrm{v}}}\right)}^{-1}sin{\theta }_{\mathrm{v}}\delta (E_{\mathrm{p}}-\eta ({\theta }_{\mathrm{v}})E_{\mathrm{p},\mathrm{c}})]}_{{\theta }_{\mathrm{v}}={\ell }^{-1}(L/L_{\mathrm{c}})}\hfill \\ \hfill & ={[{\left(E_{\mathrm{p},\mathrm{c}}\frac{\partial \eta }{\partial {\theta }_{\mathrm{v}}}\right)}^{-1}sin{\theta }_{\mathrm{v}}\delta (L-\ell ({\theta }_{\mathrm{v}})L_{\mathrm{c}})]}_{{\theta }_{\mathrm{v}}={\eta }^{-1}(E_{\mathrm{p}}/E_{\mathrm{p},\mathrm{c}})},\hfill \end{array}\end{array}$$(11)

where ℓ⁻¹ is the inverse function of ℓ and η⁻¹ is the inverse function of η. In what follows, we show results using the first of the above equalities, which highlights the dependence on L, but the results using the second equality are entirely analogous and can be obtained by exchanging L ↔ E_p, θ_c ↔ θ_c, Ep and L_c ↔ E_p, c. In the Gaussian apparent structure case this yields, for L < L_c and E_p < E_p, c,

$$\begin{array}{cc}\hfill P(L,E_{\mathrm{p}})& =\frac{{\theta }_{\mathrm{c}}^2}{L}\frac{sin\left({\theta }_{\mathrm{c}}\sqrt{2ln(L_{\mathrm{c}}/L)}\right)}{{\theta }_{\mathrm{c}}\sqrt{2ln(L_{\mathrm{c}}/L)}}\delta (E_{\mathrm{p}}-E_{\mathrm{p},\mathrm{c}}{\left(\frac{L}{L_{\mathrm{c}}}\right)}^{{\theta }_{\mathrm{c}}^2/{\theta }_{\mathrm{c},E_{\mathrm{p}}}^2})\hfill \\ \hfill & \sim \frac{{\theta }_{\mathrm{c}}^2}{L}\delta (E_{\mathrm{p}}-E_{\mathrm{p},\mathrm{c}}{\left(\frac{L}{L_{\mathrm{c}}}\right)}^{{\theta }_{\mathrm{c}}^2/{\theta }_{\mathrm{c},E_{\mathrm{p}}}^2}),\hfill \end{array}$$(12)

where the last approximate equality is valid for θ_v ≪ π/2, which corresponds to $L\gg L_{\text{c}}\exp (-1/2{\theta }_{\text{c}}^2)$ (or $E_{\text{p}}\gg E_{\text{p},\text{c}}\exp (-1/2{\theta }_{\text{c},E_{\text{p}}}^2)$). For typical values θ_c  ≪  1 (or θ_c, Ep  ≪  1), the exponential factor is tiny and hence the approximation applies to essentially all relevant luminosities and E_p’s. The luminosity function induced by a Gaussian universal apparent structure is therefore uniform in log(L), and the same applies to the E_p distribution.

The effect of a non-zero dispersion in the core quantities is equivalent to that of convolving the zero-dispersion P(L, E_p) probability density distribution with the probability density distribution of the core luminosity and peak SED photon energy P(L_c, E_p, c). In that sense, the probability density distribution of the core quantities acts essentially as a smoothing kernel: in the Gaussian apparent structure case, it introduces a smooth transition to a power law fall-off above $L_{\text{c}}^{\star }$ and $E_{\text{p},\text{c}}^{\star }$. In more physical terms, and for typical parameters relevant to the quasi-universal jet scenario, the high-end of the luminosity function is shaped by the distribution of core luminosities, while below $L_{\text{c}}^{\star }$ the luminosity function is set by the jet (apparent) structure. The left-hand panel in Fig. 1 shows the luminosity function ϕ(L) = dP/dln(L) = L∫P(L, E_p) dE_p for an example Gaussian apparent structure case, demonstrating the effect of introducing a core luminosity dispersion with three different values of the slope A.

Fig. 1. Example luminosity functions ϕ(L) = dP/dln(L) induced by a Gaussian apparent structure model (Eqs. (10), left-hand panel) and a power law apparent structure model (Eqs. (13)) with a luminosity profile slope α = 3 (right-hand panel), both with a core half-opening angle θc = 0.1 rad and $L_{\text{c}}^{\star }$ = 1053 erg s−1. In each panel, the luminosity function corresponding to a population with zero core dispersion (A → ∞ in Eq. (2)) is shown in grey, and three cases with different core luminosity distributions (slopes reported in the legend) are shown with coloured lines. The corresponding Lc distributions (arbitrarily scaled) are shown with dashed lines. In the right-hand panel, the dashed black line shows the trend ϕ(L)∝L−2/αL expected at intermediate luminosities.

We can get some further insight by adopting a power law apparent structure model with a ‘uniform core’ within θ_v ≤ θ_c, namely

$$\begin{array}{c}\hfill \begin{array}{cc}& \ell =min[1,{\left(\frac{{\theta }_{\mathrm{v}}}{{\theta }_{\mathrm{c}}}\right)}^{-{\alpha }_L}]\hfill \\ \hfill & \eta =min[1,{\left(\frac{{\theta }_{\mathrm{v}}}{{\theta }_{\mathrm{c}}}\right)}^{-{\alpha }_{E_{\mathrm{p}}}}],\hfill \end{array}\end{array}$$(13)

again with no dispersion. The change of variables approach can be applied to θ_v > θ_c, where the structure is monotonic; for θ_v ≤ θ_c it is sufficient to note that all observers see L = L_c and E_p = E_p, c, and the probability of having a viewing angle in this range is 1 − cos θ_c. Hence, the L − E_p distribution is

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill P(L,E_{\mathrm{p}})& =\frac{{\theta }_{\mathrm{c}}{L}_{\mathrm{c}}^{1/{\alpha }_L}}{{\alpha }_L}\frac{sin\left[{\theta }_{\mathrm{c}}{\left(\frac{L_{\mathrm{c}}}{L}\right)}^{1/{\alpha }_L}\right]}{L^{1+1/{\alpha }_L}}\delta (E_{\mathrm{p}}-E_{\mathrm{p},\mathrm{c}}{\left(\frac{L}{L_{\mathrm{c}}}\right)}^{{\alpha }_{E_{\mathrm{p}}}/{\alpha }_L})\hfill \\ \hfill & \sim \frac{{\theta }_{\mathrm{c}}^2{L}_{\mathrm{c}}^{2/{\alpha }_L}}{{\alpha }_L}{L}^{-1-2/{\alpha }_L}\delta (E_{\mathrm{p}}-E_{\mathrm{p},\mathrm{c}}{\left(\frac{L}{L_{\mathrm{c}}}\right)}^{{\alpha }_{E_{\mathrm{p}}}/{\alpha }_L})\hfill \end{array}\end{array}$$(14)

for $L_{\text{min}}^{\star }\le L<L_{\text{c}}$ and $E_{\text{p},\text{min}}^{\star }\le E_{\text{p}}<E_{\text{p},\text{c}}$, where $L_{\text{min}}^{\star }=L_{\text{c}}{(2{\theta }_{\text{c}}/\pi )}^{{\alpha }_L}$ and $E_{\text{p},\text{min}}^{\star }=E_{\text{p},\text{c}}{(2{\theta }_{\text{c}}/\pi )}^{{\alpha }_{E_{\text{p}}}}$, and the last approximate equality is valid for $L\gg L_{\text{min}}^{\star }$ and $E_{\text{p}}\gg E_{\text{p},\text{min}}^{\star }$. Again, the alternate form, which highlights the dependence on E_p, can be obtained by the substitutions L ↔ E_p, α_L ↔ α_Ep and L_c ↔ E_p, c. For L ≥ L_c, E_p ≥ E_p, c, we have

$$\begin{array}{c}\hfill P(L,E_{\mathrm{p}})=(1-cos{\theta }_{\mathrm{c}})\delta (L-L_{\mathrm{c}})\delta (E_{\mathrm{p}}-E_{\mathrm{p},\mathrm{c}}).\end{array}$$(15)

While in this limiting zero-dispersion case the jet core clearly extends over a zero-measure region of the (L, E_p) plane, a non-zero dispersion spreads this over a more physically sound, finite region of the plane. Hence, a power law apparent structure induces a power law luminosity function with a slope L^−2/α_L in the logarithm of L (Pescalli et al. 2015). Similar conclusions apply to the induced E_p distribution. The slope parameter α_L, together with the core half-opening angle θ_c, control the extent of the luminosity function as a consequence of the limited physical viewing angle range 0 ≤ θ_v ≤ π/2. The right-hand panel in Fig. 1 shows the luminosity function in an example power law case with α_L = 3.

The above derivation also shows that in a quasi-universal structured jet scenario ϕ(L) can never attain a positive slope, except at the low-luminosity end, or in a narrow luminosity range close to the core for small dispersions (unless α_L < 0, which however does not seem a likely physical possibility). This is a feature of quasi-universal jet models.

2.3.2. L − E_p correlation induced by the jet structure

Since both ℓ and η are functions of the viewing angle, the quasi-universal structured jet model inherently implies a Yonetoku correlation (Yonetoku et al. 2004) between L and E_p within the observed population (e.g. Salafia et al. 2015), in addition to any intrinsic correlation that may hold between the core quantities L_c and E_p, c. The shape of the induced correlation can be obtained by eliminating θ_v in the average apparent structure functions, to obtain η(ℓ), and is already apparent in the Dirac-delta functions in Eqs. (12) and (14). In the Gaussian average apparent structure case, it yields $(E_{\text{p}}/E_{\text{p},\text{c}})={(L/L_{\text{c}})}^{{\theta }_{\text{c}}^2/{\theta }_{\text{c},E_{\text{p}}}^2}$, so that the slope of the correlation is set by the ratio of the scale parameters θ_c and θ_c, Ep over which L and E_p decay with the viewing angle. Similarly, in the power law case the relation is (E_p/E_p, c) = (L/L_c)^α_Ep/α_L. In general, the induced correlation is a power law (or a collection of power law branches) whenever the functional forms of ℓ and η are the same, and its slope is set by the ratio of the decay rates of ℓ and η.

2.3.3. Average apparent structure model adopted in this work

In this study, in order to endow the average apparent structure model with a high degree of flexibility (in the absence of strong theoretical constraints on the expected shape), we adopted a double smoothly broken power law model with a nearly constant core within θ_v < θ_c and a break at a wider angle θ_w, that is,

$$\begin{array}{c}\hfill \begin{array}{cc}& \ell ({\theta }_{\mathrm{v}})={[1+{\left(\frac{{\theta }_{\mathrm{v}}}{{\theta }_{\mathrm{c}}}\right)}^s]}^{-{\alpha }_L/s}{[1+{\left(\frac{{\theta }_{\mathrm{v}}}{{\theta }_{\mathrm{w}}}\right)}^s]}^{-({\beta }_L-{\alpha }_L)/s}\hfill \\ \hfill & \eta ({\theta }_{\mathrm{v}})={[1+{\left(\frac{{\theta }_{\mathrm{v}}}{{\theta }_{\mathrm{c}}}\right)}^s]}^{-{\alpha }_{E_{\mathrm{p}}}/s}{[1+{\left(\frac{{\theta }_{\mathrm{v}}}{{\theta }_{\mathrm{w}}}\right)}^s]}^{-({\beta }_{E_{\mathrm{p}}}-{\alpha }_{E_{\mathrm{p}}})/s},\hfill \end{array}\end{array}$$(16)

where the ‘smoothness’ parameter is set to s = 4, which makes the transitions between the power law branches relatively sharp, to compensate the fact that the intrinsic dispersion of core quantities tends to smooth out the induced breaks in the luminosity function. The implied L − E_p correlation has two branches, with slopes α_Ep/α_L and β_Ep/β_L, the break being around $L\sim L_{\text{w}}^{\star }=L_{\text{c}}^{\star }{({\theta }_{\text{w}}/{\theta }_{\text{c}})}^{-{\alpha }_L}$ and $E_{\text{p}}\sim E_{\text{p},\text{w}}^{\star }=E_{\text{p},\text{c}}^{\star }{({\theta }_{\text{w}}/{\theta }_{\text{c}})}^{-{\alpha }_{E_{\text{p}}}}$. Figure 2 demonstrates the features of such a model in detail, including the induced L and E_p probability distributions, based on an example choice of parameters.

Fig. 2. Example double smoothly broken power law jet apparent structure model and induced L and Ep probability distributions. Panels a.1 and a.2 show the apparent jet structure functions $L_{\text{c}}^{\star }\mathcal{l}({\theta }_{\text{v}})$ and $E_{\text{p},\text{c}}^{\star }\eta ({\theta }_{\text{v}})$ (solid red and green line, respectively) adopting the parametrisation defined in Eq. (16) and using the parameter values reported in the middle panel. Panel a.0 shows the induced Ep(L) relation due to the structure (purple line). In panels a.1 and a.2, the dotted black lines show the transition angles θc and θw. The dashed lines demonstrate the analytical trends discussed in the text. Panel b.1 shows the induced luminosity function ϕ(L) = dP/dlnL = L∫P(L, Ep) dEp (red line), with dashed lines demonstrating approximate analytical trends (significant departures are due to the core dispersion), and with the core luminosity distribution shown with a dashed dark grey line, centred at Lc and with its integral normalised to 1 − cos θc. Panel b.2 is the analogue of b.1, but showing the dP/dlnEp probability distribution (green line).

With this choice of average apparent jet structure functions, the model features a total of 14 free parameters. We list these parameters in Table 1, along with brief definitions and with information on the priors adopted on each of them in the analysis, discussed later in the text.

2.4. Sample definition and selection effect modelling

In order to compare a population model with an observed sample, the selection effects that shape the latter must be taken into account. Here we describe our sample choice and the procedure that we employed to model the underlying selection effects. We adopted a description of the SGRB photon spectrum as a cut-off power law (Ghirlanda et al. 2004), dṄ/dE(E,E_p,obs,α) ∝ E^αexp[−(2+α)E/E_p,obs], where Ṅ represents the rate of photons hitting the detector, E is the photon energy, and E_p, obs = E_p/(1 + z). We set the low-energy photon index to α = −0.4, which is the median of the values reported in the Fermi/GBM online catalogue for SGRBs. The peak photon flux in the E₀ − E₁ keV band is then defined as

$$\begin{array}{c}\hfill p_{[E_0,E_1]}(L,E_{\mathrm{p}},z)=\frac{1}{4\pi d_{\mathrm{L}}^2(z)}\frac{L}{{\mathcal{E}}_{[E_0,E_1]}(E_{\mathrm{p}},z)},\end{array}$$(17)

where d_L is the luminosity distance and

$$\begin{array}{c}\hfill {\mathcal{E}}_{[E_0,E_1]}(E_{\mathrm{p}},z)=\frac{{\int }_{0.1\mathrm{keV}/(1+z)}^{{10}^7\mathrm{keV}/(1+z)}\mathrm{d}\dot{N}/\mathrm{d}E(E,E_{\mathrm{p},\mathrm{obs}},\alpha )E\mathrm{d}E}{{\int }_{E_0\mathrm{keV}}^{E_1\mathrm{keV}}\mathrm{d}\dot{N}/\mathrm{d}E(E,E_{\mathrm{p},\mathrm{obs}},\alpha )\mathrm{d}E}.\end{array}$$(18)

We stress here again that we extended the customary 1 − 10⁴ keV pseudo-bolometric rest-frame band to the wider 0.1 − 10⁷ keV to include possible cases with very high E_p. In practice, we pre-computed³ ℰ_{[50 − 300]} for Fermi/GBM and ℰ_{[15 − 150]} for Swift/BAT over a uniformly spaced two-dimensional grid in (log E_p, log z) and then used two-dimensional linear interpolation to recover it and obtain the photon flux from Eq. (17) (or equivalently to obtain L from p and E_p, obs, by inverting the equation).

We considered three reference SGRB samples: (i) SGRBs detected by Fermi/GBM, with available spectral information in the public catalogue; (ii) SGRBs detected by both Swift/BAT and Fermi/GBM, with a number of additional cuts to reach a high completeness in redshift; and (iii) SGRBs detected by Fermi/GBM with a GW counterpart, which currently includes only GRB 170817A/GW170817. In what follows, we describe in detail the selection cuts of each sample and our approach to the modelling of selection effects.

2.4.1. Observer-frame sample: Fermi/GBM SGRBs with spectral information

Our ‘observer-frame’ sample includes GRBs detected by the GBM on board Fermi from the start of the mission up to mid July 2018, after which no spectral information is present in the online catalogue (von Kienlin et al. 2020) at the time of writing. Among these, we selected 367 events that are nominally short, that is, their duration T₉₀ < 2 s, as our initial raw sample. As a first approach, we considered working with the sub-sample of bursts whose peak photon flux is larger than a ‘completeness’ threshold p_{[50 − 300]} > p_lim, GBM, where the GBM flux in the 50–300 keV band is measured with 64 ms binning to best approximate the actual peak photon flux. By visual inspection we found that above p_lim, GBM = 3.5 cm⁻² s⁻¹ the observed N(> p_{[50 − 300]}) distribution looks like a single power law (see Fig. 3), and hence we selected this as our completeness threshold. This is a practical approach historically employed to construct flux-complete samples. The corresponding detection probability can be modelled simply as

$$\begin{array}{c}\hfill P_{\det ,\mathrm{GBM}}(L,E_{\mathrm{p}},z)=\mathrm{\Theta }(p_{[50-300]}(L,E_{\mathrm{p}},z)-p_{\lim ,\mathrm{GBM}}),\end{array}$$(19)

Fig. 3. Inverse cumulative distribution of Fermi/GBM SGRB peak photon fluxes. The solid red line shows the inverse cumulative number of SGRBs detected by GBM with spectral information available in the catalogue, as a function of the peak photon flux measured on a 64 ms timescale in the 50–300 keV band. The pink band shows the one-sigma-equivalent Poisson error. The dashed black line shows a power law $p_{[50-300]}^{-3/2}$ for reference. The vertical dotted line shows the chosen value of the flux cut for the flux-limited sample, plim, GBM = 3.5 cm−2 s−1.

where Θ is the Heaviside step function, that is, Θ(x) = 1 if x > 0 and Θ(x) = 0 otherwise.

This selection criterion significantly reduces the sample size: out of a total of 367 SGRBs in the raw sample, only 212 have p_{[50 − 300]} > 3.5 cm⁻² s⁻¹. Most importantly, the discarded events possibly probe the luminosity function down to lower luminosities, which is where most of the useful information on the jet structure resides in a quasi-universal jet scenario. Last, but not least, the only SGRB with reliable viewing angle information, that is, GRB 170817A, is not included in this flux-limited sample. In order to access the flux-incomplete part of the Fermi/GBM SGRB sample, we carefully constructed a detection probability P_det(L, E_p, z) by simulating the response of the GBM NaI detectors to SGRBs with a broad range of characteristics, as described in detail in Appendix C. In this study, we compare the results obtained by using either of the strategies in modelling the selection effects: hereafter, we refer to the reduced sample of Fermi/GBM SGRBs with p_{[50 − 300]} > 3.5 cm⁻² s⁻¹ as the ‘flux-limited sample’, and to the analysis adopting the associated simplified selection effect model as the ‘flux-limited sample analysis’. Conversely, the analysis performed adopting the simulated Fermi/GBM detection efficiency is referred to as the ‘full sample analysis’.

As a further countermeasure against possible biases, we applied an additional quality cut to both the above samples: we removed events with best-fit E_p, obs > 10 MeV or E_p, obs < 50 keV, which fall outside the spectral range where the effective area of the GBM detectors is optimal: this removes two events whose uncertainty on E_p, obs is very large, reducing the flux-limited sample to 210 events. In the full sample analysis, we also removed another ten events with a low best-fit peak 64 ms photon flux p_{[50 − 300]} < 1 cm⁻² s⁻¹, all of which have very large errors on both p_{[50 − 300]} and E_p, obs. This reduces the ‘full’ sample to 355 events. To reflect these further quality cuts, we updated our model detection probabilities by multiplying them by Θ(p_{[50 − 300]} − 1 cm⁻² s⁻¹)Θ(E_p, obs − 50 keV)Θ(10 MeV − E_p, obs).

2.4.2. Rest-frame sample: Flux-complete sample of Swift/BAT SGRBs also observed by Fermi/GBM

For some of the events in our Fermi/GBM observer-frame sample, redshift information is available and can be used to constrain the population parameters better than can be done using the observer-frame information (p_{[50 − 300]}, E_p, obs) only. In order to avoid biases, on the other hand, any additional selection effects at play in the sub-sample with known redshift must be accounted for in the inference, that is, in the P_det model. The redshift determination is a complex process that involves multiple facilities and depends not only on the prompt emission properties, but also on those of the afterglow. Therefore, modelling the associated selection effects is prohibitive. On the other hand, Salvaterra et al. (2012) and D’Avanzo et al. (2014) showed that it is possible to construct a sample with a selection that is easier to model, but that leads to a high redshift completeness. The selection involves two cuts that do not bias the redshift distribution, namely (i) a cut on the foreground interstellar dust extinction A_V and (ii) a cut on the Swift X-Ray Telescope (XRT) slew time; plus a cut on the BAT 64 ms peak flux in the 15–150 keV band, p_{[15 − 150]} > 3.5 cm⁻² s⁻¹, to ensure flux completeness. The original SGRB sample constructed in this way, known as the S-BAT4 (D’Avanzo et al. 2014), included 16 events, 11 of which had a measured redshift. Thanks to a considerable effort spent by the community, and in particular by Fong et al. (2022) and Nugent et al. (2022), in identifying SGRB host galaxies and measuring their redshifts, it has been recently possible to construct an extended ‘S-BAT4ext’ sample (Ferro et al. 2023; D’Avanzo et al., in prep.) with a more than doubled size and an increased redshift completeness. For this work, we adopt the S-BAT4ext sub-sample of 18 events that have been jointly detected by Fermi/GBM with T₉₀ < 2 s (for consistency with our observer-frame sample). Of these, 16 have either a spectroscopic (12 events) or photometric (4 events) redshift measurement, as listed in Table 2. For these SGRBs we performed an independent analysis of their peak spectra, based on publicly available Fermi/GBM data, as described in Appendix B, obtaining posterior samples of their bolometric flux F and observed peak photon energy E_p, obs, adopting a prior π(F, E_p, obs)∝F⁻¹. For the 12 events with a spectroscopic redshift, these were converted into samples of P(L, E_p, z | d_i) by simply fixing the redshift at the best fit value. For the four events with photometric redshift, we obtained samples of P(L, E_p, z | d_i) by computing $L_{l,m}=4\pi d_{\text{L},l}^2{F}_m$ and E_p, l, m = (1 + z_l)E_{p, obs, m}, where ${\left\{z_l\right\}}_{l=1}^{N_{\text{s}}}$ are N_s photometric redshift posterior samples from the Broad-band Repository for Investigating Gamma-ray burst Host galaxies Traits (BRIGHT) catalogue (Nugent et al. 2022), ${\left\{d_{\text{L},l}\right\}}_{l=1}^{N_{\text{s}}}$ are the corresponding luminosity distances under our assumed cosmology and ${\{(F_m,E_{\text{p},\text{obs},m})\}}_{m=1}^{M_{\text{s}}}$ are M_s posterior samples from the spectral analysis. In both cases, the effective prior on the source parameters is π(L, E_p, z)∝(1 + z)⁻¹L⁻¹. The two events with unknown redshift were not included in the analysis.

The careful selection adopted to construct this sample allowed us to model the underlying selection effects by multiplying the GBM detection efficiency, P_det, GBM, by the simple BAT detection efficiency:

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

where p_lim, BAT = 3.5 cm⁻² s⁻¹ is numerically identical to p_lim, GBM by pure chance.

2.4.3. Viewing angle sample: GRB 170817A/GW170817

The third and last sample we considered is that of Fermi/GBM SGRBs that have a GW counterpart produced by the inspiral of a BNS merger, from which a measurement of θ_v can be obtained under the assumption that the jet is launched along the direction of the total angular momentum. At present, the sample clearly consists of the single event GRB 170817A with its counterpart GW170817. We use d_G17 to indicate the available data regarding the prompt emission and the GW signal of the event, and with d_HOST the available data of the host galaxy NGC4993 (Coulter et al. 2017; Hjorth et al. 2017; Cantiello et al. 2018). We took the host galaxy spectroscopic redshift z_H = 0.009783 (Coulter et al. 2017; Hjorth et al. 2017) as the redshift of the source, neglecting the small uncertainty on the actual cosmological redshift (i.e. corrected for the galaxy proper motion), and hence P(z | d_HOST) = δ(z − z_H). We obtained the posterior P(L, E_p | d_G17, z_H), whose 50% and 90% containment contours are shown in the left-hand panel of Fig. 4, from our own analysis of the peak spectrum (Appendix B). For what concerns the viewing angle, in order to break the distance-inclination degeneracy inherent in the BNS inspiral GW analysis, we proceeded similarly as in Mandel (2018, but keeping the cosmological parameters fixed, differently from them), as follows: we approximated the host galaxy luminosity distance r uncertainty as

$$\begin{array}{c}\hfill P(r|d_{\mathrm{HOST}})=\frac{exp[-\frac{1}{2}{\left(\frac{r-{\mu }_r}{{\sigma }_r}\right)}^2]}{\sqrt{2\pi {\sigma }_r^2}},\end{array}$$(21)

Fig. 4. L − Ep contours for our sub-sample of Fermi/GBM SGRBs detected also by Swift/BAT with a 64 ms peak flux p[15 − 150] > 3.5 cm−2 s−1 (see the main text for the full sample selection criteria). The main panel shows the contours that contain 50% (thick lines) and 90% (thin lines) of the posterior probability on (L, Ep) for bursts with a spectroscopically measured redshift, plus GRB 170817A (which does not belong to the GBM+BAT sample); the right-hand inset shows the corresponding contours for four events with a photometric redshift measurement (from Fong et al. 2022); the left-hand inset shows the contours for the remaining two SGRBs with an unknown redshift, constructed assuming a uniform prior on z in the range (10−4, 4). These two events were not included in the analysis.

with μ_r = 40.7 Mpc and σ_r = 2.36 Mpc (mean and square sum of statistical and systematic errors from the measurements performed by Cantiello et al. 2018). We then computed the posterior on the viewing angle as

$$\begin{array}{c}\hfill P({\theta }_{\mathrm{v}}|d_{\mathrm{G}17},d_{\mathrm{HOST}})={\displaystyle {\int }_0^{\infty }}P({\theta }_{\mathrm{v}},r|d_{\mathrm{G}17})P(r|d_{\mathrm{HOST}})\mathrm{d}r,\end{array}$$(22)

where P(θ_v, r | d_G17) is the joint posterior on θ_v and r from the GW analysis. In practice, we obtained samples of the posterior probability density in Eq. (22) by re-sampling the publicly available posterior samples from the low-spin prior analysis of GW170817 performed by Abbott et al. (2019) with a weight equal to the right-hand side of Eq. (21) evaluated at the luminosity distance of each sample. The viewing angle posterior probability density obtained from a kernel density estimate on the resulting samples is shown in Fig. 5.

Fig. 5. GRB 170817A viewing angle posterior probability distribution, assuming the jet to be aligned with the GW170817 binary total angular momentum. The blue line shows the posterior probability constructed using the posterior samples from the low-spin-prior GW analysis (Abbott et al. 2019), while the red line shows the result of conditioning on the host galaxy distance (Cantiello et al. 2018), as explained in the main text.

We modelled the selection effects acting on this sample as the product of P_det, GBM(L, E_p, z) times a GW detection efficiency P_det, GW(θ_v, z). Since the time-volume surveyed by aLIGO and Advanced Virgo so far is dominated by that of the third observing run, O3 (LIGO Scientific Collaboration 2021), we constructed the detection efficiency assuming, for the sake of simplicity, the GW network sensitivity of O3, neglecting periods with a lower sensitivity. To do so, we retrieved the dataset made publicly available by the LIGO Scientific Collaboration et al. (2023) that contains information on the response of online GW search pipelines to a large number of simulated signals injected into O3 data. We re-weighted the BNS merger injections in that dataset to reflect a population with a primary mass distribution $p(m_1)\propto m_1^{-2}$ between m_1, min = 1.2 M_⊙ and m_1, max = 2.1 M_⊙ (similar to the preferred distribution from the GWTC-3 population analysis, Abbott et al. 2023) and a flat secondary mass distribution p(m₂ | m₁)∝Θ(m₁ − m₂)Θ(m₂ − m_1, min). Assuming the inclination ι and the jet viewing angle θ_v to be related by

$$\begin{array}{c}\hfill {\theta }_{\mathrm{v}}(\iota )=\{\begin{array}{cc}\iota \hfill & \hfill 0\le \iota <\pi /2\\ \pi -\iota \hfill & \hfill \pi /2\le \iota <\pi \end{array},\end{array}$$(23)

we binned the injected signals into a number of two-dimensional bins in (θ_v, z) space centred at (θ_v, i, z_j). Calling w_k and ρ_k the weight and network signal-to-noise ratio (S/N) associated with the k-th injected signal, we estimated P_det, GW(θ_v, z) at the centre of each bin as

$$\begin{array}{c}\hfill P_{\det ,\mathrm{GW}}({\theta }_{\mathrm{v},i},z_j)\sim \frac{{\sum }_{k\in {\mathcal{I}}_{i,j}}{w}_k\mathrm{\Theta }({\rho }_k\ge 12)}{{\sum }_{k\in {\mathcal{I}}_{i,j}}{w}_k},\end{array}$$(24)

where ℐ_i, j represents the set of indices k of injections whose viewing angle and redshift fall into the bin centred at (θ_v, i, z_j). The GW detection efficiency on the rest of the (θ_v, z) space was obtained by two-dimensional linear interpolation of the resulting values. The relatively high cut ρ ≥ 12 includes GW170817 (Abbott et al. 2017a) and ensures that the detection can be represented by a simple cut in S/N, in analogy with the flux completeness cuts discussed previously.

2.5. Inference on the population properties

Within a Bayesian hierarchical approach, the posterior probability on the population parameters λ_pop can be written as (Eqs. (7) and (8) in Mandel et al. 2019, hereafter M19)

$$\begin{array}{c}\hfill P({\mathit{\lambda }}_{\mathrm{pop}}|\mathit{d})=\frac{\pi ({\mathit{\lambda }}_{\mathrm{pop}})}{P(\mathit{d})}{\displaystyle \prod _{i=1}^N}\frac{\int P(d_i|{\mathit{\lambda }}_{\mathrm{src}})P_{\mathrm{pop}}({\mathit{\lambda }}_{\mathrm{src}}|{\mathit{\lambda }}_{\mathrm{pop}})\mathrm{d}{\mathit{\lambda }}_{\mathrm{src}}}{\int P_{\mathrm{pop}}({\mathit{\lambda }}_{\mathrm{src}}|{\mathit{\lambda }}_{\mathrm{pop}})P_{\det }({\mathit{\lambda }}_{\mathrm{src}})\mathrm{d}{\mathit{\lambda }}_{\mathrm{src}}},\end{array}$$(25)

where the index i runs over the N events in the sample, d_i (the i-th element of d) represents the corresponding data in the detectors (i.e. Fermi/GBM, plus either Swift/BAT or the aLIGO and Advanced Virgo detectors in our case), π(λ_pop) is the prior on the population parameters, and $P(\mathit{d})=\prod _{i=1}^{N}P(d_i)$ is a normalisation factor. We indicate with 𝒩_i = 𝒩_i(d_i, λ_pop) the numerator of the fraction after the product symbol in the above equation, that is,

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

and with 𝒟 = 𝒟(λ_pop) the denominator, namely

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

This equation contains the detection efficiency P_det(λ_src), which must be chosen consistently with the selection effects acting on each of the considered samples.

2.5.1. Numerator 𝒩_i for observer-frame sample events

For events in our observer-frame sample, which have unknown redshifts, we approximated the posterior probability on the source parameters P(L, E_p, z | d_i) = P(d_i | λ_src)π(L, E_p)π(z)/P(d_i) by neglecting the uncertainty on the measured peak photon flux p_i = p_{[50 − 300],i} and peak photon energy E_{p, obs, i} (which is justified by the large sample size and by the quality cuts discussed in the previous section), and by assuming the posterior to simply reflect the prior π(z) along the z axis. In other words, we assumed P(L, E_p, z | d_i)∝π(z) when keeping L and E_p fixed: this is equivalent to stating that no redshift information is available. This leads to

$$\begin{array}{c}\hfill P(d_i|{\mathit{\lambda }}_{\mathrm{src}})\sim \frac{\delta (L-4\pi d_{\mathrm{L}}^2\mathcal{E}{p}_i)\delta (E_{\mathrm{p}}-(1+z)E_{\mathrm{p},\mathrm{obs},i})}{\pi (L,E_{\mathrm{p}})}P(d_i).\end{array}$$(28)

The appropriate prior π(L, E_p) was obtained by applying a coordinate transform to a uniform prior in both p_i and E_p, obs for these events, that is,

$$\begin{array}{c}\hfill \pi (L,E_{\mathrm{p}})=\overline{)\begin{array}{cc}\frac{\partial p_i}{\partial L}& \frac{\partial p_i}{\partial E_{\mathrm{p}}}\\ \frac{\partial E_{\mathrm{p},\mathrm{obs},i}}{\partial L}& \frac{\partial E_{\mathrm{p},\mathrm{obs},i}}{\partial E_{\mathrm{p}}}\end{array}}\pi (p_i,E_{\mathrm{p},\mathrm{obs},i})\propto \frac{p_i}{(1+z){\widehat{L}}_i(z)},\end{array}$$(29)

where |X| represents the determinant of matrix X, and ${\widehat{L}}_i=4\pi d_{\mathrm{L}}^2(z){\mathcal{E}}_{[50-300]}((1+z)E_{\mathrm{p},\mathrm{obs},i},z)p_i$. The integral in the numerator of Eq. (25) is then

$$\begin{array}{c}\hfill {\mathcal{N}}_i\sim P(d_i){\displaystyle {\int }_0^{\infty }}(1+z)\frac{{\widehat{L}}_i(z)}{p_i}{P}_{\mathrm{pop}}({\widehat{L}}_i(z),(1+z)E_{\mathrm{p},\mathrm{obs},i},z|{\mathit{\lambda }}_{\mathrm{pop}})\mathrm{d}z,\end{array}$$(30)

and we note that the P(d_i) term eventually cancels out in Eq. (25) (as in M19).

2.5.2. Numerator 𝒩_i for rest-frame sample events

For events in the rest-frame sample, in order to carry out the integral over z, L and E_p in 𝒩_i, we employed the same Monte Carlo approximation as M19, that is,

$$\begin{array}{c}\hfill {\mathcal{N}}_i\sim \frac{P(d_i)}{N_{\mathrm{s}}}{\displaystyle \sum _{j=1}^{N_{\mathrm{s}}}}\frac{P_{\mathrm{pop}}(L_{i,j},E_{\mathrm{p},i,j},z_{i,j}|{\mathit{\lambda }}_{\mathrm{pop}})}{\pi (L_{i,j},E_{\mathrm{p},i,j})},\end{array}$$(31)

where ${\{L_{i,j},E_{\text{p},i,j},z_{i,j}\}}_{j=1}^{N_{\text{s}}}$ represent a total of N_s samples of the luminosity, peak photon energy, and redshift posterior of the GRB. In cases with a photometric redshift, the (L, E_p) posterior was obtained by combining the results from our analysis of the spectrum at the peak of the GRB with the redshift posterior samples from Nugent et al. (2022), obtained from the BRIGHT online catalogue⁴. In cases with a spectroscopic redshift, we simply kept z_i, j fixed at the best-fit redshift. The prior π(L, E_p) in Eq. (31) is proportional to (1 + z)⁻¹L⁻¹, as discussed in Sect. 2.4.2.

2.5.3. Numerator 𝒩_i for viewing angle sample events

For the viewing angle sample, that is, for GRB 170817A, it was necessary to proceed differently for the full sample and flux-limited sample analyses: the peak photon flux of the GRB, p_{[50 − 300]} = 2.02 cm⁻² s⁻¹, is below the completeness cut p_lim, GBM = 3.5 cm⁻² s⁻¹ adopted in the flux-limited sample analysis; hence, in this case, the simple treatment of GBM selection effects is not adequate. In the full sample analysis, this is not a problem since the GBM detection efficiency model from Appendix C accounts for the smooth decrease in the detection efficiency at low peak photon fluxes. In the full-sample case, the correct form of 𝒩_i can be obtained by reformulating our population model by including θ_v among the source parameters, ${\mathit{\lambda }}_{\text{src}}^{\star }=(L,E_{\text{p}},z,{\theta }_{\text{v}})$. The population probability then becomes $P_{\text{pop}}^{\star }(L,E_{\text{p}},z,{\theta }_{\text{v}}|{\mathit{\lambda }}_{\text{pop}})=P(L,E_{\text{p}}|{\theta }_{\text{v}},{\mathit{\lambda }}_{\text{pop}})sin{\theta }_{\text{v}}P(z|{\mathit{\lambda }}_{\text{pop}})$, and the integral over θ_v is deferred to 𝒩_i, namely

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathcal{N}}_i& =\int \int \int \int P(d_i|L,E_{\mathrm{p}},z,{\theta }_{\mathrm{v}})\hfill \\ \hfill & \times P(L,E_{\mathrm{p}}|{\theta }_{\mathrm{v}},{\mathit{\lambda }}_{\mathrm{pop}})sin{\theta }_{\mathrm{v}}P(z|{\mathit{\lambda }}_{\mathrm{pop}})\mathrm{d}L\mathrm{d}{E}_{\mathrm{p}}\mathrm{d}z\mathrm{d}{\theta }_{\mathrm{v}}.\hfill \end{array}\end{array}$$(32)

It is straightforward to verify that whenever d_i does not contain information on the viewing angle this leads to exactly the same definition of 𝒩_i as before.

For GRB 170817A, we have (see Sect. 2.4.3)

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill P(d_{\mathrm{G}17},d_{\mathrm{HOST}}|L,E_{\mathrm{p}},z,{\theta }_{\mathrm{v}})& =\frac{P({\theta }_{\mathrm{v}}|d_{\mathrm{G}17},d_{\mathrm{HOST}})}{\pi ({\theta }_{\mathrm{v}})}\hfill \\ \hfill & \times \frac{P(L,E_{\mathrm{p}}|d_{\mathrm{G}17},z)}{\pi (L,E_{\mathrm{p}},z)}\delta (z-z_{\mathrm{H}})P(d_{\mathrm{G}17},d_{\mathrm{HOST}}),\hfill \end{array}\end{array}$$

where the prior π(θ_v) = sin θ_v corresponds to that used in the GW analysis. This leads to

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \frac{{\mathcal{N}}_{\mathrm{G}17}}{P(d_{\mathrm{G}17},d_{\mathrm{HOST}})}& ={\displaystyle \int \int \int }P({\theta }_{\mathrm{v}}|d_{\mathrm{G}17},d_{\mathrm{HOST}})\frac{P(L,E_{\mathrm{p}},z_{\mathrm{H}}|d_{\mathrm{G}17})}{\pi (L,E_{\mathrm{p}},z_{\mathrm{H}})}\hfill \\ \hfill & \times P(L,E_{\mathrm{p}}|{\theta }_{\mathrm{v}},{\mathit{\lambda }}_{\mathrm{pop}})P(z_{\mathrm{H}}|{\mathit{\lambda }}_{\mathrm{pop}})\mathrm{d}L\mathrm{d}{E}_{\mathrm{p}}\mathrm{d}{\theta }_{\mathrm{v}}.\hfill \end{array}\end{array}$$

Approximating again the integral with a Monte Carlo sum, we obtained

$$\begin{array}{c}\hfill \frac{{\mathcal{N}}_{\mathrm{G}17}}{P(d_{\mathrm{G}17},d_{\mathrm{HOST}})}\sim \frac{P(z_{\mathrm{H}}|{\mathit{\lambda }}_{\mathrm{pop}})}{N_{\mathrm{g}}{N}_{\mathrm{s}}}{\displaystyle \sum _{j=1}^{N_{\mathrm{g}}}\sum _{k=1}^{N_{\mathrm{s}}}}\frac{P(L_k,E_{\mathrm{p},k}|{\theta }_{\mathrm{v},j},{\mathit{\lambda }}_{\mathrm{pop}})}{\pi (L_k,E_{\mathrm{p},k},z_{\mathrm{H}})},\end{array}$$(33)

where ${\{(L_k,E_{\text{p},k})\}}_{k=1}^{N_{\text{s}}}$ are a total of N_s samples from our analysis of the peak spectrum of GRB 170817A and ${\left\{{\theta }_{\text{v},j}\right\}}_{j=1}^{N_{\text{g}}}$ are N_g posterior samples of the probability in Eq. (22). Clearly, the term P(d_G17, d_HOST) eventually cancels out in Eq. (25) as before.

In the flux-limited sample analysis, the information on the viewing angle, luminosity and peak photon energy of GRB 170817A/GW170817 can still be used to condition the apparent structure parameters before the flux-limited sample analysis is performed, because the structure must be consistent with what has been observed in that event. The starting point is

$$\begin{array}{c}\hfill \begin{array}{cc}& P({\mathit{\lambda }}_{\mathrm{pop}}|d_{\mathrm{G}17},d_{\mathrm{HOST}})=\hfill \\ \hfill & {\displaystyle \int \int \int }P({\mathit{\lambda }}_{\mathrm{pop}}|L,E_{\mathrm{p}},{\theta }_{\mathrm{v}})P(L,E_{\mathrm{p}},{\theta }_{\mathrm{v}}|d_{\mathrm{G}17},d_{\mathrm{HOST}})\mathrm{d}L\mathrm{d}{E}_{\mathrm{p}}\mathrm{d}{\theta }_{\mathrm{v}}.\hfill \end{array}\end{array}$$(34)

Application of Bayes’ theorem gives

$$\begin{array}{c}\hfill P({\mathit{\lambda }}_{\mathrm{pop}}|L,E_{\mathrm{p}},{\theta }_{\mathrm{v}})=\frac{P(L,E_{\mathrm{p}},{\theta }_{\mathrm{v}}|{\mathit{\lambda }}_{\mathrm{pop}})\pi ({\mathit{\lambda }}_{\mathrm{pop}})}{\pi ({\theta }_{\mathrm{v}})\pi (L,E_{\mathrm{p}})},\end{array}$$(35)

but we have P(L, E_p, θ_v | λ_pop) = P(L, E_p | θ_v, λ_pop) sin θ_v and π(θ_v) = sin θ_v, which therefore leads to P(λ_pop | L, E_p, θ_v) = P(L, E_p | θ_v, λ_pop)π(λ_pop)/π(L, E_p). Substitution of this into Eq. (34) leads to

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill P({\mathit{\lambda }}_{\mathrm{pop}}|d_{\mathrm{G}17},d_{\mathrm{HOST}})& =\pi ({\mathit{\lambda }}_{\mathrm{pop}}){\displaystyle \int \int \int }P(L,E_{\mathrm{p}},{\theta }_{\mathrm{v}}|d_{\mathrm{G}17},d_{\mathrm{HOST}})\hfill \\ \hfill & \times \frac{P(L,E_{\mathrm{p}}|{\theta }_{\mathrm{v}},{\mathit{\lambda }}_{\mathrm{pop}})}{\pi (L,E_{\mathrm{p}})}\mathrm{d}L\mathrm{d}{E}_{\mathrm{p}}\mathrm{d}{\theta }_{\mathrm{v}}.\hfill \end{array}\end{array}$$(36)

This is similar to 𝒩_G17, but the redshift information is not used. This can be again approximated by Monte Carlo integration over samples drawn from the posterior P(L, E_p, θ_v | d_G17, z_H), namely

$$\begin{array}{c}\hfill P({\mathit{\lambda }}_{\mathrm{pop}}|d_{\mathrm{G}17},d_{\mathrm{HOST}})\sim \frac{\pi ({\mathit{\lambda }}_{\mathrm{pop}})}{N_{\mathrm{s}}{N}_{\mathrm{g}}}{\displaystyle \sum _{k=1}^{N_{\mathrm{s}}}\sum _{j=1}^{N_{\mathrm{g}}}}\frac{P(L_k,E_{\mathrm{p},k}|{\theta }_{\mathrm{v},j},{\mathit{\lambda }}_{\mathrm{pop}})}{\pi (L_k,E_{\mathrm{p},k},z_{\mathrm{H}})}.\end{array}$$(37)

This result can then be used as a prior in the analysis of the observer-frame and rest-frame samples. The final posterior on λ_pop then becomes

$$\begin{array}{c}\hfill \begin{array}{cc}& P_{\mathrm{flux}-\mathrm{limited}}({\mathit{\lambda }}_{\mathrm{pop}}|\mathit{d},d_{\mathrm{G}17},d_{\mathrm{HOST}})\sim \frac{\pi ({\mathit{\lambda }}_{\mathrm{pop}})}{P(\mathit{d})}\left({\displaystyle \prod _{i=1}^N}\frac{{\mathcal{N}}_i}{\mathcal{D}}\right)\hfill \\ \hfill & \times \frac{1}{N_{\mathrm{s}}{N}_{\mathrm{g}}}{\displaystyle \sum _{k=1}^{N_{\mathrm{s}}}\sum _{j=1}^{N_{\mathrm{g}}}}\frac{P(L_k,E_{\mathrm{p},k}|{\theta }_{\mathrm{v},j},{\mathit{\lambda }}_{\mathrm{pop}})}{\pi (L_k,E_{\mathrm{p},k},z_{\mathrm{H}})}.\hfill \end{array}\end{array}$$

Hence, the posterior takes a similar form as in the full-sample analysis case, the difference being in a missing P(z_H | λ_pop)/𝒟 factor.

2.5.4. Denominator 𝒟 for events in the three samples

The denominator 𝒟 in the above expressions represents the fraction of events in the population that pass the sample selection criteria (M19), that is, the ‘accessible’ events according to the selection effects model. For the observer-frame sample events, the denominator takes the form

$$\begin{array}{c}\hfill \mathcal{D}({\mathit{\lambda }}_{\mathrm{pop}})={\displaystyle \int \int \int }{P}_{\mathrm{pop}}(L,E_{\mathrm{p}},z|{\mathit{\lambda }}_{\mathrm{pop}})P_{\det ,\mathrm{GBM}}(L,E_{\mathrm{p}},z)\mathrm{d}L\mathrm{d}{E}_{\mathrm{p}}\mathrm{d}z.\end{array}$$(38)

For the rest-frame sample, it is

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \mathcal{D}({\mathit{\lambda }}_{\mathrm{pop}})& ={\displaystyle \int \int \int }{P}_{\mathrm{pop}}(L,E_{\mathrm{p}},z|{\mathit{\lambda }}_{\mathrm{pop}})\hfill \\ \hfill & \times P_{\det ,\mathrm{GBM}}(L,E_{\mathrm{p}},z)P_{\det ,\mathrm{BAT}}(L,E_{\mathrm{p}},z)\mathrm{d}L\mathrm{d}{E}_{\mathrm{p}}\mathrm{d}z.\hfill \end{array}\end{array}$$(39)

For the viewing-angle sample, it is a four-dimensional integral, namely

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \mathcal{D}({\mathit{\lambda }}_{\mathrm{pop}})& =P_{\mathrm{pop}}^{\star }(L,E_{\mathrm{p}},z,{\theta }_{\mathrm{v}}|{\mathit{\lambda }}_{\mathrm{pop}})\hfill \\ \hfill & \times P_{\det ,\mathrm{GBM}}(L,E_{\mathrm{p}},z)P_{\det ,\mathrm{GW}}({\theta }_{\mathrm{v}},z)\mathrm{d}L\mathrm{d}{E}_{\mathrm{p}}\mathrm{d}z\mathrm{d}{\theta }_{\mathrm{v}}.\hfill \end{array}\end{array}$$(40)

2.6. Choice of priors

Our Bayesian inference approach requires the definition of a prior probability density over the ‘hyper’ parameter space the λ_pop vector belongs to. We chose independent priors on all parameters, except for the pair (θ_c, θ_w) for which we enforce θ_w > θ_c, and hence

$$\begin{array}{c}\hfill \pi ({\mathit{\lambda }}_{\mathrm{pop}})=\pi ({\theta }_{\mathrm{c}},{\theta }_{\mathrm{w}})\pi (L_{\mathrm{c}}^{\star })\pi ({\alpha }_{\mathrm{L}})\dots \pi (z_{\mathrm{p}}).\end{array}$$(41)

As shown in Table 1, we selected mildly informative priors on most parameters, typically adopting a uniform or uniform-in-log prior over a relatively wide range that encompasses what we consider reasonable values. There are two exceptions: for the core half-opening angle θ_c and the break angle θ_w we adopted a prior that is uniform in the subtended solid angle, π(θ)∝ sin θ, since their role is that of setting the probability for a jet to be observed within the core or within the ‘wing’ where the break occurs; for the low-end cutoff of the core luminosity $L_{\text{c}}^{\star }$, we set a relatively high lower limit $L_{\text{c}}^{\star }$ ≥ 3 × 10⁵¹ erg s⁻¹, consistently with our quasi-universal structured jet scenario, where relatively low luminosities are produced by off-axis jets and not by intrinsically under-luminous jets.

2.7. Stochastic sampling of the population posterior

In order to realise our inference in practice, we sampled the posterior probability density on λ_pop in Eq. (25) using the publicly available python package emcee (Foreman-Mackey et al. 2013), which constitutes an efficient and flexibile implementation of the Goodman & Weare (2010) affine-invariant ensemble sampler⁵. Numerical integrals were performed using the trapezoidal rule over a four-dimensional grid with a resolution of 1000 linearly spaced points in θ_v between 0 and π/2 and 60 logarithmically spaced points in each of the remaining axes, in the domain L ∈ (10⁴⁴, 10⁵⁶) erg s⁻¹, E_p ∈ (10⁻¹, 10⁷) keV and z ∈ (10⁻³, 10¹). We ensured that this resolution was sufficient by comparing the value of the posterior at a number of points in the parameter space with those obtained with a higher resolution of 100 points over each axis, and found a negligible difference as long as A < 10 and σ_c ≳ 0.05 dex (our prior limits fall well within these requirements).

We set up the sampler to employ 14 × 2 = 28 walkers, and ran it for 10⁴ iterations for each analysis, for a total of 2.8 × 10⁵ posterior samples. The autocorrelation length of the resulting chains, averaged over all walkers, is around 650. A corner plot showing the density of the samples in the 14-dimensional parameter space is shown in Fig. A.1, while summary statistics for each parameter are reported in Table 3. The results presented in the next section are constructed using 1000 random posterior samples from these chains, after discarding the first half as burn-in⁶.

3. Results

Generally speaking, the full sample and flux-limited sample analyses yielded quite similar results, as demonstrated by the detailed posterior probability density distributions (Fig. A.1) and the summary in Table 3. The main notable differences in the full sample analysis versus the flux-limited sample analysis are a preference for a slightly larger on-axis SED peak photon energy E_p, c; a steeper slope, a, of the rate density evolution before peak; and a slightly better constrained redshift z_p of the peak of the rate density evolution. The posterior distributions from the two analyses show large overlaps for all parameters. In what follows, we present a thorough description of the results, further highlighting the differences in the results from the two analyses when relevant.

3.1. Apparent structure

First of all, we focus on the apparent jet structure. Figure 6 shows the constraint we obtained on the ‘average apparent jet structure’ $\overline{L}=L_{\mathrm{c}}^{\star }\ell$ (central panel) and ${\overline{E}}_{\mathrm{p}}=E_{\mathrm{p},\mathrm{c}}^{\star }\eta$ (bottom panel) from the two analyses, with the insets showing the posterior distribution on the core dispersion parameters A and σ_c. In this and the following figures, solid lines refer to the full-sample analysis and dashed lines to the flux-limited sample analysis. The small panel at the top shows the posterior probability density on the logarithm of θ_c (in grey) and of θ_w (in green). The figure shows that the preferred apparent structure features a narrow core of ∼​​2 − 3 deg, outside of which the luminosity falls off approximately as ${\theta }_{\mathrm{v}}^{-4.7}$ and the SED peak photon energy as ${\theta }_{\mathrm{v}}^{-2}$. The posterior on the transition angle θ_w rails against the upper boundary of its physical range, and the slopes β_L and β_Ep at larger viewing angles are not well constrained by the available data (see Table 3). This indicates that the data are consistent with an apparent structure described by a single power law outside the core: this may disfavour somewhat the presence of a distinct dissipation mechanism that dominates the gamma-ray emission at viewing angles in the range 20 ≲ θ_v/deg ≲ 30 (e.g. cocoon shock breakout as proposed by Gottlieb et al. 2018) in the majority of SGRBs, unless the apparent structure of the combined emission components follows the same power law. A caveat here is that there are most likely not more than a handful of events in our samples that are observed at such large viewing angles (see Sects. 4.2 and 4.4), and therefore this result must be taken with a grain of salt.

Fig. 6. Average apparent jet structure. The two larger panels show the constraints on the SGRB average apparent structure obtained from our analysis: solid lines show the medians of the projected posterior distributions of the jet structure functions at fixed θv from the full sample analysis (red: $\overline{L}=L_{\mathrm{c}}^{\star }\ell ({\theta }_{\mathrm{v}})$; blue: ${\overline{E}}_{\mathrm{p}}=E_{\mathrm{p},\mathrm{c}}^{\star }\eta ({\theta }_{\mathrm{v}})$), while dashed lines refer to the flux-limited sample (orange: $\overline{L}$; cyan: ${\overline{E}}_{\mathrm{p}}$). The shaded region around each line encompasses the corresponding 90% credible range at each fixed θv. Dotted lines show reference power law trends. The light green contours show the 50% and 90% credible regions of (L, θv) (middle panel) and (Ep, θv) (bottom panel) for GRB 170817A (Sect. 2.4.3). The insets show the posterior distributions on the parameters A and σc that determine the dispersion around the average structure within the population. The top smaller panel shows the posterior distributions on the logarithms of the transition angles, ln(θc) (grey) and ln(θw) (green), with the results for the full sample shown as filled areas, and those for the flux-limited sample shown as dashed lines.

The structure is by construction consistent with the GRB 170817A luminosity and E_p at the relevant viewing angle, as shown by the light green contours in the figure. These contours represent the 50% and 90% credible regions of (θ_v, L) and (θ_v, E_p) for GRB 170817A, constructed by combining the viewing angle information from the GW analysis conditioned on the host galaxy distance (Sect. 2.4.3) and our GRB 170817A spectral analysis at peak (Appendix B).

The uncertainty band on the typical luminosity $L_{\text{c}}^{\star }\mathcal{l}$ at θ_v ≪ θ_c is particularly narrow because of the chosen prior on $L_{\text{c}}^{\star }$ (see the discussion in Sect. 4.5). The larger uncertainty at θ_v ∼ θ_c arises because of the combined uncertainties on $L_{\text{c}}^{\star }$, θ_c and α_L.

The jet structure functions can also be projected onto the (L, E_p) plane: Fig. 7 shows the median $(\overline{L},{\overline{E}}_{\mathrm{p}})$ relation and the corresponding 90% credible region for the two analyses. In both cases, the bulk of the SGRBs with known redshift (shown by red crosses) are close to the upper-right end of the relation, which indicates that they are observed close to on-axis in the model (more precisely, close to the edge of the core; see also Sect. 4.2). Interestingly, the relation lies above most of the SGRBs in the known-redshift sub-sample, indicating that selection effects play a major role in how the plane is populated, according to the model. We expand on this in Sect. 4.3.

Fig. 7. Average apparent jet structure in the L − Ep plane. The solid turquoise line connects the points $(⟨\overline{L}({\theta }_{\mathrm{v}})⟩,⟨{\overline{E}}_{\mathrm{p}}({\theta }_{\mathrm{v}})⟩)$ at varying viewing angles, where ⟨ ⋅ ⟩ represents the median, from the full sample analysis. The turquoise shaded area encompasses the corresponding 90% credible range. Grey lines show 100 posterior samples of the lines $(\overline{L}({\theta }_{\mathrm{v}}),{\overline{E}}_{\mathrm{p}}({\theta }_{\mathrm{v}}))$ from the same analysis. The dashed green line and yellow area show the corresponding results for the flux-limited sample analysis. The crosses mark the positions of the Fermi/GBM SGRBs with known redshifts (all in red except for GRB 170817A, which is shown in blue). The inset shows the posterior probability density on the y parameter that sets the slope of the correlation between the on-axis luminosity and the on-axis SED peak photon energy, with the results of the two analyses colour coded as in the main panel.

The inset in Fig. 7 shows the posterior probability density on the y parameter that sets the correlation between the on-axis luminosity and peak SED photon energy. Despite our model allowing for such a correlation, the posterior is fully compatible with y = 0, that is, the absence of an intrinsic correlation between L_c and E_p, c.

3.2. Luminosity function

Our analysis does not directly constrain the local rate density, R₀, of SGRBs, because in our inference framework it is effectively only a normalisation factor as long as its prior is uniform in the logarithm (M19; Fishbach et al. 2018). In order to derive the local rate density implied by our results, we required the observed rate of events with p > p_lim to be equal to⁷R_obs(> p_lim, GBM) = N(> p_lim, GBM)/η_GBMT_GBM ≈ 27.6 yr⁻¹, where N(> p_lim, GBM) = 212 is the number of SGRBs with p > p_lim, GBM in our sample, η_GBM = 0.59 is a factor that corrects for the accessible field of view and the duty cycle of GBM (Burns et al. 2016), and T_GBM = 13 yr is the Fermi mission duration at the time of the last SGRB in the observer-frame sample. Given λ_pop, the local rate density is then

$$\begin{array}{c}\hfill R_0=\frac{R_{\mathrm{obs}}(>p_{\lim ,\mathrm{GBM}})}{{\displaystyle \int \int \int }\frac{\dot{\rho }/R_0}{1+z}\frac{\mathrm{d}V}{\mathrm{d}z}P(L,E_{\mathrm{p}}|{\mathit{\lambda }}_{\mathrm{pop}})P_{\det ,\mathrm{GBM}}(L,E_{\mathrm{p}},z)\mathrm{d}L\mathrm{d}{E}_{\mathrm{p}}\mathrm{d}z},\end{array}$$(42)

where P_det, GBM is the flux-limited form from Eq. (19) and $\dot{\rho }(z)/R_0$ is independent of R₀ (see Eq. (9)). We applied the above expression to our population posterior samples ${\left\{{\mathit{\lambda }}_{\text{pop},k}\right\}}_{k=1}^{N_{\text{p}}}$ to obtain an equal number of ${\left\{R_{0,k}\right\}}_{k=1}^{N_{\text{p}}}$ samples. In turn, this allowed us to construct samples of the posterior distribution of the luminosity function, ${\{R_{0,k}\varphi (L,{\lambda }_{\text{pop},\text{k}})\}}_{k=1}^{N_{\text{p}}}$, where ϕ(L, λ_pop) = L∫P(L, E_p | λ_pop) dE_p.

Figure 8 shows the luminosity function of SGRBs obtained in this way from the full-sample analysis (solid red) and the flux-limited sample analysis (dashed orange). The corresponding results from W15 and G16 (their fiducial model ‘a’) are shown in blue and grey, respectively. The shape of the result is somewhat in between the W15 and G16 results, but the normalisation of the high-luminosity end is in better agreement with W15 than G16. The high-luminosity end slope (1 − A ∼ 2) is in agreement with both benchmark results, while the flattening below the break is reminiscent of that found by G16, but translated to a luminosity that is lower by almost one order of magnitude. At luminosities below $L_{\text{c}}^{\star }$ ∼ 3 × 10⁵¹ erg s⁻¹, the slope of the luminosity function is ∼0.4, which is flatter than W15 (but not by as much as G16) and more similar to the older results from Guetta & Piran (2005, 2006). The low-luminosity end below ∼10⁴⁹ erg s⁻¹ remains highly uncertain.