A Bayesian Inference Framework for Gamma-Ray Burst Afterglow Properties

In the field of multi-messenger astronomy, Bayesian inference is commonly adopted to compare the compatibility of models given the observed data. However, to describe a physical system like neutron star mergers and their associated gamma-ray burst (GRB) events, usually more than ten physical parameters are incorporated in the model. With such a complex model, likelihood evaluation for each Monte Carlo sampling point becomes a massive task and requires a significant amount of computational power. In this work, we perform quick parameter estimation on simulated GRB X-ray light curves using an interpolated physical GRB model. This is achieved by generating a grid of GRB afterglow light curves across the parameter space and replacing the likelihood with a simple interpolation function in the high-dimensional grid that stores all light curves. This framework, compared to the original method, leads to a $\sim$90$\times$ speedup per likelihood estimation. It will allow us to explore different jet models and enable fast model comparison in the future.


Introduction
The gravitational wave emitted from a merging binary neutron star, GW170817 [1], followed swiftly by a Gamma-ray burst (GRB) of short duration known as GRB170817A [2,3] commenced the era of multi-messenger astrophysics (MMA). This is the first direct evidence showing neutron star mergers as the progenitor of short GRBs [4]. The union between electromagnetic (EM) and gravitational wave (GW) observations allow for astrophysical objects to be investigated with alternative probes, enabling the source properties to be studied under a new light.
The coalescence of compact objects with at least one neutron star is of particular interest. By the end of the third LIGO observing run (O3), one more binary neutron star (BNS) event, GW190425 [5], and two NS-BH mergers are reported [6]. However, GW170817 remains the only event with a confirmed GW-EM detection and has allowed a major leap in terms of its scientific implications.
A broadband search of the EM counterpart of GW170817 shows both thermal kilonova emission (AT2017gfo) generated from the radioactive decay of newly synthesized elements in the ejecta and a non-thermal synchrotron component from the relativistic jet [7]. The latter is the GRB afterglow that typically shines in X-ray, optical, and radio wavelengths and can last for months to years following the GRB. For the case of GRB170817A, the first confirmed X-ray counterpart was observed by the Chandra X-ray Observatory and announced at around 9 days post-merger [8]. Late-time monitoring of this event unveils the X-ray afterglow emission reaching its peak luminosity at around 155 days [9,10] and continues to shine at more than 1000 days after the GRB [11][12][13].
Modeling the afterglow light curves would shed lights on the physical environment of the system, as the afterglow peak time is dependent on the observers' viewing angle and the angular profile of the jet. Studies of the temporal behavior of GRB170817A afterglow suggest a structured jet seen off-axis, see, e.g., in [13][14][15][16][17][18]. Understanding of the jet structure is essential to explain the observations of off-axis BNS events such as GRB170817A and uncovering the underlying physical mechanisms behind their production. With better modelling of the jet structure, we are able to more confidently infer the rate of GRB production [19], and the inference of the viewing angle of individual MMA events improves, allowing for tighter constraints on inference of the Hubble constant [20,21].
In this paper, we describe a semi-analytical GRB afterglow model that takes into account the effects of lateral spreading in Section 2. Based on this model, we develop a Bayesian framework that would allow quick parameter estimation and explain it in Section 3. We present the results with this methodology and discuss its future application in Section 4.

