Comparing Short Gamma-Ray Burst Jet Structure Models

A structured gamma-ray burst jet could explain the dimness of the prompt emission observed from GRB$\,170817$A but the exact form of this structure is still ambiguous. However, with the promise of future joint gravitational wave and gamma-ray burst observations, we shall be able to examine populations of binary neutron star mergers rather than a case-by-case basis. We present an analysis that considers gravitational wave triggered binary neutron star events both with and without short gamma-ray burst counterparts assuming that events without a counterpart were observed off-axis. This allows for Bayes factors to be calculated to compare different jet structure models. We perform model comparison between a Gaussian and power-law apparent jet structure on simulated data to demonstrate that the correct model can be distinguished with a log Bayes factor of $>5$ after less than 100 events. Constraints on the apparent structure jet model parameters are also made. After 25(100) events the angular width of the core of a power-law jet structure can be constrained within a $90\%$ credible interval of width $ \sim9.1(4.4)^{\circ} $, and the outer beaming angle to be within $\sim19.9(8.5)^{\circ}$. Similarly we show the width of a Gaussian jet structure to be constrained to $\sim2.8(1.6)^{\circ}$.


INTRODUCTION
Gamma-ray bursts (GRBs) are pulses of high energy electromagnetic (EM) radiation from astrophysical sources. There is strong evidence that there are two populations of GRBs -long GRBs and short GRBs (sGRBs) (Kouveliotou et al. 1993). sGRBs are generally shorter in duration and have harder spectral energies. It has been long believed that the merger of binary neutron star (BNS) systems are the progenitors of sGRBs, and in August 2017 this was confirmed by the detection of GRB 170817A (Goldstein et al. 2017;Savchenko et al. 2017). The observation was accompanied by a gravitational wave (GW) compact binary coalescence signal, GW170817, by the LIGO and VIRGO scientific collaboration (Abbott et al. 2017a) and therefore became the first joint EM and GW observation (Abbott et al. 2017b) -a landmark in multimessenger astronomy. However, while GRB 170817A was the closest detected sGRB, at a distance of ∼ 40 Mpc, it was also ∼ 10 3 times less energetic than any of the weakest previously observed sGRBs with known redshift. While there are numerous proposed explanations, it is thought that either sGRB energies fall below the current assumed energy range, that GRB 170817A was viewed off-axis or that sGRBs have a non-uniform energy distribution within the solid angle subtended by their jets.
There had been much work on the jet structure of long GRBs prior to the detection of GRB 170817A; Rossi et al. (2002), Zhang & Meszaros (2002) and further work by Rossi et al. (2004), Zhang et al. (2004a) and Lloyd-Ronning et al. (2004) consider a quasi-universal long GRB jet structure, where discrepancies in observed isotropic energies across observations is due to the varying viewing angle of events described by the same underlying jet structure and show that a Gaussian or power-law jet structure can each fully account for all variations in observed long GRBs until that date. Lazzati & Begelman (2005) consider the effect of surrounding shocked matter about the jet, known as a jet cocoon, as a possible cause for a quasi-universal jet structure. Relativistic magnetohydrodynamic simulations have shown a structured jet can naturally form without a cocoon (Aloy et al. 2005;Kathirgamaraju et al. 2019). It has also been shown that radial and angular structuring can evolve from jets that are initially uniform in energy distribution through hydrodynamical interactions with surrounding material in numerical simulations (Xie et al. 2018;Gill et al. 2019).
Determining the true jet-energy distribution will allow both an explanation for the luminosity of GRB 170817A, as well as insight into the astrophysics behind the jet formation caused by the BNS merger. This energy distribution may be described by a shape function y(θ), where θ is the polar angle from the jet axis such that 0 ≤ θ ≤ π/2. Although the electromagnetic flux at the detector is known (which is proportional to y), in order to infer the distribution we need some way of measuring the viewing angle θ v . For simplicity, we may assume that a universal jet structure exists so that by observing sGRBs with various θ v and y values we can build up a picture of the jet structure. θ v may be measured by radio observations of the superluminal motion of the jet over the months following the event (Mooley et al. 2018), by measuring the source size over time (Ghirlanda et al. 2019), as well as from the spectral features of kilonovae (Kasen et al. 2015;Metzger 2017). GW signals can be used to infer the joint probability distribution of source parameters, including both the luminosity distance d L and inclination angle cos ι = J · N where J is the angular momentum vector of the system and N is between the system and the observer. Under the assumption that all sGRBs are generated by BNS mergers, this will give us the information needed to measure the energy emission (from the EM flux and luminosity distance) and the viewing angle θ v = min(ι, π − ι). By observing multiple BNSs with both EM and GW channels, and assuming a common jet model, we can gradually build up a picture of the apparent jet energy distribution function. This is even the case for BNSs that are only detectable from their GWs and without detectable sGRBs; if we can deduce that the sGRB was undetectable due to being viewed outside the confines of the jet the inferred viewing angle can be used to constrain the width of the jet structure. We derive a method to implement this idea, demonstrate it in a simulated dataset and investigate the number of events required to measure jet parameters and perform model selection. Using future joint GW and prompt sGRB detections to investigate the underlying jet structure has been explored in previous work.  and Gupte & Bartos (2018) consider different jet structures to infer the possible rates of joint GW and prompt GRB detections, an approach that has recently been expanded upon by Farah et al. (2019). These studies, as well as this work, consider the apparent jet structure of sGRBs rather than the intrinsic jet structure (see Section 2 for details on the difference between intrinsic and apparent jet structures). Recent work by Biscoveanu et al. (2019) also considers a similar analysis to this work, but uses the joint GW and sGRB detections to probe the intrinsic jet structure, however we present the first use of this method to compare different jet structure models.
In Section 2, we specify the Gaussian and power-law apparent jet structure models that are compared in this work. In Section 3, we discuss the model and state the likelihood and priors used in analyzing the EM and GW data. Section 5 contains the results when given sets of 100 simulated BNS events of either jet structure model. Lastly, we further discuss the results and conclude in Section 6.

JET STRUCTURES
Jet structuring describes the sGRB jet energy distribution over solid angle. We assume that the distribution is axisymmetric and therefore can be described only in terms of a function of the polar angle, θ. Salafia et al. (2015) introduce the concept of apparent and intrinsic jet structure which is adopted here. The intrinsic structure is the angular energy distribution in the frame of the event while the apparent structure is dependent on the observer's frame and so are related by a Lorentz transformation. The difference between these two distributions is explained by relativistic beaming: the apparent area of the source that the observer receives radiation from depends on the Lorentz factor Γ of the jet and the viewing angle, θ v . When Γ 1, the observer receives emission from the whole visible emitting surface of the source (the head of the jet). This can lead to emission being observed even when θ v is outside the confines of the intrinsic jet structure -often termed an off-axis observation (Granot et al. 2002). In the ultrarelativistic limit, with Γ 1, the observed prompt radiation is mostly from the emitting jet material traveling along the observers line of sight θ = θ v , and little from any surrounding jet material (Rhoads 1997). Therefore the intrinsic jet structure depends on both the angular Lorentz factor and energy distribution in the frame of the event. Additionally Beniamini & Nakar (2018) demonstrate that even with a narrow Lorentz factor distribution low gamma-ray production efficiency at large angles suppresses emission and causes narrow jets for long GRBs.
For sources at negligible redshifts, the apparent isotropic equivalent energy of an event viewed at angle θ v can be written as: where F is the observed fluence. For BNS mergers close enough for GW detection ( 400 Mpc), cosmological effects are small. Therefore we assume there is negligible redshifting and that the source spectrum is fully apparent to the detectors, removing the necessity for a cosmological k-correction (Bloom et al. 2001). We describe E iso (θ v ) in terms of the face-on isotropic equivalent energy E iso,0 = E iso (0) and shape function of the apparent structure y(θ): (2) In this proof-of-principle analysis we concentrate on two simple inhomogeneous models to describe y(θ): a Gaussian jet, and a power-law jet, whose apparent structure functions are shown in Figure 1. These models are often used to describe the intrinsic jet structure of long and short GRBs (Zhang & Meszaros 2002;Rossi et al. 2004;Kumar & Granot 2003;Lloyd-Ronning et al. 2004;Salafia et al. 2015;Lamb & Kobayashi 2017;Oganesyan et al. 2019;Mogushi et al. 2019). Assuming the jet is ultrarelativistic with a bulk Lorentz factor Γ > 100 with little variation across θ and constant gamma-ray production efficiency, emission from off-axis observations is less significant, and the apparent structure closely resembles the intrinsic structure . Therefore we adopt these models and use them to approximate the apparent jet structure in this analysis, which is appropriate for the dominant on-axis emission.
Power law jet -The power-law jet describes a structure with a uniform energy distribution until angle θ c , after which it decays as a power law (θ/θ c ) −k with gradient k. This model was initially proposed by Mészáros et al. (1998) to explain the power-law fit to the decay in the afterglow light curve. Initial work in Rossi et al. (2002); Zhang & Meszaros (2002); Lazzati & Begelman (2005) suggested a value of k ∼ 2 for long GRBs to hold with the observed relation Pescalli et al. (2015) shows that k > 2 allows for a better fit to data. However, like Lamb & Kobayashi (2017) and Oganesyan et al. (2019) we assume k = 2 for simplicity. The power-law jet can be further parameterized by including a sharp cut-off at opening angle θ j where θ c < θ j . This is supported by simulations (Aloy et al. 2005;Rezzolla et al. 2011).
The power-law shape function is then: (3) Gaussian jet -The Gaussian jet depends on a width parameter θ σ such that: The Gaussian jet structure was initially proposed in Zhang & Meszaros (2002) as an alternative quasi-universal jet structure model that could explain the relation found in Frail et al. (2001). Since then, Gaussian-like jet structures have been reproduced in simulations (McKinney 2006). Lyman et al. (2018); Troja et al. (2018); Lamb & Kobayashi (2018) later found that a Gaussian jet is preferable to a power-law jet structure in fitting the broad-band GRB 170817A afterglow data.

THE MODEL
We consider a dataset D consisting of the GW strain H and corresponding measured fluence F measurements of N BNS merger events, where D = {F , H}. The analysis is described by parameters Θ and assumes a jet structure model M while I denotes information used to set the priors and form of the likelihood functions. The parameters describing the EM emission in the jet are Θ = {θ M , E iso,0 }, which is the set of structure parameters of model M and the E iso,0 for each event respectively. While ignoring the evidence, the posterior of Θ can be written: where Φ are the set of nuisance parameters, consisting of the luminosity distance and viewing angle of each event, Φ = {d L , θ v } that are marginalized over. This integral can be approximated from S samples of the joint posterior distributions on d L and θ v produced from GW inference:

GRB Likelihood
For each event, it is assumed that the measured fluence F is equal to the fluence from the sGRB F µ with some added background fluence b and Gaussian noise, which is assumed to be drawn from a normal distribution with standard deviation σ n , which is not constant over events. Therefore the GRB likelihood is the product of all events: Values b and σ n are assumed to be known. F µ depends on each event parameters and the assumed model and can be calculated using Equations 1 and 2.

Priors
Priors are placed on E iso,0 for each event and the model dependent structure parameters θ M . The prior on E iso,0 is independent of the assumed jet model and is a log-normal distribution Salafia et al. 2015): with a mean µ E = ln 10 49 and standard deviation σ E = 1. This is a rather narrow energy prior, but we assume most of the variation in the GRB isotropic equivalent energy is due to the jet structuring. For the Gaussian jet, p(θ σ |M, I) = (π/2) −1 with 0 < θ σ < π/2. For the power-law jet, we use a uniform prior p(θ j |M, I) = (π/2) −1 and p(θ c |θ j , M, I) = θ −1 j , with θ c < θ j .

DATA SIMULATIONS
The model is tested on simulated BNS mergers data. A value for d L , θ v and E iso,0 is assigned to each event. Events are assumed to be uniformly distributed within the space constrained by the detector horizon, d max = 400 Mpc, such that p(d L ) = 3d 2 L /d 3 max . The θ v of each event is assumed to be distributed isotropically, and therefore p(cos θ v ) = U (0, 1). Each E iso,0 is sampled from the prior distribution described in Equation 7. The observed GRB fluence of each event is then generated from these variables from the likelihood in Equations 6 by assuming a jet structure model and assigning values to each θ M . The injected model parameters are θ σ = 20 • , θ c = 8 • and θ j = 32 • . The posteriors on d L and cos θ v are approximated by bivariate normal distributions with a covariance found from averaging fitted normal distributions to GW posteriors analyzed with LALINFERENCE -a software library for performing Bayesian inference on compact binary coalescence signals (Veitch et al. 2015). The GW posteriors are produced from simulated BNS injections given to the advanced LIGO and advanced Virgo network at 2019 sensitivity with a network signal-to-noise ratio > 8.
The fluence parameters σ n and b would vary over detections as well between different detectors when considering real events, where they would then have to be determined on a case-by-case basis. For example, given the merger time of a BNS merger signal from a GW detection, we can search for a corresponding sGRB in data measured by Fermi's Gamma-ray Burst Monitor (GBM) using a technique similar to that in Blackburn et al. (2015) and Burns et al. (2019) where the background can be fitted using the method described in Goldstein et al. (2016) and the fluence can be determined by fitting spectral models as described in Gruber et al. (2014) and Yu et al. (2016). For these simulations, the standard deviation on F scales with the signal strength and background so that σ n = k(F µ + b), to approximate Poisson statistics where we have assumed that the photon count scales with the fluence linearly with constant k which is dependent on the detector's effective area and energy band. For these simulations we assume all photons in the 10-300 keV energy band have energies ∼ 150 keV and are received by a detector with surface area of ∼ 300 cm 2 so that k = 7 × 10 −10 erg cm −2 photon −1 . The background is set to be constant across events to b = 3.7 × 10 −6 erg cm −2 , which is an approximation of the background for GRB 170817A integrated over the burst duration for energies ranging from 10 − 300 keV (Goldstein et al. 2017). In reality, the background count rate of the GBM fluctuates over short timescales and a large background count rate increases uncertainties on the fluence measurements of weak sGRB signals.

RESULTS
We test the analysis on two datasets, each of N = 100 BNS events that are simulated as described in Section 4, each with S = 50 samples from the joint d L and θ v posteriors. One dataset D GJ consists of sGRBs produced from the Gaussian jet structure model and the other dataset D PL that from a power-law jet structure model. This is then repeated for 4 realizations of D GJ and D PL .

PyMC3's NUTS sampler (Hoffman & Gelman 2014) is used to calculate the posteriors in Equation 5 and the marginal likelihoods in Equation A8
.
For each dataset ln B PL,GJ is calculated according to Appendix A, as well as for subsets of these datasets. These are shown in Figure 2 where the two lines plot the mean ln B PL,GJ over the 4 different realizations for the subsets of D PL in blue and D GJ in purple. The error bars show 1σ standard deviation between the realizations ln B PL,GJ is defined by Equation A1 so that a positive value indicates the data is best described by the power-law jet structuring model, while a negative value implies the Gaussian jet structuring model. The magnitude of ln B PL,GJ is proportional to the certainty the analysis has of this decision. Figure 2 demonstrates the analysis can successfully identify the correct underlying model with increasing confidence as more events are considered.
The posteriors of the jet structure model parameters are also determined.  Table 1.

CONCLUSION AND DISCUSSION
Modeling the jet structure of sGRBs will allow us to further understand events like GW170817/GRB 170817A. We demonstrate a comprehensive and fully Bayesian way of incorporating the distance and viewing angle posteriors from GW inference with the prompt emission from any accompanying on-axis or off-axis sGRB to examine the apparent jet structure of sGRBs. Model comparison is performed between jet structure models and within 100 such events the power-law and Gaussian jet structure models can be distinguished with a significant log Bayes factor > 5 as shown in Figure 2. Despite the large widths of the GW posteriors, jet structure model parameters can be well constrained after just 25 events for both the Gaussian jet structure and the power-law shown in Figure 3 to within 90% credible intervals of widths ∆θ σ ∼ 2.8 • , ∆θ c ∼ 9.1 • and ∆θ j ∼ 19.9 • . After 100 events the constraints become tighter still with bounds of ∆θ σ ∼ 1.6 • , ∆θ c ∼ 4.2 • and ∆θ j ∼ 8.5 • as seen in Figure 3 and Table 1.
The current BNS rate is predicted to be 110 − 3840 Gpc −3 yr −1 (Abbott et al. 2019a). With a search volume of 2.5 × 10 6 Mpc 3 yr, the number of BNS detections during the third observing run is predicted to be 2 +8 −2 yr −1 and increase to 8 +42 in the fourth observing run after KAGRA (Somiya 2012) is operational and with further detector design improvements, increasing the search volume to 1.3 × 10 7 Mpc 3 yr (Abbott et al. 2019b). The analysis does not require every detection to be accompanied by a sGRB, however the sky localization of the BNS from the GW detection must be within the field of view of an operational GRB detector. A fraction of these events will be shielded by the Earth or occur while satellites transition over the South Atlantic Anomaly and therefore would not qualify for this analysis. Therefore a dataset of ∼ 25 BNS events could be obtainable within this period but it may take until the third generation detectors are built until a dataset of ∼ 100 BNS events is amassed. With the installment of third generation GW detectors (eg. Einstein Telescope (Punturo et al. 2010)), the BNS detection rate will likely increase into the thousands per year. This will be complemented by developments in GRB detectors (eg. THESEUS (Amati et al. 2018)) which will provide improved sensitivity in GRB observations and allow for a more detailed study of the jet structure when viewed outside the jet confines.
This work serves as a proof-of-principle and a starting point for a more inclusive analysis which can be applied to real data. This would include a sGRB likelihood that accounts for the spectral models used in the fluence calculation such as described in Gruber et al. (2014); Yu et al. (2016). Currently the GW distance and inclination posteriors are approximated as normal distributions and any selection effects that are present in real data are neglected. When real non-Gaussian GW posteriors are considered these effects will need to be accounted for. A similar analysis is performed by Farah et al. (2019), where detection statistics are considered to account for such effects.
There are various ways in which the model described here could be expanded and assumptions could be removed. The variation in observed isotropic equivalent energies of GRBs may not be to do with the jet structuring but instead due to large variations of the intrinsic energy of each sGRB's central engine. A wider prior on the face-on isotropic equivalent energy would account for this possibility, as Fan et al. (2017) demonstrated for estimating the luminosity of sGRBs with uniform apparent jet structures. However this would require a much larger dataset to discern jet structure models of sGRBs with non-uniform apparent jet structure. The current analysis could also be modified to account for variations in sGRB jet structure model parameters by placing priors on them and instead inferring the hyperparameters, such as performed by Biscoveanu et al. (2019), who perform a similar analysis. This would also allow for possible correlations between parameters to be distinguished such as the possible anticorrelation between E iso and θ c (Nakar 2019). We assume that the sGRB apparent jet structure can be modeled by a power-law or Gaussian jet structure. However the Lorentz factor and gamma-ray production efficiency are likely much more dynamic than as assumed in this proof-of-principle analysis and may cause the apparent structure to differ from these simple models especially at large angles. In future work we aim to perform model comparison on the intrinsic jet structure and to also model the varying Gamma factor of the source over viewing angle would give deeper insight into the underlying sGRB astrophysics. This can be done by fully accounting for beaming effects by considering a model such as that described in Ioka & Nakamura (2019). Whether or not every BNS merger produces a sGRB is still unknown. Lamb & Kobayashi (2016) suggest that the majority of BNS mergers fail to produce a sGRB if their Lorentz factors are as low as other high-energy astrophysical phenomena such as blazars and active galactic nuclei. The sGRB production rate from BNS mergers can be parameterized and incorporated similarly to Williams et al. (2018). While increasing the model complexity can help answer a number of interesting questions, it increases the required amount of data necessary to be able to perform model comparison. Such a dataset will not be available until third generation GW detectors are operational. Therefore until more data is available it is only feasible to make comparisons between specific models.
We focus on a power-law and Gaussian jet structure for the sGRBs, this analysis could expand to allow for comparison between any pairs of models. Often emission from the surrounding lower energy and Lorentz factor cocoon (Ramirez-Ruiz et al. 2002;Zhang et al. 2004b) has been modeled as a separate component from the central jet in a two-component jet structure (Peng et al. 2005;Lazzati & Perna 2019). There has also been much work into inferring the jet structure from the rise and decay of the observed afterglow (eg. Lamb & Kobayashi (2017) (2018)), which could be used in tandem with this study to provide further evidence when comparing the different jet structure models. This can include off-axis sGRB observations (Lamb & Kobayashi 2018;Beniamini et al. 2020), specifically through features of off-axis afterglows such as X-ray plateaus (Oganesyan et al. 2019;Beniamini et al. 2019).
This integral can be calculated using thermodynamic integration, a technique inspired by elements from statistical mechanics (Gelman & Meng 1998;Lartillot & Philippe 2006). Let p(D|M, I) equate to the partition function of a model M , Z = p(D|M, I), and we denote p(D|Θ, Φ, M, I)p(Θ, Φ|M, I) to be state q. Therefore the probability of the system being in state q is: which is equivalent to Bayes theorem with p = p(Θ, Φ|D, M, I), the joint posterior on Θ and Φ. Now define: q β = p(D|Θ, Φ, M, I) β p(Θ, Φ|M, I), where β is analogous to the thermodynamic definition β ∝ 1/T where T is the temperature. From this definition, q 0 = p(Θ, Φ|M, I) and q 1 = p(D|Θ, Φ, M, I)p(θ|M, I). From Equation A3, Z 0 = 1 and Z 1 = p(D|M, I). Increasing β then represents the transition from the 'uninformed' prior state to the 'informed' posterior state. The partition function holds the same properties as in thermodynamics, and is related to q β by the potential energy U : An expression for ln p(D|M, I) can be made: from the definition of q β . This integral can be approximated by calculating E β [ln p(D|Θ, Φ, M, I)] at a range of discrete 0 ≤ β ≤ 1, where β modifies the likelihood according to Equation A4. We consider 25 temperatures, where β is spread by (i/24) −5 where i = 0, ..., 24 according to Calderhead & Girolami (2009).