Unpacking Merger Jets: A Bayesian Analysis of GW170817, GW190425 and Electromagnetic Observations of Short Gamma-Ray Bursts

The joint detection of both gravitational wave (GW) GW170817 (Abbott et al. 2017a), and counterpart short gamma-ray burst (sGRB) GRB 170817A (Goldstein et al. 2017; Savchenko et al. 2017), followed by the detection of kilonova AT 2017gfo (e.g., Evans et al. 2017; McCully et al. 2017) not only solidified the belief that many sGRBs are produced from the merger of binary neutron star (BNS) systems, but also began the era of GW multimessenger astronomy (Abbott et al. 2017d). The combination of both the electromagnetic (EM) and GW data gave insight into various problems that a detection through a single data channel could not provide, ranging from cosmology (Abbott et al. 2017c), the origin of the abundance of heavy elements in the Universe (Tanvir et al. 2017), tests of general relativity, and the speed of gravity (Abbott et al. 2017b) among others.

However the detection of the event not only provided answers but also provoked questions, as it was inferred through GW parameter inference that the event was exceptionally nearby at only 40 Mpc (Abbott et al. 2019), giving an observed isotropic luminosity of the event of 10⁴⁷ erg s⁻¹, three orders of magnitude lower than that observed for any other sGRB (Abbott et al. 2017b). It was also inferred through GW parameter inference that the event was viewed at a wide angle of 14°–41° from the central axis (Abbott et al. 2019). This led to the hypothesis that the jet of GRB 170817A exhibited some wide-angle structure to produce the observed flux, and that it may still have had a typically luminous central jet component. The long-duration observations of the event's afterglow across the EM spectrum provided further evidence for this claim (Troja et al. 2017, 2018a, 2020; Alexander et al. 2018; D'Avanzo et al. 2018; Lyman et al. 2018; Margutti et al. 2018; Mooley et al. 2018; Ruan et al. 2018).

While evidence for this wide-angle structure is provided by the multitude of afterglow observations, the functional form of the luminosity profile over viewing angle remains in question (e.g., Granot et al. 2017; Duffell et al. 2018; Gottlieb et al. 2018; Lamb & Kobayashi 2018; Mooley et al. 2018; Troja et al. 2018a; Beniamini et al. 2019; Fraija et al. 2019; Ioka & Nakamura 2019; Lamb et al. 2019, 2020; Salafia et al. 2019; Biscoveanu et al. 2020; Takahashi & Ioka 2021). Accurate modeling of this jet structure is important in both understanding the astrophysics of the event, and in preventing systematic biases in any multimessenger analysis that must make assumptions about the jet geometry (Lamb et al. 2021; Nakar & Piran 2021), which include constraints on BNS merger rate (Wanderman & Piran 2015), and the Hubble constant(Abbott et al. 2017c)—where it has been shown that without a reliable model of the jet structure the Hubble constant constraints are limited (Mastrogiovanni et al. 2021). An analogous observation to GRB 170817A and GW170817 was GRB 980425, which had a similar isotropic energy to GRB 170817A, and when considering the spectral peak energy of the prompt emission, it sits as an outlier from the Amati relation (Amati et al. 2002) for long GRBs (as does GRB 170817A for the equivalent relation for sGRBs). In the case of GRB 980425, the afterglow was not detected—there was a supernova, SN1998bw, at both optical and radio (Galama et al. 1998; Kulkarni et al. 1998); at the time, this was claimed to be distinct from the cosmological GRB population.

With the promise of future GW BNS merger detections (Abbott et al. 2020b), a jet structure model that best represents the data of joint detection events should be discerned (Hayes et al. 2020). However, current analyses are limited to the data provided by GW170817/GRB 170817A (including AT 2017gfo), GW BNS merger event GW190425 (Abbott et al. 2020a), as well as the population of sGRBs detected independently of GW detection (Lien et al. 2016; Poolakkil et al. 2021). Previous similar work has considered constraining jet structure models with joint GW and EM detections using a Bayesian analysis (Biscoveanu et al. 2020; Farah et al. 2020; Hayes et al. 2020). The historical rate of detected sGRBs and GWs has also been used in Bayesian analyses to constrain the jet structure (Williams et al. 2018; Sarin et al. 2022). Work by Mogushi et al. (2019) and Tan & Yu (2020) combine the detection rates of sGRBs and GWs with the prompt emission data of GRB 170817A and the parameter inference results of GW170817, along with assuming a luminosity function fitted by short GRB events with known redshift. The jet structure has been constrained given the weak gamma-ray emission detected by the INTEGRAL SPI-ACS detector coincident with the GW event GW190425 in the work by Saleem et al. (2020). In this work, we put forward a comprehensive Bayesian framework that combines the prompt emission data from GRB 170817A, GW parameter inference posteriors of GW170817 and GW190425, GW-informed BNS merger rate, as well as the detection rate of sGRBs by the Neil Gehrels Swift Observatory (Swift) detector. This analysis is then further combined with the information provided by a luminosity function fitted by sGRB events with known redshifts to provide tighter constraints, with the caveat of also introducing bias into the analysis. We provide parameter constraints and model comparison results between a classical top-hat (TH) jet structure and three different jet structures with wide-angle structuring: a Gaussian, power-law, and double-Gaussian jet. These results are presented alongside constraints on the intrinsic luminosity (when it is not fitted by the luminosity function), and the merger rate for both cases.

The physical model assumed is detailed in Section 2, before the analysis method and data are laid out in Section 3. In Section 4 we report the results of the analysis, both for when the fitted luminosity function is incorporated and when it is not. In Section 5 the implications of the results are discussed and a conclusion provided in Section 6.

The sGRB data consist of the observed T90 integrated flux F as well as the number of observed sGRBs N_EM. The average T90 integrated flux $\hat{F}$ is related to the isotropic equivalent luminosity L_iso at a given viewing angle θ_v , as well as a redshift z-dependent luminosity distance d_L and k-correction k:

$$\begin{eqnarray}&&\hat{F}=\displaystyle \frac{{L}_{\mathrm{iso}}({\theta }_{{\rm{v}}})}{4\pi {d}_{L}{\left(z\right)}^{2}k(z)}.\end{eqnarray} \tag{ 1 }$$

The mean number of sGRBs ${\hat{N}}_{\mathrm{EM}}$ observed by a detector within a duration T that covers an area of sky equal to ΔΩ depends on the redshift, viewing angle, and intrinsic luminosity Λ = {z, θ_v , L₀} through the number density:

$$\begin{eqnarray}&&{\hat{N}}_{\mathrm{EM}}=\int \displaystyle \frac{d{\hat{N}}_{\mathrm{EM}}}{d{\rm{\Lambda }}}\left[\int p({D}_{\mathrm{EM}}| {\rm{\Lambda }}){{dD}}_{\mathrm{EM}}\right]d{\rm{\Lambda }},\end{eqnarray} \tag{ 2 }$$

where D_EM denotes the selection effects of the detector, such that D_EM = 1 if a detection is made and D_EM = 0 otherwise.

We consider sGRBs detected by the Swift instrument, which has a detector response determined empirically in Lien et al. (2014) to fit

$$\begin{eqnarray}\int p({D}_{\mathrm{EM}}| {\rm{\Lambda }}){{dD}}_{\mathrm{EM}}=\left\{\begin{array}{ll}0 & \mathrm{if}\ \hat{F}\lt {F}_{\mathrm{thr}},\\ a\tfrac{b+c\hat{F}/{F}_{0}}{1+\hat{F}/{{dF}}_{0}} & \mathrm{else},\end{array}\right.\end{eqnarray} \tag{ 3 }$$

where a = 0.47, b = −0.05, c = 1.46, d = 1.45, F₀ = 6 × 10⁻⁶ erg s⁻¹ cm⁻², and F_thr = 5.5 × 10⁻⁹ erg s⁻¹ cm⁻².

The relation between the number density and the physical parameters of Λ is

$$\begin{eqnarray}&&\displaystyle \frac{d{\hat{N}}_{\mathrm{EM}}}{d{\rm{\Lambda }}}=T\displaystyle \frac{{\rm{\Delta }}{\rm{\Omega }}}{8\pi }\displaystyle \frac{{{ \mathcal R }}_{\mathrm{GRB}}(z)}{1+z}\displaystyle \frac{{dV}(z)}{{dz}}p({L}_{0}| {\rm{\Sigma }})\sin {\theta }_{v},\end{eqnarray} \tag{ 4 }$$

where V(z) is the co-moving volume, ${{ \mathcal R }}_{\mathrm{GRB}}$ is the rate of sGRBs and p(L₀∣Σ) is the intrinsic luminosity function given the hyperparameter Σ.

2.1. Short Gamma-Ray Burst Rate

The rate of sGRBs is assumed to be in the form

$$\begin{eqnarray}&&{{ \mathcal R }}_{\mathrm{GRB}}(z)={R}_{\mathrm{BNS}}{R}_{\mathrm{GRB}}(z),\end{eqnarray} \tag{ 5 }$$

where R_BNS is the local rate of BNS mergers and R_GRB(z) is defined so that R_GRB(0) = 1. This assumes that every BNS merger results in an sGRB, and that the number of sGRBs produced by neutron star–black hole mergers is negligible.

The form of R_GRB(z) can be assumed to follow the star formation rate R_*(z) convolved with the probability distribution of the delay time between the system formation and the eventual merger that leads to the sGRB P(t) (Wanderman & Piran 2015):

$$\begin{eqnarray}&&{R}_{\mathrm{GRB}}(z)\propto {\int }_{{t}_{\min }}^{T(\infty )-T(z)}\displaystyle \frac{{R}_{* }({z}_{* })}{1+{z}_{* }}P(t){dt},\end{eqnarray} \tag{ 6 }$$

where z_* = z(T(z) + t) is the redshift when the system was formed, T(z) is the look-back time and ${t}_{\min }$ is the minimum delay time. This minimum delay time is set to 20 Myr and P(t) ∝ 1/t according to Guetta & Piran (2006). The star formation rate is assumed to be of the form (Cole et al. 2001)

$$\begin{eqnarray}&&{R}_{* }(z)\propto \displaystyle \frac{a+{bz}}{1+{\left(z/c\right)}^{d}}H(z),\end{eqnarray} \tag{ 7 }$$

where the parameter values are taken from Hopkins & Beacom (2006) to be a = 0.017, b = 0.13, c = 3.3, and d = 5.3.

2.2. Cosmology

A flat, vacuum-dominated Universe is assumed. The co-moving volume distribution over redshift is defined as

$$\begin{eqnarray}&&\displaystyle \frac{{dV}(z)}{{dz}}=4\pi \displaystyle \frac{c}{H(z)}{\left(\displaystyle \frac{{d}_{L}(z)}{1+z}\right)}^{2}.\end{eqnarray} \tag{ 8 }$$

For a flat cosmology, the luminosity distance is related to the redshift by

$$\begin{eqnarray}&&{d}_{L}(z)=(1+z)\displaystyle \frac{c}{{H}_{0}}{\int }_{0}^{z}\displaystyle \frac{{H}_{0}}{H(z^{\prime} )}{dz}^{\prime} ,\end{eqnarray} \tag{ 9 }$$

where H₀ is the Hubble constant and H(z) is equal to

$$\begin{eqnarray}&&H(z)={H}_{0}\sqrt{{{\rm{\Omega }}}_{m}{\left(1+z\right)}^{3}+{{\rm{\Omega }}}_{{\rm{\Lambda }}}}.\end{eqnarray} \tag{ 10 }$$

Here Ω_m = 0.308 and Ω_Λ = 0.692 are the matter density and dark energy density, respectively (Hogg 1999), with values taken from Adam et al. (2016) along with H₀ = 67.8 km s⁻¹ Mpc⁻¹.

The look-back time, defined as the time between when a source emits light at redshift z and the time it is detected, is then

$$\begin{eqnarray}&&T(z)=\displaystyle \frac{1}{{H}_{0}}{\int }_{0}^{z}\displaystyle \frac{{H}_{0}}{(1+z^{\prime} )H(z^{\prime} )}{dz}^{\prime} .\end{eqnarray} \tag{ 11 }$$

For a flat, vacuum-dominated Universe, the inverse function has an analytical expression (Petrillo et al. 2013)

$$\begin{eqnarray}&&z(T)={\left(\displaystyle \frac{{{\rm{\Omega }}}_{{\rm{\Lambda }}}}{{{\rm{\Omega }}}_{m}}\right)}^{1/3}{\left[{\left(\displaystyle \frac{1+W(T)}{1-W(T)}\right)}^{2}-1\right]}^{1/3}-1,\end{eqnarray} \tag{ 12 }$$

where

$$\begin{eqnarray}&&W(T)=\exp \left[\mathrm{ln}\left(\displaystyle \frac{1+\sqrt{{{\rm{\Omega }}}_{{\rm{\Lambda }}}}}{1-\sqrt{{{\rm{\Omega }}}_{{\rm{\Lambda }}}}}\right)-3{H}_{0}\sqrt{{{\rm{\Omega }}}_{{\rm{\Lambda }}}}T\right].\end{eqnarray} \tag{ 13 }$$

The k-correction accounts for the cosmological redshifting in the intrinsic sGRB spectrum with respect to the detector's spectrum (Bloom et al. 2001):

$$\begin{eqnarray}&&k(z)=\displaystyle \frac{{\int }_{{\nu }_{s,1}/(1+z)}^{{\nu }_{s,2}/(1+z)}\nu f(\nu )d\nu }{{\int }_{{\nu }_{1}}^{{\nu }_{2}}\nu f(\nu )d\nu },\end{eqnarray} \tag{ 14 }$$

where the source spectrum is bounded between ν_s,1 = 1 keV and ν_s,2 = 10 MeV. The energy range assumed for the Swift detector is ν₁ = 15 keV and ν₂ = 150 keV while ν₁ = 50 keV and ν₂ = 300 keV are assumed for the Fermi detector. Here we assume the form of f(ν) follows the Band function described in Band et al. (1993) with a peak spectral energy of 800 keV, and power-law indices of −1 and −2.25 for the low- and high-energy regimes, respectively. We note that while a comptonized spectral model is more often preferred for sGRBs (e.g., see Poolakkil et al. 2021), the peak spectral energy of sGRBs typically exceeds the high-energy limit of the Swift detector and therefore the model choice makes little impact.

2.3. Intrinsic and Isotropic Equivalent Luminosity

The luminosity structure of a gamma-ray burst is defined to be

$$\begin{eqnarray}&&L(\theta )={L}_{0}{y}_{L}(\theta ),\end{eqnarray} \tag{ 15 }$$

where L₀ is the intrinsic luminosity at θ = 0. It is assumed that the distribution of intrinsic luminosity L₀ follows a Schechter function:

$$\begin{eqnarray}&&p({L}_{0}| {\rm{\Sigma }}=\{{L}_{0}^{* },\gamma \})\propto {\left(\displaystyle \frac{{L}_{0}}{{L}_{0}^{* }}\right)}^{-\gamma }{e}^{-{L}_{0}/{L}_{0}^{* }},\end{eqnarray} \tag{ 16 }$$

for ${L}_{\min }\leqslant {L}_{0}\leqslant {L}_{\max }$, where ${L}_{\min }={10}^{-3}{L}_{0}^{* }$ and ${L}_{\max }\,={10}^{2}{L}_{0}^{* }$.

Similarly, the Lorentz factor's dependence over angle follows

$$\begin{eqnarray}&&{\rm{\Gamma }}(\theta )=({{\rm{\Gamma }}}_{0}-1){y}_{{\rm{\Gamma }}}(\theta )+1,\end{eqnarray} \tag{ 17 }$$

with Γ₀ being the Lorentz factor of the jet at θ = 0. Given these definitions, both y_L and y_Γ are defined to equal 1 at θ = 0.

The Lorentz factor determines the degree of relativistic beaming, which for the luminosity is governed by

$$\begin{eqnarray}&&\begin{array}{l}{{ \mathcal B }}_{L}({\theta }_{v},\theta )=\displaystyle \frac{1}{4{\rm{\Gamma }}{\left(\theta \right)}^{6}{\left({A}^{2}-{B}^{2}\right)}^{2}}\\ \ \ \times \ \left(5{\left(\displaystyle \frac{A}{\sqrt{{A}^{2}-{B}^{2}}}\right)}^{3}+3\left(\displaystyle \frac{A}{\sqrt{{A}^{2}-{B}^{2}}}\right)\right),\end{array}\end{eqnarray} \tag{ 18 }$$

with $A=1-\beta \cos \theta \cos {\theta }_{v}$ and $B=-\beta \sin \theta \sin {\theta }_{v}$ where $\beta =\sqrt{1-{\rm{\Gamma }}{\left(\theta \right)}^{-2}}$.

The apparent isotropic equivalent luminosity, for an observer at θ_v from the jet axis, can then be related to the intrinsic luminosity via the beaming function by combining Equation (15) and Equation (18):

$$\begin{eqnarray}&&{L}_{\mathrm{iso}}({\theta }_{v})={L}_{0}{\int }_{0}^{{\theta }_{j}}\displaystyle \frac{1}{2}{ \mathcal B }({\theta }_{v},\theta ){y}_{L}(\theta )\sin \theta d\theta ,\end{eqnarray} \tag{ 19 }$$

where θ_j is the maximum angle for which the beaming of gamma-rays occurs. A conservative maximum outer jet angle for the emission of gamma-rays is approximated by considering scattering by electrons accompanying baryons within the jet. The condition is given by Matsumoto et al. (2019) and Lamb et al. (2022):

$$\begin{eqnarray}&&{\rm{\Gamma }}({\theta }_{j})\simeq 19.1{\left(\displaystyle \frac{L({\theta }_{j})}{{10}^{51}\,\mathrm{erg}\,{{\rm{s}}}^{-1}}\right)}^{1/6},\end{eqnarray} \tag{ 20 }$$

Beyond this limit, θ > θ_j , the jet becomes opaque to gamma-rays.

2.4. Jet Structures

The implications of structuring within compact stellar merger jets for the EM counterparts from GW detected systems has been highlighted in the literature (Lazzati et al. 2017b; Lamb & Kobayashi 2017; Kathirgamaraju et al. 2018; Beniamini et al. 2020); here, we choose a sample of fiducial jet structure models that are representative of the literature diversity.

The TH jet is the simplest structure, where the beam is uniform until the jet opening angle θ_j where the jet sharply cuts off:

$$\begin{eqnarray}\begin{array}{rcl}{y}_{L}(\theta ) & = & \left\{\begin{array}{ll}1 & \ \mathrm{if}\ 0\leqslant \theta \leqslant {\theta }_{j},\\ 0 & \ \mathrm{else}.\end{array}\right.\\ {y}_{{\rm{\Gamma }}}(\theta ) & = & \left\{\begin{array}{ll}1 & \ \mathrm{if}\ 0\leqslant \theta \leqslant {\theta }_{j},\\ 0 & \ \mathrm{else}.\end{array}\right.\end{array}\end{eqnarray} \tag{ 21 }$$

We note that the condition expressed in Equation (20) is not enforced for this case, as Γ₀ is above the Equation (20) limit at all points within the jet.

Wide-angle structure can be introduced with a Gaussian jet (GJ) structure, described by a single-width parameter θ_σ (e.g., Rossi et al. 2002, 2004; Zhang & Meszaros 2002; Kumar & Granot 2003):

$$\begin{eqnarray}&&{y}_{L}(\theta )={e}^{-\tfrac{1}{2}{\left(\tfrac{\theta }{{\theta }_{\sigma }}\right)}^{2}},\quad {y}_{{\rm{\Gamma }}}(\theta )={e}^{-\tfrac{1}{2}{\left(\tfrac{\theta }{{\theta }_{\sigma }}\right)}^{2}}.\end{eqnarray} \tag{ 22 }$$

An alternative to the Gaussian profile has the wide-angle emission expressed as a three-parameter power-law (PL) jet structure (e.g., Kumar & Granot 2003; Rossi et al. 2004; Zhang et al. 2004), where the jet can be described by some uniform core out to width θ_c , and then the intrinsic luminosity structure falls off at wide angles according to power s and the Lorentz factor with a:

$$\begin{eqnarray}\begin{array}{rcl}{y}_{L}(\theta ) & = & \left\{\begin{array}{ll}1 & \ \mathrm{if}\ 0\leqslant \theta \leqslant {\theta }_{c},\\ {\left(\tfrac{\theta }{{\theta }_{c}}\right)}^{-s} & \ \mathrm{else}.\end{array}\right.\\ {y}_{{\rm{\Gamma }}}(\theta ) & = & \left\{\begin{array}{ll}1 & \ \mathrm{if}\ 0\leqslant \theta \leqslant {\theta }_{c},\\ {\left(\tfrac{\theta }{{\theta }_{c}}\right)}^{-a} & \ \mathrm{else}.\end{array}\right.\end{array}\end{eqnarray} \tag{ 23 }$$

Finally, let us consider a two-component, or double-Gaussian jet (DG), with emission from both an inner core described by a Gaussian structure of width θ_in and an outer cocoon described by width θ_out (Salafia et al. 2020):

$$\begin{eqnarray}\begin{array}{rcl}{y}_{L}(\theta ) & = & (1-{ \mathcal C }){e}^{-\tfrac{1}{2}{\left(\tfrac{\theta }{{\theta }_{\mathrm{in}}}\right)}^{2}}+{ \mathcal C }{e}^{-\tfrac{1}{2}{\left(\tfrac{\theta }{{\theta }_{\mathrm{out}}}\right)}^{2}},\\ {y}_{{\rm{\Gamma }}}(\theta ) & = & \displaystyle \frac{(1-{ \mathcal C }){e}^{-\tfrac{1}{2}{\left(\tfrac{\theta }{{\theta }_{\mathrm{in}}}\right)}^{2}}+{ \mathcal C }{e}^{-\tfrac{1}{2}{\left(\tfrac{\theta }{{\theta }_{\mathrm{out}}}\right)}^{2}}}{(1-\tfrac{{ \mathcal C }}{{ \mathcal A }}){e}^{-\tfrac{1}{2}{\left(\tfrac{\theta }{{\theta }_{\mathrm{in}}}\right)}^{2}}+\tfrac{{ \mathcal C }}{{ \mathcal A }}{e}^{-\tfrac{1}{2}{\left(\tfrac{\theta }{{\theta }_{\mathrm{out}}}\right)}^{2}}}.\end{array}\end{eqnarray} \tag{ 24 }$$

The luminosity of the outer cocoon is equal to ${ \mathcal C }{L}_{0}$ and the Lorentz factor ${ \mathcal A }({{\rm{\Gamma }}}_{0}-1)+1$.

We do not consider hollow-cone jet structure models in our study (see, e.g., Nathanail et al. 2021; Takahashi & Ioka 2021); we expect that the combination of our intrinsic luminosity distribution (Equation (16)) and the beaming (Equation (19)) will wash out the effect of any hollow-cone structuring within the core. For this study, the structure outside of the jet's core is the critical component.

Constraints are placed on model parameters λ of a model M when given data ${ \mathcal D }$ in Bayesian data analysis by determining the posterior distribution using the Bayes theorem:

$$\begin{eqnarray}&&p(\lambda | { \mathcal D },M)=\displaystyle \frac{{ \mathcal L }({ \mathcal D }| \lambda )\pi (\lambda )}{p({ \mathcal D }| M)},\end{eqnarray} \tag{ 25 }$$

where ${ \mathcal L }$ is the likelihood, π is the prior, and the normalization term $p({ \mathcal D }| M)$ is the evidence. Consider comparing two models M₁ and M₂ when given data ${ \mathcal D }$. In the context of Bayesian data analysis, the statistic used to compare two models is the posterior odds, defined as

$$\begin{eqnarray}&&{{ \mathcal O }}_{12}=\displaystyle \frac{p({M}_{1}| { \mathcal D })}{p({M}_{2}| { \mathcal D })}=\displaystyle \frac{p({M}_{1})}{p({M}_{2})}\displaystyle \frac{p({ \mathcal D }| {M}_{1})}{p({ \mathcal D }| {M}_{2})}.\end{eqnarray} \tag{ 26 }$$

Normally we are interested in cases where the a priori probability of either model being correct is comparable, and therefore the posterior odds is dominated by the Bayes factor:

$$\begin{eqnarray}&&{{ \mathcal B }}_{12}=\displaystyle \frac{p({ \mathcal D }| {M}_{1})}{p({ \mathcal D }| {M}_{2})},\end{eqnarray} \tag{ 27 }$$

which quantifies the contribution to the posterior odds given by the data ${ \mathcal D }$. A value of $\mathrm{ln}{{ \mathcal B }}_{12}\gt 0$ favors M₁, while $\mathrm{ln}{{ \mathcal B }}_{12}\lt 0$ favors M₂.

The analysis is performed by applying the model described in Section 2 with an assumed jet structure from Section 2.4 given both GW and sGRB prompt emission data.

Table 1 lists the notation used in the following section. The data can be split into that produced by a GW-triggered event, denoted with the subscript "GW," and that produced from an EM trigger, denoted with the subscript "EM." Figure 1 displays the relations between the observable variables that are represented by shaded nodes and the latent variables represented by clear nodes through the edges. The part of the network contained within the panel is dependent on each GW event, while those outside of the panel describe the population of events.

Table 2. Summary of the Data Used in the Analysis

Data set	Data
${{ \mathcal D }}_{170817}$	{x170817, F170817}
${{ \mathcal D }}_{190425}$	{x190425, F190425}
${{ \mathcal D }}_{R}$	{NEM, NGW}
${{ \mathcal D }}_{170817+R}$	$\{{{ \mathcal D }}_{170817},{{ \mathcal D }}_{R}\}$
${{ \mathcal D }}_{\mathrm{all}}$	$\{{{ \mathcal D }}_{170817},{{ \mathcal D }}_{190425},{{ \mathcal D }}_{R}\}$

Download table as:  ASCIITypeset image

The data from the N_GW GW-triggered events consist of the GW strain x_GW and the flux of the counterpart F_GW. The GW-triggered events may not necessarily require a counterpart to be considered for the analysis. If the sky localization of the source coincides with the sky coverage of GRB detectors then we can assume that it was not detected due to its distance and orientation to us. The current events that meet these criteria are both GW170817 with GRB 170817A as well as GW190425 and the nondetection of its counterpart, under the assumption that an sGRB was produced, given the Fermi detector covered 50% of the sky localization, INTEGRAL SPI-ACS had a near-complete sky coverage (Martin-Carillo et al. 2019), and Konus-Wind covered the entire sky (Hosseinzadeh et al. 2019).

The EM-triggered events are simply the number of sGRB detections that Swift made within a 9.8 yr operational period (after assuming a 78% duty cycle over a 12.6 yr observation (Mogushi et al. 2019)), N_EM.

The likelihood can be decoupled into two terms, one of which considers GW-triggered events and the other EM-triggered:

$$\begin{eqnarray}&&{ \mathcal L }\,=\,{{ \mathcal L }}_{\mathrm{EM}}{{ \mathcal L }}_{\mathrm{GW}}.\end{eqnarray} \tag{ 28 }$$

The likelihood of the EM-triggered events is a Poisson distribution with a mean given in Equation (2):

$$\begin{eqnarray}&&{{ \mathcal L }}_{\mathrm{EM}}=\displaystyle \frac{{\hat{N}}_{\mathrm{EM}}{\left({\rm{\Theta }},{R}_{\mathrm{BNS}},{\rm{\Sigma }}\right)}^{{N}_{\mathrm{EM}}}{e}^{-{\hat{N}}_{\mathrm{EM}}({\rm{\Theta }},{R}_{\mathrm{BNS}},{\rm{\Sigma }})}}{{N}_{\mathrm{EM}}!}.\end{eqnarray} \tag{ 29 }$$

The mean is evaluated over a regular grid of shape $(z,{(2\theta /\pi )}^{1/3},{\mathrm{log}}_{10}{L}_{0})=(\times 50,\times 100,\times 1000)$. The angular grid points were chosen to be distributed over a power law so as to populate low-θ areas of the parameter space with grid points, while also maintaining a relatively high density of points at wider angles where emission from some jet structures is still significant.

The GW-triggered events likelihood is the product of each of the N_GW events:

$$\begin{eqnarray}&&{{ \mathcal L }}_{\mathrm{GW}}\propto \displaystyle \prod _{i=1}^{{N}_{\mathrm{GW}}}\displaystyle \sum _{j=1}^{S}p({F}_{i,\mathrm{GW}}| {\rm{\Theta }},{{\rm{\Phi }}}_{i,j},{L}_{j,0}),\end{eqnarray} \tag{ 30 }$$

where S samples are taken of Φ_i,j and L_j,0 from p(Φ_i , L₀∣x_i,GW, Σ). The parameters can be sampled from separate distributions p(Φ_i ∣x_i,GW) and p(L₀∣Σ), respectively, where p(Φ_i ∣x_i,GW) are samples from the posteriors produced from GW parameter estimation for each event. The likelihood of the prompt emission of the GW-triggered events is assumed to be a Gaussian distribution of width σ_F about a mean described in Equation (1):

$$\begin{eqnarray}&&p({F}_{\mathrm{GW}}| {\rm{\Theta }},{\rm{\Phi }},{L}_{0})=\displaystyle \frac{1}{\sqrt{2\pi {\sigma }_{F}^{2}}}\exp \left(-\displaystyle \frac{{\left({F}_{\mathrm{GW}}-F\right)}^{2}}{2{\sigma }_{F}^{2}}\right).\end{eqnarray} \tag{ 31 }$$

The priors for the model are specified in Table 3 for each of the GRB rate, luminosity function, and jet structure parameters. A normal distribution is denoted ${ \mathcal N }(\mu ,\sigma )$ with a mean of μ and standard deviation of σ, a uniform distribution as ${ \mathcal U }(A,B)$ with lower bound of A and upper bound of B, a Gamma distribution and inverse Gamma distribution as Γ(α, β) and Γ⁻¹(α, β) with a shape of α and a scale of β (not to be confused with the Lorentz factor Γ(θ)). We assume that every BNS merger results in a GRB so that _BNS = 1.

Table 3. Assumed Prior Distributions for Each Parameter

Model	Parameter	Prior
	${\mathrm{log}}_{10}{R}_{\mathrm{BNS}}^{\prime}$	${ \mathcal N }(-6.6,0.77)$
	${\mathrm{log}}_{10}{L}_{0}^{* }$'	${ \mathcal N }(51.6,1)$
	γ	${ \mathcal U }(0,1)$
	Γ0	Γ−1(2, 2.5 × 10−3)
TH	θj	${ \mathcal U }(0,\pi /2)$
GJ	θσ	${ \mathcal U }(0,\pi /2)$
PL	θc	${ \mathcal U }(0,\pi /2)$
	s	Γ(2, 4)
	a	Γ(2, 1)
DG	θin	${ \mathcal U }(0,\pi /2)$
	θout	${ \mathcal U }(0,{\theta }_{\mathrm{out}})$
	${\mathrm{log}}_{10}{ \mathcal C }$	${ \mathcal U }(-6,0)$
	${\mathrm{log}}_{10}{ \mathcal A }$	${ \mathcal U }(-6,0)$

Notes. The analysis assumes one jet structure out of the top-hat (TH), Gaussian (GJ), power-law (PL), and double-Gaussian (DG) models. Some parameters are made unit-less so that ${R}_{\mathrm{BNS}}^{\prime} ={R}_{\mathrm{BNS}}/{\mathrm{Mpc}}^{-3}\,{\mathrm{yr}}^{-1}$ and ${L}_{0}* ^{\prime}$ $\,=\,{L}_{0}^{* }/\mathrm{erg}\,{{\rm{s}}}^{-1}$.

Download table as:  ASCIITypeset image

The T90 integrated flux of GRB 170817A in the Fermi detector's 50–300 keV band is set at F₁₇₀₈₁₇ = 1.4 × 10⁻⁷ erg s⁻¹ cm⁻² with an uncertainty of σ₁₇₀₈₁₇ = 3.64 × 10⁻⁸ erg s⁻¹ cm⁻² (Goldstein et al. 2017). For the unobserved counterpart of GW190425, it is assumed that the T90 integrated flux takes a value of zero with an uncertainty of σ₁₉₀₄₂₅ = 10⁻⁸ erg s⁻¹ cm⁻² as a conservative upper bound to the Fermi detector's detection threshold (Tan & Yu 2020). We note that the detector sensitivity limit of INTEGRAL SPI-ACS approaches that of the Fermi detector, which had a more complete sky coverage in comparison. The distance and viewing angle posteriors of GW170817 and GW190425 are each represented by 500 samples taken from the their respective parameter estimation data releases. In this work we consider an observing period of approximately 9.8 yr by the Swift detector in which it observed N_EM = 107 sGRBs as recorded by Lien et al. (2016), given a sky coverage of ΔΩ = 0.1. The log prior on the rate of BNS mergers of ${ \mathcal N }(-6.6,0.77)$ is roughly chosen to reflect the constraints imposed to the rates by GWTC-2 (Abbott et al. 2020c). The log prior on ${L}_{0}^{* }$ is centered around the fitted value taken from Mogushi et al. (2019) with a standard deviation set to span one order of magnitude. This is chosen to reflect some prior information in the allowed luminosity from prior observations, but with a width to allow for flexibility into higher- or lower-luminosity regimes.

Posteriors and Bayes factors are calculated from Equations (25) and (27) by assigning ${ \mathcal D }$ and λ as the variables in Table 1. We collect the data into three sets: one only given the number of Swift detections N_EM and GW detections N_GW called ${{ \mathcal D }}_{R}$, another with the combined GW170817 GW x₁₇₀₈₁₇ and EM data F₁₇₀₈₁₇ called ${{ \mathcal D }}_{170817}$, and the other with GW190425 GW data x₁₉₀₄₂₅ and the flux from the nondetection F₁₉₀₄₂₅ called ${{ \mathcal D }}_{190425}$. The analysis is performed on five combinations of these data sets as shown in Table 2: ${{ \mathcal D }}_{170817}=\{{x}_{170817},{F}_{170817}\}$, ${{ \mathcal D }}_{190425}=\{{x}_{190425},{F}_{190425}\}$, ${{ \mathcal D }}_{R}\,=\{{N}_{\mathrm{EM}},{N}_{\mathrm{GW}}\}$, ${{ \mathcal D }}_{170817+R}=\{{{ \mathcal D }}_{170817},{{ \mathcal D }}_{R}\}$, and ${{ \mathcal D }}_{\mathrm{all}}\,=\{{{ \mathcal D }}_{170817},{{ \mathcal D }}_{190425},{{ \mathcal D }}_{R}\}$. We allow that $\lambda =\{{\rm{\Sigma }},{{ \mathcal R }}_{0},{{\rm{\Theta }}}_{M}\}$, where Θ_M are the jet structure model parameters dependent on jet structure model M. The three analyses are repeated for each of the jet structure models: M = TH, GJ, PL, and DG.

The posterior samples and evidence for each case are calculated via the nested sampling algorithm Nessai, which utilizes machine-learning techniques to drastically reduce the number of evaluations of the expensive likelihood function (Williams et al. 2021).

The analysis that is described in the previous section is applied to all three sets of data. The full corner plots for each jet structure model are shown in the Appendix, where Figure A1 shows the results for the TH jet, Figure A2 the GJ, Figure A3 the PL jet and Figure A4 the DG jet structure model. The posteriors for each of the data sets are overlaid upon one another where ${{ \mathcal D }}_{170817}$ is shown in red, ${{ \mathcal D }}_{190425}$ in orange, ${{ \mathcal D }}_{R}$ in violet, ${{ \mathcal D }}_{170817+R}$ in blue, and ${{ \mathcal D }}_{\mathrm{all}}$ in black. The log evidence $\mathrm{ln}p({ \mathcal D }| M)$ that corresponds to each posterior is shown in Table A1 for each of the five data sets over the four jet structure models. The log Bayes factors between the different models given the same data set can simply be calculated by taking the difference between entries of the same row.

Much of the discussion in this section concerns the posterior constraints when all of the data are considered ${ \mathcal D }={{ \mathcal D }}_{\mathrm{all}}$; however, the outcomes given the other subsets are considered to explain these results and provide further insight.

The log Bayes factors between the jet structure models are given as the left-hand entries of Table 4 when given the ${{ \mathcal D }}_{\mathrm{all}}$ data. A positive value indicates that the data support the model of the row while a negative value supports the column model. Evidently the TH model is less favourable than the models with wide-angle jet structuring, with Bayes factors of −0.54, −0.46, and −0.55 between it and the GJ, PL, and DG jet structures, respectively. The log Bayes factors between the PL, GJ, and DG is slight, with only insignificant evidence in favor of the GJ and DG model of log Bayes factors of less than 0.1, and negligibly small log Bayes factors between the two.

Zoom In Zoom Out Reset image size

Table 4. The Log Bayes Factor $\mathrm{ln}{ \mathcal B }$ between Each Model When Given ${{ \mathcal D }}_{\mathrm{all}}$ without or with the Fitted Luminosity Function

$\mathrm{ln}{{ \mathcal B }}_{\mathrm{row},\mathrm{col}}$	Top-hat	Gaussian	Power Law	Double Gaussian
Top-hat	0	−0.54/−0.91	−0.46/−0.45	−0.55/−0.64
Gaussian	0.54/0.91	0	0.08/0.46	−0.01/0.27
Power law	0.46/0.45	−0.08/−0.46	0	−0.09/−0.19
Double Gaussian	0.55/0.64	0.01/−0.27	0.09/0.19	0

Notes. A positive value is evidence toward the model of the row while a negative is evidence that favors the column model. Slight evidence is provided for models with wide-angle structuring over the TH model while little evidence distinguishes between the PL, GJ, and DG models without the fitted luminosity function. When the fitted luminosity function is included, slight evidence is further provided in favor of the GJ over the DG and PL jet structures while the TH remains least favored. All values can be assumed to have uncertainties of ±0.05.

Download table as:  ASCIITypeset image

The constraints on the rate of BNS mergers when given ${{ \mathcal D }}_{\mathrm{all}}$ are shown in Figure 2(a) for the four jet structure models. These constraints take the form of posterior distributions that are represented in the violin plots, where the outermost solid vertical lines indicate the minimum and maximum sample value while the fill in between represents the probability density. The 90% narrowest credible intervals are shown by the vertical dashed lines which enclose the median indicated by the middle solid line. The posterior distributions are compared to samples from the prior distribution at the top of the figure. For all cases, the posterior places tighter constraints on the merger rate than the prior distribution. The cases with wide-angle structuring (GJ, Pl, and DG models) produce similar posterior distributions to one another, centered around a value of ∼10⁻⁷ Mpc⁻³ yr⁻¹ consistent with the mean of the prior. The GJ structure produces the narrowest constraints with a 90% credible interval of ${\mathrm{log}}_{10}{R}_{\mathrm{BNS}}/{\mathrm{Mpc}}^{-3}\,{\mathrm{yr}}^{-1}=-{6.7}_{-0.8}^{+0.7}$, compared to the PL and DG models of $-{6.8}_{-0.9}^{+0.7}$ and $-{6.6}_{-0.9}^{+0.9}$, respectively. The TH model favors lower rates of BNS mergers, and even pushes the lower bound on the 90% credible interval to lower values of that of the prior, constraining it between ${\mathrm{log}}_{10}{R}_{\mathrm{BNS}}/{\mathrm{Mpc}}^{-3}\,{\mathrm{yr}}^{-1}=-{7}_{-1.0}^{+0.9}$.

Zoom In Zoom Out Reset image size

The median intrinsic luminosity posteriors determined for each model when given ${{ \mathcal D }}_{\mathrm{all}}$ is shown in Figure 2(b). The mean intrinsic luminosity ${\hat{L}}_{0}$ is determined by drawing L₀ and γ from the respective posterior distribution and then drawing 1000 samples from the corresponding Schechter function of Equation (16) before finding the ensemble median. This process is then repeated for 2000 median intrinsic luminosity samples. These posteriors take a form similar to the rates posteriors in Figure 2 as violin plots where the shaded probability density is contained within the outermost maximum and minimum values indicated by the solid vertical lines, while the median is marked by the middle solid line. The median is enclosed by the narrowest 90% credible intervals displayed as dashed vertical lines. A distribution of prior samples of ${\hat{L}}_{0}$ is also plotted at the top of the figure, which is determined by sampling from the individual ${\mathrm{log}}_{10}{L}_{0}^{* }$ and γ priors defined in Table 3. The TH jet structure resembles the prior in width, but shifts to favor lower luminosity and exhibits some bimodality as the probability density pinches at the median. This is due to the bimodality of the L₀ posterior distribution in Figure A1 given ${{ \mathcal D }}_{170817+R}$, which shall be discussed later in Section 5. The models with wide-angle structure tend toward lower mean intrinsic luminosity values, with the GJ model constrained to ${\mathrm{log}}_{10}{\hat{L}}_{0}/\mathrm{erg}\,{{\rm{s}}}^{-1}={50.4}_{-1.3}^{+1.6}$, PL model ${50.1}_{-1.6}^{+1.3}$, and DG model of ${50.2}_{-1.5}^{+1.5}$. The TH model is constrained to ${\mathrm{log}}_{10}{\hat{L}}_{0}/\mathrm{erg}\,{{\rm{s}}}^{-1}={50.2}_{-1.7}^{+1.7}$ which we can compare to the prior of ${51}_{-1.8}^{+1.8}$.

The constraints from the posteriors on the jet structure parameters given each jet structure model are shown in Table 5. The median is quoted along with the upper and lower bounds placed by the 90% narrowest credible intervals.

Table 5. Constraints Placed on Each of the Variables for Each Model Given Data Set ${{ \mathcal D }}_{\mathrm{all}}$ with and without the Fitted Luminosity Function

Model	Parameter	Constraints
		w/o Fitted LF	Fitted LF
TH	θj	${14.9}_{-14.3}^{{+46.0}^\circ}$	${9.2}_{-7.9}^{{+17.1}^\circ}$
GJ	θσ	${5.9}_{-5.4}^{{+28.2}^\circ}$	${4.2}_{-3.2}^{{+5.2}^\circ}$
PL	θc	${10.6}_{-10.0}^{{+32.9}^\circ}$	${6.0}_{-4.9}^{{+9.2}^\circ}$
	s	${6.7}_{-5.2}^{+8.1}$	${6.4}_{-4.3}^{+7.9}$
	a	${1.5}_{-1.4}^{+2.2}$	${1.5}_{-1.4}^{+2.1}$
DG	θin	${6.3}_{-6.2}^{{+24.1}^\circ}$	${3.6}_{-2.8}^{{+6.2}^\circ}$
	θout	${49.3}_{-38.0}^{{+40.7}^\circ}$	${46.2}_{-36.3}^{{+43.1}^\circ}$
	${\mathrm{log}}_{10}{ \mathcal C }$	$-{3.8}_{-2.2}^{+2.9}$	$-{4.1}_{-1.9}^{+1.9}$
	${\mathrm{log}}_{10}{ \mathcal A }$	$-{3.0}_{-2.9}^{+2.8}$	$-{3.1}_{-2.7}^{+2.9}$

Note. The median of each posterior distribution is quoted along with upper and lower bounds placed by 90% credible intervals.

Download table as:  ASCIITypeset image

4.1. Fitted Luminosity Function

The analysis is repeated but instead of assuming a prior distribution on the luminosity scale and shape, a luminosity function (LF) is fitted from the observed isotropic equivalent luminosity of sGRBs with associated redshifts. The values of L₀ and γ are taken from the mean fitted Schechter function in Mogushi et al. (2019) of ${\mathrm{log}}_{10}{L}_{0}/\mathrm{erg}\,{{\rm{s}}}^{-1}=51.6$ and γ = 0.55, fitted to the isotropic equivalent luminosity of 35 sGRBs.

The log evidence between each of the jet structure models and the different data sets are shown in Table A2, while the respective Bayes factors when given ${{ \mathcal D }}_{\mathrm{all}}$ between each of the jet structure models are shown on the right-hand entries of Table 4.

The posteriors on the local rate of BNS merger when given the fitted LF and ${{ \mathcal D }}_{\mathrm{all}}$ are shown in Figure 3 in the same format as Figure 2(a) where the widths of the violin plots indicate the probability density, the maximum and minimum sample is indicated by the extreme solid vertical lines, and median with the middle solid vertical line. The narrowest 90% credible intervals are indicated by the dashed vertical lines and are ${\mathrm{log}}_{10}{R}_{\mathrm{BNS}}/{\mathrm{Mpc}}^{-3}\,{\mathrm{yr}}^{-1}=-{7.0}_{-0.6}^{+0.7}$, $-{7.2}_{-0.8}^{+0.7}$, $-{6.8}_{-0.9}^{+0.8}$ for the GJ, PL, and DG models, respectively. For the TH jet model, the rate is constrained to ${\mathrm{log}}_{10}{R}_{\mathrm{BNS}}/{\mathrm{Mpc}}^{-3}\,{\mathrm{yr}}^{-1}=-{7.3}_{-0.9}^{+1.0}$.

Zoom In Zoom Out Reset image size

The constraints on the jet structure models given the fitted LF and ${{ \mathcal D }}_{\mathrm{all}}$ are presented in Table 5 for each jet structure model, with the upper and lower bounds representing the narrowest 90% credible intervals.

5.1. Wide-angle Jet Structuring

5.2. Cocoon Emission

5.3. Rate of Binary Neutron Star Mergers

5.4. Luminosity Function

5.5. GW190425

We provide an extensive Bayesian analysis that constrains the jet structure, intrinsic LF, and rate of BNS mergers as well as providing a comparison between competing jet structure models. This is achieved by combining four data avenues: (1) the parameter inference posteriors from a GW trigger, (2) the sGRB flux when a counterpart is detected or the detector flux upper limit otherwise, (3) the observation rate of detected sGRBs, and (4) the merger rate informed from GW observations. We perform this analysis using the GW triggers GW170817 and GW190425, GRB 170817A, the nondetection of a GW190425 counterpart, the rate of sGRB detections by the Swift detector within a 9.8 yr observation period, and a merger rate consistent with the constraints imposed by GWTC-2 (Abbott et al. 2020c). This provides us with the following results:

• 1.  
  The rate of BNS mergers is constrained within 1–1300 Gpc⁻³ yr⁻¹, improving upon the results of GWTC-3.
• 2.  
  Wide-angled jet structures prove more compatible with the given model than TH jet in explaining the observed number of sGRBs in the wake of the low observed isotropic luminosity of GRB 170817A.
• 3.  
  Slight evidence is provided for a cocoon-like wide-angle jet structure when considering the observed rate of sGRBs and GRB 170817A. However, the evidence becomes awash across all wide-angle jet structures when GW190425 is included in the analysis.
• 4.  
  While providing evidence against wide-angle structuring, the hypothesis that GW190425 had a typical sGRB counterpart with a GRB 170817A–like jet structure and would remain undetectable to the Fermi detector given full-sky coverage is feasible given the model.

The analysis was extended to consider a fitted intrinsic LF to further incorporate the detected flux and estimated redshifts of past sGRB detections. This provides the following results:

• 1.  
  The rate of BNS mergers is further constrained to 2–680 Gpc⁻³ yr⁻¹.
• 2.  
  Slight evidence for the GJ structure is provided, unless GW190425 is excluded in which the cocoon-like DG jet structure is equally favored.

However, we note that the fitting of the LF requires strong assumptions about the jet structure and therefore introduces a bias toward jet structures with wide central components. Interestingly, this bias does not appear to manifest in the resulting Bayes factors where the TH jet loses favor over the wide-angle jet structures. A future analysis will work to incorporate the flux measurements and inferred redshifts of detected sGRBs directly, and therefore avoid introducing this bias. Future work would also include incorporating afterglow data into the analysis for events that coincide with a GW detection (Lin et al. 2021), as well as additional very long baseline interferometry (VLBI) imaging data for GW170817 (Ghirlanda et al. 2019) or when VLBI data are available for future events, where the simple analytic method used to approximate the inclination angle from superluminal motion has been shown to be reasonably accurate (Fernández et al. 2022). Additionally, the short GRB population includes claims of "off-axis" events such as GRB 150101B (Troja et al. 2018b), an under-luminous GRB; here the inclusion of afterglow data in any future analysis could improve the marginal structured jet preference found via the current analysis which considers only the prompt emission and redshift.

We are grateful for computational resources provided by Cardiff University, and funded by an STFC (grant No. ST/I006285/1) supporting UK Involvement in the Operation of Advanced LIGO. The authors thank Shiho Kobayashi for fruitful discussions. F.H. was supported by Science and Technology Research Council (STFC) (grant No. ST/N504075/1). J.V., I.S.H. and M.J.W. are supported by STFC (grant No. ST/V005634/1). M.J.W. was also supported by STFC (grant No. 2285031). G.P.L. is supported by a Royal Society Dorothy Hodgkin Fellowship (grant Nos. DHF-R1-221175 and DHFERE-221005). E.T.L is supported by National Science and Technology Council (NSTC) of Taiwan (grant No. 111-2112-M-007-020).

http://github.com/Fergus-Hayes/Unpacking_the_jet

The log evidence for each data set given each model is provided in Table A1 while Table A2 provides the evidence for when the fitted luminosity function is assumed.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table A1. Log Evidence $\mathrm{ln}p({ \mathcal D }| M)$ Given the Five Different Data Sets When Assuming Each of the Four Models

Data ${ \mathcal D }$	Model M
	Top Hat	Gaussian	Power Law	Double Gaussian
${{ \mathcal D }}_{170817}$	11.87 ± 0.02	11.8 ± 0.03	11.79 ± 0.03	12.08 ± 0.03
${{ \mathcal D }}_{190425}$	17.29 ± 0.01	16.3 ± 0.03	16.41 ± 0.02	16.42 ± 0.03
${{ \mathcal D }}_{R}$	−7.67 ± 0.04	−7.8 ± 0.04	−7.73 ± 0.04	−7.53 ± 0.04
${{ \mathcal D }}_{170817+R}$	4.45 ± 0.05	5.23 ± 0.05	5.15 ± 0.05	5.37 ± 0.05
${{ \mathcal D }}_{\mathrm{all}}$	22.12 ± 0.05	22.66 ± 0.05	22.58 ± 0.05	22.67 ± 0.05

Download table as:  ASCIITypeset image

Table A2. Log Evidence $\mathrm{ln}p({ \mathcal D }| M)$ Given the Five Different Data Sets When Assuming Each of the Four Models Given the Fitted Schechter LF

Data ${ \mathcal D }$	Model M
	Top Hat	Gaussian	Power Law	Double Gaussian
${{ \mathcal D }}_{170817}$	11.48 ± 0.02	11.6 ± 0.03	11.45 ± 0.03	11.62 ± 0.03
${{ \mathcal D }}_{190425}$	17.28 ± 0.01	16.23 ± 0.03	16.29 ± 0.03	16.21 ± 0.03
${{ \mathcal D }}_{R}$	−8.08 ± 0.05	−8.53 ± 0.05	−8.23 ± 0.05	−7.84 ± 0.05
${{ \mathcal D }}_{170817+R}$	3.86 ± 0.05	4.97 ± 0.05	4.74 ± 0.05	4.91 ± 0.05
${{ \mathcal D }}_{\mathrm{all}}$	21.87 ± 0.05	22.78 ± 0.05	22.32 ± 0.05	22.51 ± 0.05

Download table as:  ASCIITypeset image

The full posteriors are shown in overlaid corner plots given all data sets in Figures A1, A2, A3, and A4 for the TH, GJ, PL and DG jet structure models, respectively.