GRB Afterglow
A broadband afterglow lasting hours-to-days typically follows a GRB, where this GRB afterglow is produced in the shock system that forms as the ultra-relativistic jet that emitted the GRB decelerates. Using energy conservation for a spherical blastwave, the dynamical behaviour of the decelerating system can be modeled, and the change in Lorentz factor, dΓ, with the change in swept-up mass found, see, e.g., in [22]. The instantaneous Γ and sweptup mass of the blastwave can be used to find the synchrotron emission responsible for the broadband GRB afterglow. By following the method in [23], a single jetted outflow with a half opening angle, θ j , is split into multiple emission components and the dynamics and synchrotron flux, relative to the observers line-of-sight from each segment, can be found. Summing across the multiple components at a given observer time, t, is equivalent to integrating across the equal arrival time surface for the jet, and by allowing each component to have a unique energy and Γ, then the afterglow from complex jet structures can be found, see, e.g., in [24].
This method naturally accounts for the edge-effect in a GRB afterglow, where the jet edge becomes visible for an on-axis observer when Γ(t) < 1/θ j , and the resulting jet-break in the lightcurve. We further add the effects of lateral spreading in the jet by considering the maximum sound speed of each component and the change in radius due to this lateral spreading (see in [23,25]). With this method, the afterglow from spreading and non-spreading jets, with a defined structure profile can be found for observers at any viewing angle i.e., within the jet opening angle for cosmological GRBs, see, e.g., in [26], and at higher viewing angles for gravitational wave counterparts to neutron star mergers, see e.g., in [15,16].
The semi-analytic approach used in generating the lightcurves divides the jet emission area over its radial and angular components into a grid of a given resolution. The dynamics and flux are solved for each resolution element within the jet individually before summing to give the resultant lightcurve. The duration of each afterglow iteration depends on the resolution, where for off-axis events the resolution can be set to low values >50, while for events with low viewing angle, wide opening angles or high Γ, the required resolution becomes much higher, and typically >10,000 to maintain the same level of accuracy. Fitting algorithms can require 100s of thousands of individual model lightcurves for a reliable fit. In our test model with four free parameters, it requires at least ≥1 ×10 4 iterations for a single analysis. A model that takes several seconds, or even minutes, per iteration is therefore not ideal.

Jet Structuring
For a power-law jet, the energy and Lorentz factor distribution of the ejecta with respect to jet central axis is scaled as In this scenario, the jet is parameterized by two characteristic angles, θ c and θ j . Starting from the center of the jet, the energy is uniformly distributed within the core, θ c . Outside of θ c , the energy distribution follows a power-law decay with a sharp decline at the edge, θ j . This structure was proposed to account for the observed power-law decay in GRB afterglow light curves. Here, we assume the power-law index k = 2 as used in Lamb and Kobayashi [24].

Parameter Estimation
The model described in Section 2 is used to simulate X-ray afterglow light curves from 0.15 to 1000 days after the GRB event. Based on the X-ray fluxes, we generate random noise assuming a minimum signal-to-nose ratio (SNR) of 3. The injection data are the sum of both the signal and the noise components. We perform the parameter estimation using Bayesian inference, from which the posterior distribution of jet parameters is given by, according to Bayes theorem, where L(d|θ) is the likelihood function of the data given a set of model parameters θ, and π(θ) is the prior distribution for θ and Z is the evidence, or the observed data. A nested sampling algorithm is used to calculate the evidence and the posterior densities. We carry out the analysis with Bilby, a Python-based Bayesian inference library [27], and the Dynesty sampler [28]. However, one of the issues with stochastic sampling over a complicated physical system is that the likelihood evaluation in the process are computationally expensive. In order to increase the efficiency in investigating different jet models, we employ an interpolated GRB model and substitute the original likelihood function with it. This is done by first simulating each model parameter across a certain range in the parameter space and storing the resulting light curves in a multidimensional grid. The detailed configuration of the grid is specified below: where θ j is the jet half-opening angle, θ c is the core width of a structured jet, θ obs is the angle between the observer and the jet central axis, and 0 is the isotropic equivalent energy of the jet. For each dimension, the parameter range is linearly cut into subsets. We increase the density of splits for parameters that have higher impacts to the resulting light curves . The total number of entries in our 4D Cartesian grid is then 0 × θ j × θ c × θ obs = 40 × 50 × 40 × 50 = 4,000,000. Figure 1 shows how afterglow light curves distribute in the three-dimensional space for an observer at different viewing angles. When the inclination angle increases, the afterglow flux peaks at later times. Jet parameters in this dataset are θ j = 0.2 rad (∼11.5 • ), θ c = 0.15 rad (∼8.6 • ) and 0 = 2 × 10 52 ergs −1 . To put the emphasis on the jet structure and save computational power, we keep the rest of the afterglow parameters constant throughout the simulation. The fixed parameters include the initial bulk Lorentz factor Γ = 100, the ambient particle number density n = 0.001 cm −3 , the shock accelerated election index p = 2.15, the microphysical parameters B = 0.01 and e = 0.1 [24], and the luminosity distance d = 40 Mpc [1] for a GW170817-like source . The resolution in the jet model is set to 50 for all light curves.
We adopt uniform priors with upper and lower bounds the same as the simulation range for every parameter of interest. At each Monte Carlo sampling point, instead of calculating the likelihood from the original model, it performs a linear interpolation between the adjacent parameter values in the grid. A comparison between the interpolated and theoretical light curves with same parameter values is shown in Figure 2. The largest difference appears at observing angles near the jet edge, i.e., slightly off-axis case with~10% deviation. As the geometry between the jet and observers affects the light curves the most, we set the numbers of splits for θ j and θ obs larger than other parameters. Increasing the density of splits in the grid would raise the accuracy but the simulation time also dramatically grows. For postpeak observing epochs and the rest of off-axis scenarios, the interpolated light curves well resemble the theoretical values with a <10% deviation. As our input data has a SNR value of 3, the error caused by the lack of entries in the interpolation process is not the dominant source of uncertainty.

