The Jet Opening Angle and Event Rate Distributions of Short Gamma-ray Bursts from Late-time X-ray Afterglows

We present a comprehensive study of 29 short gamma-ray bursts (SGRBs) observed $\approx 0.8-60$ days post-burst using $Chandra$ and $XMM-Newton$. We provide the inferred distributions of SGRB jet opening angles and true event rates to compare against neutron star merger rates. We perform uniform analysis and modeling of their afterglows, obtaining 10 opening angle measurements and 19 lower limits. We report on two new opening angle measurements (SGRBs 050724A and 200411A) and eight updated values, obtaining a median value of $\langle \theta_{\rm j} \rangle \approx 6.1^{\circ}$ [-3.2$^{\circ}$,+9.3$^{\circ}$] (68\% confidence on the full distribution) from jet measurements alone. For the remaining events, we infer $\theta_{\rm j}\gtrsim 0.5-26^{\circ}$. We uncover a population of SGRBs with wider jets of $\theta_{\rm j} \gtrsim 10^{\circ}$ (including two measurements of $\theta_{\rm j} \gtrsim 15^{\circ}$), representing $\sim 28\%$ of our sample. Coupled with multi-wavelength afterglow information, we derive a total true energy of $\langle E_{\rm true, tot} \rangle \approx 10^{49}-10^{50}$\,erg which is consistent with MHD jet launching mechanisms. Furthermore, we determine a range for the beaming-corrected event rate of $\mathfrak{R}_{\rm true} \approx360-1800$ Gpc$^{-3}$ yr$^{-1}$, set by the inclusion of a population of wide jets on the low end, and the jet measurements alone on the high end. From a comparison with the latest merger rates, our results are consistent with the majority of SGRBs originating from binary neutron star mergers. However, our inferred rates are well above the latest neutron star-black hole merger rates, consistent with at most a small fraction of SGRBs originating from such mergers.

One of the most important parameters in determining the energy scales, rates and different progenitor channels of SGRBs is the opening angle of the relativistic outflow (e.g., Frail et al. 2001;Fong et al. 2015a;Mandhai et al. 2018), which can be inferred from the afterglow light curves.The opening angle is key in the determination of the true energy scale and rates of these events (e.g., Frail et al. 2001;Fong et al. 2015a;Mandhai et al. 2018), which are fundamental properties that can help differentiate between potential progenitor channels.Considering an observer aligned with the jet axis, the jet opening angles can be calculated from the detection of temporal steepenings, or 'jet breaks' in the broad-band afterglow light curves (Piran 1999;Rhoads 1999;Sari et al. 1999;Panaitescu 2005).
In the coming years, the continued upgrades of gravitational wave detectors (e.g., Advanced LIGO/Virgo/KAGRA Abbott et al. 2020a) will refine the localization of gravitational wave events and significantly increase the number of compact binary merger detections.One can use a combination of the SGRB event rate, the BNS/NS-BH merger rates constrained by the detection of gravitational waves, and the Galactic compact-binary merger rates, to place constraints on the fraction of BNS and NS-BH mergers that produce SGRBs.
Making use of all the available X-ray information collected by Chandra and XMM-Newton for Swift SGRBs detected since 2005, we perform a comprehensive study of bursts.Utilizing the broadband afterglow information of these events available in the literature, we constrain their energetics, circumburst densities and jet opening angles.We use the distribution of these opening angles to derive the true event rate of SGRBs and compare it with the observed and theoretically-predicted NS merger rates.In Section 2, we introduce the SGRB sample.In Section 3, we describe the uniform reduction and analysis of the X-ray data.In Section 4, we model the temporal afterglow behavior for each burst in our final sample of 29 SGRBs, and use the broad-band information of each event to constrain the energetics and circumburst densities.In Section 5, we present the distribution of SGRB jet opening angles, energy scales and event rates derived from our final sample.In Section 6, we discuss the implications of our findings on the energy scales and jet launching mechanisms, as well as compare our SGRB rate estimates with the BNS and NS-BH merger rates.In Section 7, we summarize the main results.In the Appendix, we present the models used for fitting the afterglow light curves, and further explain the statistical test used to determine the best-fit model.We also present a panel with the X-ray observations of SGRB 130603B that show the X-ray afterglow of this burst and the contribution of a contaminant source.
The cosmology employed in this paper is standard, with H 0 = 69.6 km s −1 Mpc −1 , Ω M = 0.286, Ω vac = 0.714 (Bennett et al. 2014).All the optical observations have been transformed to AB mag system and corrected for Galactic extinction in the direction of the burst (Schlafly & Finkbeiner 2011).Unless otherwise noted, uncertainties on median values correspond to 16th and 84th percentiles on the full distribution.