Results and Discussion
The input data, produced with the method described in Section 3, are the X-ray light curve generated from a GRB jet with a power-law profile. The width of the jet is set to θ j = 0.3 rad with the core angle being θ c = 0.15 rad. Considering an off-axis observer, we set the viewing angle θ obs = 0.36 rad (21 • ) and the isotropic energy is 0 = 2 × 10 52 erg s −1 . The remaining parameters have the same values as those used to simulate the 4D grid. This light curve is composed of the 0.4 keV fluxes at observing times from 0.16 days to 1000 days post-event and visualized as blue dots in Figure 3. After incorporating the interpolated model in the sampling process, the time required for calculating the likelihood per iteration reduces from 7.78 seconds to ∼0.09 seconds with only 1 CPU. The overall run time with this method decreases by ∼90 times. Figure 4 shows the joint posterior of the four parameters mentioned above. The injection values are well recovered within the 68% confidence interval. This serves as a proof-of-principle that this method is not only efficient, but also capable of capturing the jet features if the afterglow emission is intrinsically produced by such a system.
In this work, we perform an analysis on simulated light curves with a uniform signalto-noise ratio throughout all observing epochs. In reality, for the same patch of sky and the same observational instrument, the background fluence would remain at a similar scale. The temporal evolution of a GRB afterglow, which arises until reaching the peak luminosity at t peak and then declines with time, yields a varying observational significance. The X-ray emission from the position of GW170817, except for the first observation at 2.3 days, reports a significance of ≥3σ at all epochs. We set all SNR equaling to the detection limit for a conservative estimate. Recent studies show that the X-ray counterpart of GW170817 is still bright after 1200 days [11,12]. The 3σ excess compared to the prediction from an off-axis structured jet poses a challenge to previous agreements from multi-wavelength studies. It also makes modeling the afterglow of GRB170817 worthwhile even after four years as new evidence is still emerging.  On the other hand, we use relatively wide uniform priors for jet structure parameters in the analysis. One of the strengths of multi-messenger astronomy is that we can combine the information obtained from various approaches and place tighter constraints on the source properties. As gravitational wave observation provides an independent measure of viewing angle and luminosity distance, the posteriors from GW parameter estimation can be passed to the afterglow analysis as prior function and thus increase the accuracy in search of the GRB afterglow parameters.
In this paper, we present the validity of a Bayesian framework combined with an interpolated GRB model. It effectively lessens the computational cost while being capable of constraining jet structure parameters in low SNR circumstances. In the upcoming era of a global gravitational wave detector network (LIGO-VIRGO-KAGRA) and new EM observatories like the Gravitational Wave High-energy Electromagnetic Counterpart All-sky Monitor (GECAM), it is foreseeable that we will find more GW-EM counterparts in the near future. A fast, functional parameter estimation will be beneficial for modeling the GRB afterglow light curves. With this methodology, more jet models, e.g., Gaussian or Top-hat jets, will be investigated and its application on the real data of GW170817, the only event so far with GW-EM detection, will be included in our future works.

Conflicts of Interest:
The authors declare no conflict of interest.