THE SHORT GAMMA-RAY BURST SAMPLE
In this paper, we present the analysis of the SGRBs that meet all the following criteria: 1. Bursts classified as short, or short with extended emission (SGRB-EE).Generally, we select bursts with Swift/BAT durations T 90 2 s (15-350 keV).We also include in the sample the following verified SGRB-EE events or events with potential EE (Lien et al. 2016): GRB 050724A with an initial hard, short pulse of ∼0.25 s and a soft, long component longer than 100 s in the light curve (Krimm et al. 2005;Barthelmy et al. 2005), and GRBs 080503 (Perley et al. 2009), 101219A (Krimm et al. 2010), 150424A (Barthelmy et al. 2015), 170728B (Ukwatta et al. 2017) and 200219A (Laha et al. 2020).In addition, we include GRB 050709, with T 90 ≈ 0.07 s (30 − 400 keV), that was discovered by the High Energy Transient Explorer (HETE ; Villasenor et al. 2005;Fox et al. 2005).
2. SGRBs with Chandra and/or XMM-Newton observations.In general, these observations are on the tail end or after Swift/XRT monitoring (referred as late-time observations from now on), starting at δt ≈0.8-1.7 days in the observer frame (where δt stands for the time between the burst trigger and the mid-time of the observation; Figure 1).
3. Events with well-sampled X-ray afterglows.To model the X-ray afterglow and place any constraint on the decline rate, we select those events with sufficient Swift/XRT information (more than two bins in the XRT light curve) at δt 1 day, that subsequently have been observed by Chandra and XMM-Newton.An exception for this rule is GRB 150101B, for which we ignore the XRT data since it is highly affected by the contribution of a neighboring Active Galactic Nucleus (Campana 2015;Fong et al. 2015bFong et al. , 2016)).Therefore, we only use the Chandra and XMM-Newton data sets for this burst.
4. Available broadband information of the SGRB afterglow.We require the afterglows to have been observed in the optical and radio bands.We use detections for both bands, when not available we utilize the upper-limit information, to determine the burst energetics and environment densities that allow to obtain constraints on jet opening angles.2021), Barthelmy et al. (2021) and Laskar et al. (2022), respectively.SGRBs with extended emission are indicated with (EE) (see Sec. 2).( 6) Redshifts (z), where "sp" denotes spectroscopic, "sp,AG" indicates spectroscopic afterglow and "ph" refers to photometric redshifts.( 7)-( 9) "Y" (yes) and "N" (no) indicate whether or not the burst position was observed by that observatory.( 10 Based on the initial SGRB sample published by Fong et al. (2015a), and the SGRBs or SGRB-EE events detected by either Swift/BAT1 or Fermi /GBM2 between 2015 − 2021, we find a total number of 118 SGRBs or SGRB-EE, out of which 90 were followed by the Swift/XRT and have clear X-ray afterglow detections3 .Applying our second criterion, the group of 90 SGRBs is reduced to 32 events observed by Chandra, XMM-Newton or both observatories at late times.We exclude GRBs 050813, 161104A and 210919A since their XRT light curves were only sampled by a single bin, and even in conjunction with the late-time observations, do not allow for reliable light curve modeling.All our criteria are ultimately met by 29 SGRBs that comprise our final sample.We collect a total of 60 late-time observations across all events in the sample: 18 bursts were observed by Chandra, 4 by XMM-Newton, and 7 additional bursts were monitored by both observatories.We present a uniform analysis of the late-time X-ray afterglow information of the 29 SGRBs with 10 bursts not being covered in the published, peer-reviewed literature: GRBs 150423A, 150424A, 150831A, 170728B, 180727A, 191031D, 200411A, 201006A, 210919A and 210726A.In Figure 1, we show that Chandra and XMM-Newton observations are critical in both the characterization of the X-ray behaviour of SGRBs and the identification of jet breaks in the light curves (see Section 5) at δt 0.8 days when Swift sensitivity is not enough to follow up the faint X-ray afterglows.
We list the basic properties such as X-ray positions, durations and redshifts of the 29 SGRBs in our sample and their available X-ray afterglow information from each X-ray observatory in Table 1.Our sample has redshifts of z ∼ 0.134 − 2.211, which is representative of nearly the entire range of SGRBs (e.g., Fong et al. 2022;Nugent et al. 2022;O'Connor et al. 2022).Eighteen of the events have confirmed spectroscopic redshifts from their host galaxies while six have photometric redshifts.Two bursts have constrained redshift intervals (GRBs 111020A and 211106A), one case (GRB 150423A) has a redshift determined from its afterglow and two events have an inconclusive host galaxy association (Fong et al. 2022).We include the redshift values in our subsequent analysis to put constraints on the explosion energetics and circumburst density properties (Table 4) of each burst in the sample, as well as deter-mining their half-opening angles (θ j ; referred as opening angle hereinafter, see Section 5).

X-RAY OBSERVATIONS & ANALYSIS
We perform a systematic reduction of all the available Chandra and XMM-Newton observations in the SGRB sample in order to model their X-ray spectra and obtain the unabsorbed X-ray flux afterglow light curves.

Chandra Observations
First, we collect the 24 Chandra observations out of 47 exposures that our team obtained through dedicated target of opportunity (ToO) and director discretionary time (DDT) programs (PIs: Berger, Fong, Rouco Escorial and Schroeder; Table 2) to follow up SGRB afterglows at late times.Additionally, we retrieved archival data of the other SGRBs (23 out of 47 exposures) in our sample utilizing the Chandra Search and Retrieval interface (ChaSeR4 ).All Chandra observations in this sample were obtained using the ACIS-S detector, except for GRB 110112A which was observed with the HRC.We used the CIAO 4.12 software package (Fruscione et al. 2006) with the CALDB 4.9.0 calibration database files for reducing the data.We reprocessed the data using the chandra repro task to obtain new Level II event files, and filtered each observation to omit time intervals with high background flares (using threshold rates of ≥ 0.5 − 0.6 counts s −1 in the 0.5 − 7 keV range for the ACIS-S detector).Then we ran the CIAO source detection tool, wavdetect, to perform a blind search for X-ray sources.We used the default wavelet scale values (2 and 4) to detect small features, and a threshold (sigthresh) of 10 −6 for identifying a source pixel.
Once we obtained a list of all detected X-ray sources for a given observation, we searched the Swift/XRT enhanced afterglow positions for a spatially coincident Chandra source.In case of a positive match, we calculated the 1σ uncertainty on the afterglow position as: where σ centroid is the afterglow centroid 1σ uncertainty provided by wavdetect (typically ≈ 0.02 − 0.3 ), and σ astro is the Chandra absolute astrometric uncertainty (1σ) of 0.65 .The positions and their uncertainties are listed in Table 1.We found that 14 SGRBs have at least one afterglow detection with Chandra6 .For these bursts, we obtained the count rates and spectral information from afterglows using circular source regions centered on the Chandra afterglow positions, with varied radii between 1.5 to 3.0 depending on the extension of the source.For the background count rates and spectra, we used source-free annuli centered on the source with default inner and outer radii of 30 and 60 , respectively.In some cases, we needed to vary the sizes of these regions slightly to avoid neighboring sources.
For each detection, we extracted the afterglow net count rate using the CIAO/dmextract tool in the Chandra 0.5 − 7 keV energy range.We also generated the source and background spectra, together with their ancillary response files (arf) and redistribution matrix files (rmf) using the specextract tool.In the case of non-detections, we converted the extracted Chandra count rate (0.5 − 7 keV) from 1.5 circular regions on the XRT positions into 3σ count rate upper limits using Poissonian confidence levels according to Table 1 in Gehrels (1986).
In addition, we detected neighboring X-ray sources close to the afterglow positions of GRBs 130603B, 180418A, 200522A and 210726A.These sources could not be resolved in the source regions of the Swift/XRT or XMM-Newton/EPIC observations 7 and fall within the SGRB point spread function (PSF).To correct the SGRB afterglow light curves from the contributing Xray flux of these contaminants, we extracted their spectra and modeled them to account for their spectral behavior in our analysis (see Section 3.4).

XMM-Newton Observations
Similarly to Section 3.1, we gather our 7 XMM-Newton observations of SGRBs from our ToO dedicated program (PI: Fong; Table 2).For the remaining 6 observations of SGRBs that our group did not obtain, we searched in the XMM-Newton Science Archive (XSA 8 ) for public observations.All EPIC observations were obtained utilizing the pn and both MOS detectors.We used the XMM-Newton Science Analysis System (SAS; version 18.0.0Gabriel et al. 2004) for reducing and analyzing the observations.First, we produced calibrated 7 Even though GRBs 180418A and 200522A have already published data sets corrected for these effects (Rouco Escorial et al. 2021;O'Connor et al. 2021, respectively), we re-perform this analysis for consistency. 8http://nxsa.esac.esa.int/nxsa-web/#searchevent lists by running the emproc and epproc tasks and filtered them for any background flaring activity.We use threshold rates of ≥ 0.4 − 0.5 counts s −1 in the 10−12 keV range for the pn detector, and threshold rates of ≥ 0.25 − 0.3 counts s −1 in the >10 keV range for the MOS detectors to filter out intervals of high background.
Next, we performed a blind search for X-ray sources running the edetectchain9 routine.This SAS task performs a simultaneous search for EPIC sources on background-filtered images extracted in five energy bands (0.2 − 0.5 keV for MOS and 0.3 − 0.5 keV for pn, and 0.5 − 1 keV, 1 − 2 keV, 2 − 4.5 keV, 4.5 − 12 keV for all detectors) and generates a final list of detected sources on a background-filtered image representing the full energy range.We then searched for detected sources in coincidence with the Swift/XRT positions; for relevant detected sources, we calculated the positional uncertainty (Equation 1), where σ 2 centroid is provided by edetectchain and σ astro = 1.5 , the XMM-Newton systematic error (Rosen et al. 2016;de la Calle Pérez et al. 2021).
We detected afterglows in at least one observation for 8 SGRBs.We extracted the count rates and spectra of the sources from 20 radius circular regions centered on the XMM-Newton positions10 in the 0.3 − 10 keV energy band.To obtain the background count rate and spectral information for each detection, we used a source-free region with a size identical to the source region, randomly placed on the same detector quadrant as the afterglow.We verified that none of the data were affected by pileup.We generated the rmf and arf files by running rmfgen and arfgen, respectively.When available, we jointly fit the EPIC/pn and EPIC/MOS spectra of each SGRB (see Section 3.4).If X-ray afterglows were not detected, we obtained 3σ upper limits on the count rates (0.3 − 10 keV) fixed at the positions of the Swift/XRT afterglows using the same size of the XMM-Newton detection regions utilizing the eupper SAS task.

Swift Observations
For Swift/XRT data, we collected all of the available X-ray afterglow information from the Swift light curve repository11 (Evans et al. 2007(Evans et al. , 2009)).We obtained the unabsorbed X-ray flux light curve by applying the relevant counts-to-unabsorbed flux conversion factors to the count-rate light curves, depending on the Swift observing mode (Windowed Timing (WT) and Photon Counting (PC) modes).The automatic spectral fitting routine fits each burst's X-ray spectrum to a double-component absorbed power law model with one of the absorption components set to the Galactic column density in the direction of the burst (N H,Gal ; Willingale et al. 2013), and the other accounting for any excess in the line of sight as an intrinsic neutral hydrogen absorption column (N H,int ), while the power-law component is characterized by the X-ray photon index (Γ X ).
For the four cases (GRBs 130603B, 180418A, 200522A and 210726A) in which the contaminating sources were discovered in Chandra observations (Section 3.1), we utilized the option 'create time-sliced spectra' in the Swift repository to obtain the spectra of each bin in the dynamically12 binned XRT light curves.Performing a spectral analysis in which we account for the contamination of these unrelated astrophysical sources leads to lower SGRB afterglow fluxes at the tail of the XRT monitoring.In Section 3.4, we describe the method to obtain the corrected unabsorbed X-ray flux light curves for each SGRB.

X-ray Spectral Analysis
We analyzed all of the X-ray afterglow spectra of the 29 SGRBs in our sample utilizing Xspec (12.10.1f;Arnaud 1996).Before proceeding with the spectral analysis, we binned the spectra with grppha to guarantee at least one count per bin to avoid negative counts when accounting for the background during the analysis.For modeling the X-ray afterglow spectra in Xspec, we used WILM abundances (Wilms et al. 2000) with X-ray cross-sections set to VERN (Verner et al. 1996) and statistics to W-statistics (Cash-statistics for Poisson data and background; Wachter et al. 1979).We characterized the spectra with a two-component absorption power-law model (tbabs*ztbabs*pow) defined by N H,Gal (Willingale et al. 2013), N H,int at the redshifts13 of the SGRBs when available (see Table 1) and Γ X .
Initially, we allowed all spectral parameters to vary as free except for N H,Gal and z.However, we find that for 14 bursts, the 3σ confidence interval of N H,int was consistent with zero.Thus, we performed a revised spectral fit for these events setting N H,int = 0 (see Table 2).For SGRBs initially monitored by XRT and with more than one detection with Chandra or XMM-Newton, we performed a joint spectral fit of the entire data set tying Γ X and N H,int between spectra and leaving the normalization as a free parameter, which is valid because we did not find evidence for spectral evolution for any of the bursts.We determined the best-fit spectral parameters within the 0.5−7 keV and 0.3−10 keV energy ranges for the Chandra and XMM-Newton spectra, respectively.
For SGRB afterglows that were detected by both XMM-Newton and Chandra, we set a common energy range of 0.5 − 7 keV, which is well-covered by both observatories, and used a constant multiplicative model (const*tbabs*ztbabs*pow) to account for the crosscalibration between the different instruments of the observatories in our joint fitting.We set the XMM-Newton/EPIC-PN constant value to 1 and calculated the rest of the constants (const EPIC−MOS1 = 1.088, const EPIC−MOS2 = 1.106, const ACIS−S3 = 1.106 and const XRT−PC = 0.965) following Table 5 in Plucinsky et al. (2017).
Finally, we determined the unabsorbed X-ray fluxes (F X ) of the afterglows within the 0.3 − 10 keV energy range, fixing the spectral parameters to the bestfit values and using the cflux convolution model (const*tbabs*ztbabs*cflux*pow) in Xspec.The unabsorbed X-ray fluxes (0.3−10 keV), together with bestfit parameters and their 1σ uncertainties are listed in Table 2.We calculated the 3σ X-ray unabsorbed flux upper limits from the upper-limits of the count rates, using the Heasarc WebPIMMS tool14 with the best-fit spectral parameters available from previous afterglow detections for each case.The 3σ upper limits for the X-ray unabsorbed fluxes are listed in Table 2.
There are four cases with contaminating sources (X1) to the Swift/XRT or XMM-Newton source regions that are resolved in Chandra observations.We employed a single absorbed power-law model (tbabs*pow) to characterize their Chandra spectra.Then we took into account these spectral parameters in the Swift and XMM-Newton spectral modeling of the SGRB afterglows and calculated the unabsorbed X-ray flux from the afterglow in the 0.3 − 10 keV energy band (following the same method as described in Rouco Escorial et al. 2021).For that, we set cflux to account only for the afterglow (AG) spectral component of the model: (tbabs*zbabs*cflux*pow) AG + (tbabs*pow) X1 .(5) Elapsed time between the burst trigger and the log-centered time of the observation (δt).( 6) Final exposure time after filtering from background events.( 7)-( 8) Galactic column density, N H,Gal , and intrinsic absorption value, NH,int.The symbol "-" means that the contribution of NH,int is negligible and is set to 0 in our fits.(9) X-ray photon index (ΓX).The spectral parameter ΓX is calculated in the 0.5 − 7 keV energy range for Chandra and between 0.3 − 10 keV for XMM-Newton (when information from both observatories is available then a 0.5 − 7 keV energy range is applied).For observations with joint spectral analysis, only the first value of ΓX is shown with uncertainties, while for the rest of observations just the value of ΓX is stated.( 10

INFERRED PROPERTIES FROM BROADBAND MODELING
Here we model the broad-band afterglows of our sample considering the standard synchrotron forward shock model (Sari et al. 1998;Granot & Sari 2002).In this scenario, the forward shock originates from the interaction of a relativistic blast wave with a constant density medium, which is the expected environment for a NS merger (e.g., Paczynski 1986;Eichler et al. 1989).
In what follows, we first determine the temporal and spectral evolution, parameterized by F ν ∝ t α ν β , where α and β are the temporal and spectral power-law indices.We then use these indices and the closure relations given by the synchrotron model (Granot & Sari 2002) to determine the power-law index of the input distribution of accelerated electrons (p), providing a framework for determining the locations of the observing bands with respect to the three break frequencies of the synchrotron spectrum: the self-absorption frequency (ν sa ), peak frequency (ν m ) and cooling frequency (ν c ).We then fit the available observational data using the mapping from Granot & Sari (2002) to determine the burst physical properties.These properties include the isotropic-equivalent energy of the jet (E K,iso ), the density of the circumburst environment (n 0 ), and the fractions of the post-shock energy transmitted to electrons ( e ) and magnetic field ( B ).In Section 5, we apply these properties to determine the jet opening angles of the bursts.

X-ray Temporal and Spectral Fitting
The temporal behavior of X-ray afterglows can be well-described with a power-law or a series of connected power laws (e.g., Zhang et al. 2006;Evans et al. 2007;Margutti et al. 2013).To ensure that we are only using data during the forward shock afterglow phase, we ignore all of the data at δt < 200 s, which often includes high-latitude emission and the tail of the γ-ray emission (Zhang et al. 2006;Racusin et al. 2009, and reference therein).Additionally, during the fits we omit portions of the light curves that are clear deviations from the main power laws15 .
To remain agnostic to the particular best-fit model, we first fit every SGRB afterglow light curve using single, broken and triple power-law models, described by Equations A1, A2 and A3 (see Appendix A for more de-tails).The main parameters of our models are the series of temporal decay indices (α 1,2,3 ), the flux normalization (C), the break times (t b1,2 ), and the smoothness parameters (s 1,2 ).For every fit, we set all model parameters as free, except for the smoothness values which are set as constants in the broken and triple power-law models.
To perform the light curve fitting, we use the Python emcee package (Foreman-Mackey et al. 2013) in which we apply an affine-invariant Markov chain Monte Carlo (MCMC) ensemble sampler (Goodman & Weare 2010;Foreman-Mackey et al. 2019) to find the solution with the highest log-likelihood for each model together with the model parameters.We initiate each fit with 15000 steps, each with 100 walkers, and discard the initial 500 steps as "burn-in" when the average log-likelihood across the chains reaches a stable value.We then calculate the medians and 1σ credible regions from the posteriors for each parameter, and use these values to determine the reduced-χ 2 value for each model.
Based on the results from our light curve fitting procedure for each burst, we use the F-test to determine whether a simpler model provides a statistically better fit to our data than a more complex model.In effect, we first test the null hypothesis (H 0 ) that a single power-law model is sufficient to describe the X-ray afterglow light curve.If H 0 is not accepted, then we test a second null hypothesis (H' 0 ) that a broken power law is the best-fit model.If H' 0 is not accepted, the afterglow light curve is best-described then by a triple power law model (see Appendix B).We find that twelve SGRBs are best-fit by a single power-law model, eleven by a broken powerlaw, and six by a triple power-law model.The medi-ans and 1σ uncertainties of the best-fit parameters for each SGRB are shown in Table 3.For afterglows best described by a single power-law, we find a median of α X ≈ −1.00, which is in agreement with the weighted mean of α X found by Fong et al. (2015a).
For SGRBs with light curves best described by the broken power-law model, we find two scenarios: (i) α 1 < α 2 in which the afterglow light curve transitions from a steep power-law of α 1 −1 to a shallower power-law of α 2 ≈ [−1.5, −0.3], or (ii) α 1 > α 2 in which the afterglow light curve steepens.Cases in which α 2 −2 can be associated, in general, with jet breaks (Section 5).
For the triple power-law model, we identify a range of morphologies in terms of afterglow behavior.Similar to the broken power-law afterglows, bursts with their final segments following the condition, α 3 −2 can be viable jet breaks (see Section 5 for further criteria).We note that our analysis also uncovers several instances of shallow power-law segments of α X ≈ −0.6, which have been proposed as phases related to an extra injection of energy, for example driven by a potential magnetar central engine (Rowlinson et al. 2013;Gompertz et al. 2014;Lü et al. 2017;Stratta et al. 2018).
We use where k bol is the bolometric correction factor to convert the Swift/BAT fluence (fγ ) in the 15 − 150 keV energy band to the 1 − 10000 keV energy band (k bol = 5), and d L is the luminosity distance (for SGRBs with unknown redshifts in Table 1 we assign the fiducial value of z = 0.5).The isotropic-kinetic energy and circumburst density are calculated with the fraction of the post-shock energy transmitted to electrons ( e) fixed to 0.1 for the different values of B and redshifts (Table 1).We define The E γ,iso value of GRBs 050709 is given in the 1 − 10000 keV energy band by Villasenor et al. (2005).† Values obtained from the analysis of Laskar et al. (2022) considering z = 1 and IC/KN corrections.

Isotropic-equivalent Kinetic Energies and Circumburst Densities
Here we constrain the values of p, E K,iso and n 0 as described in Section 4 using any available X-ray, optical, NIR and radio data in the framework of the standard synchrotron forward shock model (Sari et al. 1998;Granot & Sari 2002); these are key parameters which feed into the determination of the opening angles.For the SGRBs detected during 2005-2015, we utilize the broadband analysis from Fong et al. (2015a) for the values of p, E K,iso and n 0 .We also use the analysis of Schroeder et al. (2020) to determine these parameters for the particular cases of GRBs 050509B, 150424A and 160821B, and the latest analysis of GRB 211106A from Laskar et al. (2022).For the remaining SGRBs discovered beyond 2015, we follow the methods of Fong et al. (2015a) to constrain these shock and environment parameters.
First, we constrain the position of ν X with respect to ν c (e.g., ν X > ν c or ν m < ν X < ν c ), assuming a constant density medium.We determine under which scenario the values of p, independently determined from α X and β X , adhere to the closure relations given by Granot & Sari (2002).For each SGRB, we require that the values of p derived from α X and β X to be consistent within 1σ, calculate the weighted mean (Table 4) and assume this value as p.For two bursts, GRBs 180727A and 210726A, we obtain p = 1.92 ± 0.20 and p = 1.97 ± 0.03, respectively, which yield divergent total integrated electron energies.Thus we use p derived from β X for these cases as this parameter is generally more reliable since it remains constant over time, which both give p > 2. Additionally, in the case of GRB 200522A we find a weighted mean of p = 1.78 ± 0.30; therefore we adopt p = 2.05 ± 0.03 reported in Fong et al. (2021) for our analysis.Deriving from the standard closure relations, we find ν X > ν c for 10 cases and ν X < ν c for 19 cases, with an inferred weighted mean of p = 2.30±0.32 for the total sample of 29 SGRBs.This value is in agreement with the p = 2.40 +0.28  −0.36 found by Fong et al. (2015a).
Next, we use the broad-band data to calculate the physical burst properties, E K,iso and n 0 , following the relations of the synchrotron model (Granot & Sari 2002), which provides a mapping from the light curve and SEDs to physical properties.We assume that the radio band follows ν sa < ν rad < ν m , and that the optical band is located at ν m < ν opt < ν c .We utilize the location of ν c as an additional constraint.If ν c < ν X , we set the maximum value to ν c,max = 7.3 × 10 16 Hz, at the lower edge of the X-ray band (0.3 keV).If ν m < ν X < ν c , we constrain the minimum value of ν c,min = 2.4 × 10 18 Hz corresponding to the upper edge of the X-ray band (10 keV).
In general, we first fix e and B to the equipartition values of e = 0.1 and B = 0.117 (Table 4).This choice of microphysical parameters is based on results from GRB afterglow fitting (Ryan et al. 2015;Beniamini & van der Horst 2017) and theoretical studies which show that part (typically ∼ 10%) of the kinetic energy may be deposited into non-thermal particles (Spitkovsky 2008).For eight cases, we need lower values of B = 10 −4 − 0.01 in order to find E K,iso and n 0 solution pairs (e.g., Panaitescu & Kumar 2002;Sironi & Spitkovsky 2011).In particular, these cases are 150424A, 160821B (Schroeder et al. 2020), 191031D, 200219A and200522A (Fong et al. 2021) with B = 0.01, GRB 140903A with B = 10 −3 (Fong et al. 2015a), GRBs 050724A with B = 10 −4 (Fong et al. 2015a) and 211106A with B = 10 −5 (assuming z = 1 and IC/KN corrections Laskar et al. 2022).Additionally, we set z = 0.5 for SGRBs with no redshift information (Table 1).We point out that the derived physical parameters will be degenerate with respect to the fraction of particles accelerated into a non-thermal distribution (Eichler & Waxman 2005).
To calculate the distributions of E K,iso and n 0 , we consider a parameter space with 1000 logarithmicallyspaced steps for each parameter, with ranges of E K,iso = 10 46 − 10 56 erg and n 0 = 10 −6 − 10 3 cm −3 , where n min = 10 −6 cm −3 corresponds to the typical value of the intergalactic medium.We then constrain the E K,iso − n parameter space by combining the probability distributions from each wavelength for each SGRB, marginalize in each parameter, and normalize the 1D distributions.

JET OPENING ANGLES, BEAMING-CORRECTED ENERGETICS AND EVENT RATES
As a consequence of the relativistic beaming, the observed temporal behavior from a spherical expansion for an on-axis observer is initially similar to that of a highly collimated relativistic outflow.As the jet interacts with the circumburst medium, the value of its bulk Lorentz factor (Γ) declines over time to reach a value equal to the inverse of the jet opening angle (Γ ∼ θ −1 j ; Piran 1995) and consequently the edge of the highly collimated relativistic outflow becomes visible to the observer.This moment is known as the time of the jet break (t j ) and observationally manifests as an achromatic and temporal steepening in the afterglow light curve (Sari et al. 1999;van Eerten & MacFadyen 2013), after which the outflow may follow lateral expansion (Granot & Piran 2012).At this point, the jet undergoes a lateral expansion and the afterglow decline rate can be analytically described by t −p (Rhoads 1999;Sari et al. 1999).By constraining the parameter t j in the SGRB afterglow light curves, one can additionally derive θ j .

Determination of Jet Opening Angles
To identify jet breaks from the best-fit models, we define a jet break as the time when α significantly steepens.Quantitatively, for α X,i > α X,i+1 , we classify light curves as having jet breaks if they exhibit: (i) α X,i+1 < −1.75, and (ii) 0.75 < |∆α X | < 3.For instance, considering an ISM-like medium for the closure relations given by Granot & Sari (2002) and assuming 2 < p < 3 (Fong et al. 2015a), we know that α X can range between −1.5 < α X < −0.75 if ν X < ν c and −1.75 < α X < −1 if ν c < ν X , in case of a spherical outflow.Therefore, if the post-break segment of a light curve is steeper than α X = −1.75, the only real justification is a jet break instead of the transition of ν c across the X-ray band.Additionally, the transition of ν c across the Xray band predicts a maximum value of the change in slope of |∆α X | = 0.25 (Sari et al. 1998), so that any larger ∆α is most naturally explained by a jet break.
Considering an ISM-like medium and a top-hat geometry for the cosmological SGRB jets with the observer along the jet axis 18 , we combine the information from the X-ray afterglow light curves with constraints on E K,iso and n 0 (Section 4) to calculate θ j or their lower 18 One of the lessons learnt from GRB 170817A is that jets can be more complex and structured.Wu & MacFadyen (2019) showed that cosmological SGRBs share a similar jet structure to GRB 170817A.Likewise, Urrutia et al. (2021) demonstrated that structured jets can have significantly more injected energy at large observing angles than top-hat jets.However, for our SGRB sample, the emission is completely dominated by the jet core making the detection of any structure difficult.Therefore, the assumption of a SGRB top-hat jet geometry is a reasonable approximation for our study.
limits.We use the following equation given by Sari et al. (1999);Frail et al. (2001): where t j is time in days, E K,iso,52 is in units of 10 52 erg and n 0 is in units of cm −3 .For bursts with detected jet breaks we set t j to the best-fit time of the break in our calculation of θ j (Table 3).In the case of SGRBs with no detected jet breaks (e.g., well-described by single powerlaw declines), we fix t j to the time of the last X-ray detection (considering this time as a lower limit on t j ) to obtain lower limits on the jet opening angles (Figure 3).

Distribution of SGRB Jet Opening Measurements
Using the aforementioned criteria, we find that nine SGRBs (SGRBs 050724A, 051221A, 111020A, 130603B, 140903A, 150424A, 160821B, 200411A and 200522A; Table 3) exhibit temporal steepenings in their X-ray afterglow light curves that we attribute to jet breaks.We also add to this group the case of GRB 211106A, for which the jet break was uncovered in the radio afterglow light curve (Laskar et al. 2022).Thus, our final subsample is composed of 10 SGRBs with jet breaks.We confirme the existence of the jet breaks for eight cases that have been already published: GRB 051221A (Burrows et al. 2006;Soderberg et al. 2006), GRB 111020A (Fong et al. 2012), GRB 130603B (Fong et al. 2014), 140903A (Troja et al. 2016), GRB 150424A (Jin et al. 2018), GRB 160821B (Lamb et al. 2019;Troja et al. 2019), GRB 200522A (O'Connor et al. 2021) and GRB 211106A (Laskar et al. 2022).However, we update the time of the jet breaks on the X-ray afterglow light curves for all these cases based on our uniform modeling (see Table 4), except the jet break time of GRB 211106A that was detected in the radio band.Additionally, we report on two new jet break detections in X-rays that have not been found to date: GRBs 050724A and 200411A (see Table 4).For these ten SGRBs, we find that their median jet break times span a range of t j ≈ 0.1−30 days.We use the posteriors of t j for each case to build the total probability distribution of jet break times (Figure 3), and derive a median value of t j = 1.5[−1.2,+7.7] days.
To build the total distributions of jet opening angle measurements, we perform 5000 random draws from the probability distributions in E K,iso , n 0 and t j posterior for each jet break detection case.Since the distributions of E K,iso and n 0 are correlated, we select the appropriate value of n 0 for every drawn value of E K,iso .We then use Equation 2 to obtain individual posterior distributions of opening angles for the ten cases with detected jet breaks (Figure 4, left).The median and 1σ confidence intervals for each jet measurement are reported in Table 4.  .7]days.The distribution at tj > 20 days is clearly dominated by a peak which represents the late-time detection of the GRB 211106A jet break in the radio band.For the 19 lower limits, the times of the last X-ray detections are indicated with color-coded lower limits which roughly corresponds to the time out to which we are sensitive to jet breaks for each event.A significant fraction (10/19) of the lower limits are found at times of tj > tj , and specifically there are 6 events with last X-ray detections at tj > 10 days.This shows the existence of a SGRB group with jet breaks located at very late times (and potential wide opening angles) as the case of GRB 211106A.
We find that the angles for GRBs 051221A, 111020A, 140903A, 160821B and 200522A are in agreement within the errors with those reported by Burrows et al. 2006;Soderberg et al. 2006;Fong et al. 2012;Troja et al. 2019;O'Connor et al. 2021, respectively.However, we find an opening angle of ∼ 3 • for GRB 140903A narrower than θ j ≈ 5 • (Troja et al. 2016), and a narrower opening angle measurement of ≈ 4 • for GRB 150424A instead of θ j ∼ 7 • (Jin et al. 2018).For GRB 050724A, we now measure an opening angle of θ j ≈ 34 • , which is consistent with the previous wide opening angle lower limit of θ j 25 • (Grupe et al. 2009).The detection of a new jet break here is driven by the last Chandra detection at δt ≈ 21.6 days for which we find a ∼ 0.8 times fainter flux (see Table 2) than that reported by Grupe et al. 2006.This measurement is a limiting case of jet opening angles since theoretical studies of postmerger outflows derive a maximum value of θ j,max ≈ 30 • (Ruffert & Janka 1999;Aloy et al. 2005;Rosswog 2005;Rezzolla et al. 2011;Lazzati et al. 2021).
In the case of GRB 130603B, Fong et al. (2014) found a jet break at δt ≈ 0.47 days in the optical and radio afterglow light curves, but not in the X-ray band.However, our new MCMC modeling shows a jet break at ≈ 0.11 days post-trigger in the X-ray light curve.This earlier jet-break time is driven by the correction of the contribution of the X-ray contaminant (Figure 8 in Appendix C) to the afterglow light curve.Therefore, one does not need to invoke extra energy injection from a magnetar to explain the X-ray flux level of late-time observations.We constrain the opening angle measurement of GRB 130603B to θ j ≈ 4 • , which is in agreement with the lower end of the opening angle range (≈ 4 • −8 • ) for this event reported by Fong et al. (2014).
For the jet opening angle measurements.Overall, we find a median value of θ j,det = 6.1 for bursts with jet opening angle measurements, which is consistent with previous studies (Fong et al. 2015a;Jin et al. 2018) and in line with the opening angle estimates of jets for SGRBs produced by BNS mergers (Ghirlanda et al. 2016).In the case of SGRB 170817A, the first bona-fide detection of an SGRB jet launched by a BNS merger, the range of values inferred for the core is which is narrower than the median value we find for the SGRB jet opening angle measurements.It is clear from Figure 4 (left) that the probability beyond θ j 15 • is dominated by two events, GRBs 050724A and 211106A.

Distribution of Jet Lower Limits: A Population of Wider Jets
We also build probability distributions for the 19 SGRBs with inferred lower limits on the opening angles.We again perform 5000 random draws from the correlated distributions of E K,iso and n 0 .For these cases, we instead fix the jet break time to the (log-centered) time of the last X-ray detection.Using Equation 2, we obtain the posterior distribution of the opening angle lower limit for each event, and report the median value of the minimum opening angle and 1σ confidence intervals (Table 4).We find that the range of minimum opening angles varies between θ j,min ≈ 0.3 • -26 • .
We find that the minimum opening angles of 6 events (SGRBs 050709, 100117A, 111117A, 120804A, 150101B, and 180418A) are 1σ larger than θ j,det , indicative of a population of wider jets.To investigate if the distribution of lower limits is statistically distinct from the distribution of jet measurements, we compare the CDFs at their 68% credible regions.If the median of each distribution lies within the other's 68% credible region, then the incorporation of the lower limit distribution does not add new information.Indeed, we find that the median values are within the 16th and 84th percentiles of each other's distribution.Thus, at first glance the distribution of jet measurements appear to be a fair approximation for the observed opening angles of SGRBs.
However, we acknowledge that this may be the consequence of bias due to the limited number of events observed at t j 2 − 3 days, for which we expect to have wider jet opening angles.Indeed, we have already detected two events with wide opening angle measurements and found six cases with minimum opening angles wider than θ j,min 10 • .If more SGRBs were able to be monitored at late times (which would require more sensitive X-ray facilities than currently available), this could result in measurements or lower limits of broader jet opening angles.Therefore, for our subsequent analysis, we build a mock sample including a fiducial population of wider jets.In particular, we combine the events with opening angle measurements and six wide jets (representing the number of events with wider opening angle lower limits) with measurements centered around θ j ∼ 20 • (motivated by the findings in Murguia-Berthier et al. 2017) to investigate how an SGRB sample with constraining wider opening angles affects our results.

Energy Scale and Event Rate
The jet opening angle distribution affects the calculation of the burst true energy scales as well as the true event rate.The beaming-corrected energy is given by, E K = f b E K,iso , where f b ≡ 1 − cos(θ j ) is the beaming correction factor and E K,iso is the isotropic-equivalent kinetic energy.At the same time, emission from the relativistic jet is only immediately visible if the observer's line-of-sight is inside or intercepts the cone of the outflow.Thus, we expect the true event rate to be larger by a factor of f −1 b , where R true = f −1 b R obs .For the SGRBs with opening angle measurements in our sample, we build the CDF of f b by calculating the beaming correction factors using every value that composes the distributions of opening angles (Section 5.2).We find a median of f b,det = 0.6[−0.5,+3.0] × 10 −2 for the 10 bursts with well-measured opening angles.This value is consistent within errors to that quoted in Fong et al. (2015a).The f b value is essentially a minimum value on the median beaming correction factor as the sample used to derive it comprises jet measurements which likely represent the narrowest jets.For the mock sample including wide jets, the median value expectedly shifts to a larger value of f b = 2.8[−2.6,+4.1] × 10 −2 .
Using bursts with opening angle measurements, we find population median isotropic-equivalent kinetic and γ-ray energies of E K,iso ≈ 3.7 × 10 51 erg and E γ,iso ≈ 9.7 × 10 50 erg (1 − 10000 keV), respectively, with a γray efficiency of η γ ≈ 0.2.Incorporating f b,det , we obtain median beaming-corrected kinetic and γ-ray energies of E K ≈ 2.1×10 49 erg and E γ ≈ 5.0×10 48 erg, respectively.This results in a total beaming-corrected energy of E true,tot ≡ E K + E γ ≈ 2.6 × 10 49 erg.If we assume similar isotropic-equivalent kinetic and γ-ray energies for the mock sample including wide jets, the median total beaming-corrected energy of the SGRB population increases to E true,tot ≈ 1.3 × 10 50 erg.These values are similar to the energy scales found in previous works where E K ∼ 10 49 erg and E true ∼ 10 50 erg (Berger 2014;Fong et al. 2015a;Jin et al. 2018), respectively.In addition, even though the jet core of GRB 170817A is more collimated than most of our SGRB measurements, its total inferred energetics are similar with E K ∼ 10 49 − 10 51 erg (Margutti & Chornock 2021).
In the case of event rates, we correct the observed local rate (R obs ) derived from the SGRB luminosity function using the distribution of beaming factors.We note that several works which derive R obs via the minimum γray luminosity (e.g., Guetta & Piran 2006;Nakar et al. 2006;Coward et al. 2012;Wanderman & Piran 2015;Ghirlanda et al. 2019;Liu & Yu 2019), have converged on R obs ≈ 5 − 10 Gpc −3 yr −1 .At the lowest end,   (Margutti & Chornock 2021), is indicated as a green bar.right: The total probability distribution for the sample of 19 lower limits (red dashed histogram) in comparison with the total distribution of detected opening angles (dark grey histogram).The red arrow indicates that this distribution represents lower limit values.There is a distinctive group of SGRB lower limits composed of 6 events with wider minimum opening angles of θj,min 10 • .Thus, this indicates a possible population of SGRBs with wider opening angles than described from the distribution of measurements alone.Ghirlanda et al. (2019) derived R obs ≈ 0.5 Gpc −3 yr −1 based on available observational constraints of the Fermi SGRB population.However, this lower value is inferred from bright events with L min,iso ≈ 10 51 erg s −1 , ignoring a fraction of SGRBs that exhibit fainter luminosities.To account for the contribution of these fainter events with L min,iso ≈ 10 49 erg s −1 (e.g., Guetta & Stella 2009;Wanderman & Piran 2015), we consider a local rate of R obs ≈ 10 Gpc −3 yr −1 (Nakar et al. 2006) going forward.
Based on events with well-measured opening angles and R obs ≈ 10 Gpc −3 yr −1 , we find a median value of R true = 1789[−1510, +6334] Gpc −3 yr −1 (Figure 5).This derived true event rate is consistent within errors to previously published values based on SGRBs, albeit is on the high end (Coward et al. 2012;Fong et al. 2014Fong et al. , 2015a;;Jin et al. 2018;Dichiara et al. 2020) (Figure 5).This can be naturally explained because this rate comprises jet measurements which likely represent the narrower end of the population and ignores the existence of wider jets and lower limits.Thus, R true ≈ 1790 Gpc −3 yr −1 can actually be considered a high estimate on the true event rate.If instead we use the SGRB mock sample which includes wider jets and is likely a better representation of the parent distribution, we find a lower value of R true,mock = 362[−217, +4376] Gpc −3 yr −1 (Figure 5).Therefore, one would expect the real SGRB event rate to lie between those derived from both the events with jet opening angle measurements and mock distribution.In this section, we discuss the implications of our findings in the context of the event energy scales and potential mechanisms to launch jets.We also compare the SGRB true event rate derived from our study with other observational and theoretical rates published in the literature for these events and BNS/NS-BH mergers.Additionally, we investigate what fraction of the SGRB population can be explained by these mergers.
6.1.Jet Opening Angles: Implications on Energy scales and Potential Launch Mechanism The jet opening angles uniquely enable a determination of the true, beaming-corrected energy scale of SGRBs, which can be used to probe the energy extraction mechanism to power the jets.The calculation of the true energy for these events is key to discerning between the potential mechanisms to launch relativistic outflows (Shibata & Hotokezaka 2019), either by neutrino pair (ν ν) annihilation (Jaroszynski 1993;Mochkovitch et al. 1993;Rosswog & Ramirez-Ruiz 2002) or as a magnetically driven jet (Blandford & Znajek 1977;Rosswog et al. 2003;Rezzolla et al. 2011;Siegel & Metzger 2017) since one expects different released energy ranges.From our analysis, we find median beaming-corrected total energy releases of E true,tot ≈ (0.3 − 1) × 10 50 erg.
In the ν ν annihilation scenario, launched jets have expected opening angles of θ j ∼ 5 − 30 • and maximum beaming-corrected energies of ∼ 10 48−49 erg (Rosswog & Ramirez-Ruiz 2002;Aloy et al. 2005;Birkl et al. 2007;Dessart et al. 2009;Murguia-Berthier et al. 2017).Liu et al. (2015) and Perego et al. (2017) demonstrate that the deposited energy by the ν ν mechanism in comparison to the median inferred energy of ≈ 10 50 erg for detected SGRBs is not sufficient to power jets in these events.Even under the assumption of smaller opening angles of θ j < 10 • , it is still not sufficient to explain SGRB jets powered by the energy extracted from the ν ν annihilation mechanism (Perego et al. 2017).Based on this, the ν ν annihilation mechanism is unlikely to be the dominant energy extraction mechanism to launch SGRBs.
On the other hand, magnetohydrodynamic (MHD) processes (i.e., Blandford-Znajek mechanism;Blandford & Znajek 1977) can easily reach larger energy scales of ∼ 10 49−52 erg (Rosswog 2005;Lee & Ramirez-Ruiz 2007;Ruiz et al. 2016;Siegel & Metzger 2017).We note that there are different predictions for the jet opening angles depending on the outflow's magnetization (Rosswog & Ramirez-Ruiz 2002;Duffell et al. 2018;Nathanail et al. 2020) with expected θ j 10 • for more magnetized jets (Nathanail et al. 2020).In addition, Christie et al. (2019) demonstrated with 3D MHD simulations that poloidal post-merger magnetic fields generate jets with θ j ∼ 6 • − 13 • and up to E K,iso ∼ 10 52 erg, while toroidal post-merger magnetic field geometry produces jets with E k,iso ∼ 10 51 erg and θ j ∼ 3 • − 5 • .Indeed, this last magnetic field configuration of the post-merger disk is consistent with the jet opening angle and energetics found for GRB 170817A (Margutti & Chornock 2021), as well as the rather low efficiency of this event (η ∼ 10 −3 ; Salafia & Giacomazzo 2021).Polarimetry studies specifically focused on the radio reverse shock of long GRBs have been essential to reveal their magnetic field configurations and hence the jet launching mechanism (e.g.; Granot & Taylor 2005;Laskar et al. 2019).However, these type of studies are particularly difficult for SGRBs as they are fainter and evolve much faster than long GRBs so in the meantime, energetics can offer strong clues.Indeed, the agreement between the jet energetics and opening angles derived from the MHD scenarios and those inferred from our observations supports this scenario as the main mechanism to launch SGRB jets.

Derived Event Rates and Implications for Compact Object Merger Progenitors
One of the most important consequences of determining jet geometries is the inference on the true event rate (R true ).In particular, the estimate of R true for SGRBs can be compared with the rates of BNS (R BNS ) or NS-BH (R NS−BH ) mergers derived from the Advanced LIGO/Virgo gravitational wave detections (Abbott et al. 2017a(Abbott et al. ,b,c, 2020b(Abbott et al. , 2021)).This provides us with basic information on the fractions of BNS and NS-BH mergers that may power SGRBs and launch jets.Coupled with ejecta masses inferred from kilonovae (Gompertz et al. 2018;Rossi et al. 2020;Rastinejad et al. 2021), the true event rate can also determine the role of BNS/NS-BH mergers in populating the universe with r-process elements (e.g., Kasen et al. 2017;Hotokezaka et al. 2018;Rosswog et al. 2018 To date, several works have derived an observed event rate of R obs ≈ 5 − 10 Gpc −3 yr −1 (e.g., Guetta & Piran 2006;Nakar et al. 2006;Coward et al. 2012;Wanderman & Piran 2015;Ghirlanda et al. 2019;Liu & Yu 2019) using the minimum γ-ray luminosity of SGRBs.However, the lower end of R obs ≈ 5 Gpc −3 yr −1 is derived from the most luminous events.Therefore, for calculating the SGRB true event rate from the jet opening angle measurements and mock distributions, we choose the value of R obs ≈ 10 Gpc −3 yr −1 (assuming L iso,min ≈ 10 49 erg s −1 ; Nakar et al. 2006), which has been used in previous literature (Metzger & Berger 2012;Fong et al. 2015a;Mandel & Broekgaarden 2022).In Figure 6, we present our distributions of R true and compare these rates with BNS and NS-BH merger rates published in the literature as well as previous SGRB studies.We use We present a more realistic range of values of Rtrue that extends between the median values of the true event rate distributions derived from the opening angle measurements (vertical black dashed line) and mock samples (vertical grey dashed line).Other published rates are also shown as comparison.In particular, these rates are derived from: the detection of gravitational waves generated by BNS (Abbott et al. 2020a) and NS-BH mergers (Abbott et al. 2021), SGRBs (Fong et al. 2015a;Jin et al. 2018;Dichiara et al. 2020), Galactic BNS (O'Shaughnessy & Kim 2010;Pol et al. 2020;Grunthal et al. 2021), and the estimates for the BNS and NS-BH merger rates derived from population synthesis simulations by Sarin et al. (2022).68% confidence levels are represented.the range of rates, R true ≈ 360−1790 Gpc −3 yr −1 , where the lower end is set by the mock sample including wider jets, while the upper end is set by jet measurements only.This range is fully consistent with the estimated range of Galactic BNS merger rates (O'Shaughnessy & Kim 2010;Pol et al. 2020;Grunthal et al. 2021) and the BNS merger rate derived from the detection of gravitational waves (Figure 6; Abbott et al. 2020b).This implies that SGRBs are predominantly the result of BNS mergers with successfully launched relativistic jets, with the most successful direct evidence being the joint detection of the BNS merger GW170817 with its SGRB (GRB 170817A; Abbott et al. 2017a).
An open question remains on the fraction of SGRBs which are derived from NS-BH mergers.Thus far, the identification of EM counterparts to the few known GWdetected NS-BH mergers, GW200105 and GW200115 (e.g., Anand et al. 2021;Dichiara et al. 2021;Rastinejad et al. 2021) have not yielded EM counterparts.Theoretically, one expects the partial or complete disruption of the NS by the BH, resulting in little or no ejected matter and electromagnetic emission (e.g., Foucart 2012).However, there are conditions for which NS-BH mergers can produce successful SGRBs and/or EM emission in simulations (e.g., Bhattacharya et al. 2019;Barbieri et al. 2020;Darbha et al. 2021).Observational studies have also hinted at the possibility of a NS-BH merger con-tribution from the population of SGRBs with extended emission (Gompertz et al. 2020).However, from Figure 6, our SGRB range of rates are 13 times larger than the rates of NS-BH mergers inferred from both GW observations (Abbott et al. 2021) and population synthesis simulations (Sarin et al. 2022).Therefore, our results only support that (at most) a small fraction of SGRBs originate from NS-BH mergers (Figure 6).
In our study, we also uncover a group of SGRBs with inferred jet opening angles of θ j 10 • , which are broader than the derived median of jet opening angle measurements of θ j ≈ 6 • .For most of these cases, the lower limits on the opening angles were inferred thanks to late-time X-ray detections obtained at t j > 10 days.As seen in our study, the inclusion of even a few wide jets pushes the inferred true event rates to lower values, which has consequences for the progenitors of SGRBs.Additionally, successful jets from NS-BH mergers are expected to be wider (θ j ≈ 25 − 30 • ), due to the lower densities of the surrounding environment.In addition, the widest jets are expected to originate from mergers containing highly spinning BHs (Murguia-Berthier et al. 2017;Ruiz et al. 2018).However, jets wider than ≈ 30 • are not expected to routinely survive, resulting in a failed or choked jet (Ghirlanda et al. 2019).Thus, continued X-ray monitoring of cosmological SGRBs to late times, in tandem with monitoring of future GWdetected SGRBs events, will be imperative in constraining the true population of wide jets.

CONCLUSIONS
We have presented a comprehensive compilation of Swift SGRBs discovered between 2005 and 2021, and observed with Chandra and XMM-Newton at late times (δt > 0.8 days).We conclude our findings below: • From the 29 SGRBs in our final sample: 18 were observed by Chandra, 4 by XMM-Newton and 7 by both observatories, resulting in 60 epochs, that we uniformily-analyzed, across all events at δt > 0.8 days.
Including the wide-angle GRB 211106A with an identified break in its radio afterglow light curve, we find a range of jet break times between ≈ 0.1 − 30 days after the bursts, translating to θ j,det = 6.1 • [−3.2 • , +9.3 • ] (68% confidence on the entire distribution).
• From the non-detection of jet breaks for 19 events, we derive lower limits on the opening angles of θ j 0.3 − 26 • , including 12 new limits.Of particular interest are six events with wide inferred opening angles of θ j 10 • .Coupled with two wide-angle events, SGRBs 050724A and 211106A with θ j ≈ 34 • and ≈ 16 • , respectively, we have unveiled a growing population of SGRBs with wide jets.
• We obtain beaming-corrected total true energies between E true,tot ≈ 10 49 − 10 50 , which are consistent with MHD processes as the mechanism of energy extraction to launch jets.
• We derive a range for beaming-corrected true event rates of R true ≈ 360 − 1790 Gpc −3 yr −1 , for which the low end is set by the inclusion of wider jets, and the upper end is set by including jet measurements alone.These rates are fully consistent with the rates of BNS mergers derived from GW events, as well as the rates of Galactic BNS mergers.This aligns with expectations that the predominant progenitor channel of SGRBs is BNS mergers.It is also plausible, although cannot be confirmed given current rate uncertainties, that most BNS mergers produce successful SGRB jets.
• The SGRB event rate is 13 times larger than the GW-derived NS-BH merger rate.Thus, we find that (at most) a small fraction of SGRBs could originate from these mergers.
Our study highlights the importance of the late-time X-ray monitoring of SGRBs in constraining the beaming angles of SGRBs.SGRBs not only exhibit narrow opening angles of ≈ 6 • , but are also capable of launching broader jets with opening angles of 10 • .In recent years, a concerted effort has been made to follow-up SGRB X-ray afterglows beyond 1 day after the burst trigger, when they are expected to be fainter and exhibit breaks.This has resulted in a substantial increase in the number of jet opening angle measurements and meaningful lower limits.Additionally, the exceptional coordination between the most sensitive space-and groundbased telescopes has allowed us to perform broadband monitoring of these events, providing constraints on the basics of bursts energetics and environmental properties.Relative to Swift/XRT, the exceptional sensitivity of Chandra and XMM-Newton has enabled monitoring up to 60 days after the burst in some cases.Moreover, the high spatial resolution of Chandra enables us to disentangle the afterglow from any contaminating X-ray source.The combined capabilities of these observatories have enabled studies to confront different jet launching mechanisms and further understand their neutron star merger progenitors.The next generation of X-ray missions like NewAthena and XRISM will be indispensable to maintain sensitivity to late jet breaks, and thus wider jets.
Our study provides a baseline for the geometries of successfully-launched jets from mergers.The first and only joint detection of GW170817/SGRB 170817A (170817A) with a successfully launched jet (Abbott et al. 2017a;Goldstein et al. 2017;Savchenko et al. 2017) enabled a tight constraint on the opening angle of ≈ 2 − 4 • , and evidence for jet structure (e.g.Lamb & Kobayashi 2017;Alexander et al. 2018;D'Avanzo et al. 2018;Margutti et al. 2018;Troja et al. 2018;Xie et al. 2018;Fong et al. 2019;Hajela et al. 2019).In tandem with the past two decades of SGRB observations, the upcoming observing run (O4) of the Advanced LIGO/Virgo/KAGRA (Abbott et al. 2020a) and beyond will undoubtedly provide a complementary view on the progenitor conditions necessary to launch jets.
description: (1) GRB name.(2) Best-fit model for the X-ray afterglow.The abbreviations "SPL", "BPL" and "TPL" represent the single, broken, or triple power-law models.(3 − 5) Temporal decay indices for each segment of the best-fit model.(6) X-ray spectral index.(7 − 8) Time (in days) of the breaks for the BPL and TPL models.In case of SGRBs with detected jet breaks in their light curves, we identify the time of the jet break (Section 5) with t b1 for the broken power-law model and t b2 for the triple power-law model.Errors are 1σ.* Value derived from Swift photon index (Γ X ).

Figure 3 .
Figure 3. Distribution of the jet break times for 10 SGRBs with jet break detections (grey histogram) derived from the individual jet break time posteriors.The grey arrow indicates the median jet break time of tj = 1.5[−1.2,+7.7] days.The distribution at tj > 20 days is clearly dominated by a peak which represents the late-time detection of the GRB 211106A jet break in the radio band.For the 19 lower limits, the times of the last X-ray detections are indicated with color-coded lower limits which roughly corresponds to the time out to which we are sensitive to jet breaks for each event.A significant fraction (10/19) of the lower limits are found at times of tj > tj , and specifically there are 6 events with last X-ray detections at tj > 10 days.This shows the existence of a SGRB group with jet breaks located at very late times (and potential wide opening angles) as the case of GRB 211106A.

Figure 4 .
Figure4.left: The normalized probability distribution for each of the 10 SGRBs with detected jet opening angles (color-coded open histograms), along with the total distribution of all 10 events (grey histogram).For this sample, we find a median value of θ j,det = 6.1• [−3.2 • , +9.3 • ] (grey arrow).The inferred core jet opening angle for SGRB 170817A, θj,core = 2 • − 4 • (Margutti & Chornock 2021), is indicated as a green bar.right: The total probability distribution for the sample of 19 lower limits (red dashed histogram) in comparison with the total distribution of detected opening angles (dark grey histogram).The red arrow indicates that this distribution represents lower limit values.There is a distinctive group of SGRB lower limits composed of 6 events with wider minimum opening angles of θj,min 10 • .Thus, this indicates a possible population of SGRBs with wider opening angles than described from the distribution of measurements alone.

Figure 5 .
Figure5.The CDFs of true event rates for the opening angle measurement sample (orange solid line) and the mock sample including wide jets (yellow dashed line) assuming an observed local rate of R obs ≈ 10 Gpc −3 yr −1 .Comparing both distributions, we see that a larger number of wider jet opening angles in our sample would lead to a more constrained true event rate median value.

Figure 7 .Figure 8 .
Figure 7.The F-test flow chart followed to the best-fit model of the SGRB X-ray afterglow light curves.

Table 1 .
Basic Properties of Short GRBs with Late-time X-ray Observations − 10 keV) afterglow light curves of the 29 SGRBs in our sample.Swift observations are represented as follows: WT-mode data with light grey thin diamond, and PC-mode data with dark grey circles.Chandra and XMM-Newton detections are depicted with dark blue squares and pink diamonds, respectively.The X-ray flux upper limits (3σ) for both observatories are shown with light blue and light pink triangles, respectively.

Table 2 .
Chandra and XMM-Newton Observations and Spectral Parameters of Short GRBs

Table 3 .
SGRB Sample Light Curves Fitting Parameters

Table 4 .
Inferred Properties, Jet Opening Angles Measurements and Lower Limits of the SGRB Sample Single power law (SPL) is the best-fit model H1: we need a more complex model : Broken power law (BPL) is the best-fit model H 1 : we need triple power law (TPL) model H0: