Investigating subphotospheric dissipation in gamma-ray bursts using joint Fermi-Swift observations

TL;DR

This study tests a subphotospheric dissipation model for gamma-ray burst emission using joint Fermi and Swift data, finding it explains half of the spectra and suggesting model improvements like including magnetization.

Contribution

It introduces a Bayesian analysis of a localized subphotospheric dissipation model with joint Fermi-Swift data, providing tighter constraints and insights into GRB jet physics.

Findings

16 out of 32 spectra fit the model well

Joint fits reduce estimates of Lorentz factor and luminosity

Model fit quality suggests need for including magnetization

Abstract

The jet photosphere has been proposed as the origin for the gamma-ray burst (GRB) prompt emission. In many such models, characteristic features in the spectra appear below the energy range of the $\textit{Fermi}$ GBM detectors, so joint fits with X-ray data are important in order to assess the photospheric scenario. Here we consider a particular photospheric model which assumes localized subphotospheric dissipation by internal shocks in a non-magnetized outflow. We investigate it using Bayesian inference and a sample of 8 GRBs with known redshifts which are observed simultaneously with $\textit{Fermi}$ GBM and $\textit{Swift}$ XRT. This provides us with an energy range of $0.3$~keV to $40$~MeV and much tighter parameter constraints. We analyze 32 spectra and find that 16 are well described by the model. We also find that the estimates of the bulk Lorentz factor, $\Gamma$, and the fireball luminosity, $L_{0,52}$, decrease while the fraction of dissipated energy, $\varepsilon_{\mathrm{d}}$, increase in the joint fits compared to GBM only fits. These changes are caused by a small excess of counts in the XRT data, relative to the model predictions from fits to GBM only data. The fact that our limited implementation of the physical scenario yields 50\% accepted spectra is promising, and we discuss possible model revisions in the light of the new data. Specifically, we argue that the inclusion of significant magnetization, as well as removing the assumption of internal shocks, will provide better fits at low energies.

Tables

Table 1. Table 1: Summary of data sample, including which GBM detectors were used for the different bursts and for how long we have overlapping data of GBM and XRT. The overlap indicates the overlap between the start of the XRT data and the Fermi T 90 subscript 𝑇 90 T_{90} . Redshifts were acquired from http://www.mpe.mpg.de/~jcg/grbgen.html. We also present the values of the column densities, N H , gal subscript 𝑁 H gal N_{\mathrm{H,gal}} and N H , intr subscript 𝑁 H intr N_{\mathrm{H,intr}} . The PC data interval indicates the time interval after the BAT trigger used when determining N H , intr subscript 𝑁 H intr N_{\mathrm{H,intr}} .

GRB	Redshift	GBM Detectors	$T_{90,\mathrm{GBM}}$	$T_{0,\mathrm{XRT}}$\tablenoteWith respect to the GBM trigger.	Overlap	$N_{H,\mathrm{gal}}$	$N_{H,\mathrm{intr}}$	PC data interval
			(s)	(s)	(s)	(${10}^{22}$ cm-2)	(${10}^{22}$ cm-2)	(s)
080928	1.692	NaI3, NaI6, NaI7, BGO0	14.3	- 27.2	$37$	$0.072$	${0.47}_{-0.16}^{+0.16}$	$4262-12612$
100728A	1.567	NaI0, NaI1, NaI2, BGO0	165.4	134.2	$44$	$0.165$	${3.16}_{-0.42}^{+0.44}$	$1333-7548$
100814A	1.44	NaI7, NaI8, BGO1	150.5	96.1	$55$	$0.018$	${0.00}^{+0.09}$	$5923-34841$
100906A	1.727	NaI8, NaI11, BGO1	110.6	83.8	$27$	$0.353$	${1.24}_{-0.64}^{+0.63}$	$10592-99797$
140206A	2.73	NaI10, NaI11, BGO1	27.3	-5.6	$34$	$0.070$	${1.08}_{-0.44}^{+0.51}$	$126574-956006$
140512A	0.725	NaI0, NaI1, NaI3, BGO0	148.0	108.9	$41$	$0.147$	${0.31}_{-0.05}^{+0.06}$	$4581-12859$
151027A	0.81	NaI0, NaI3, BGO0	123.4	93.4	$18$	$0.037$	${0.58}_{-0.13}^{+0.13}$	$11633-63739$
161117A	1.549	NaI1, NaI2, NaI10, BGO0	122.2	72.5	$53$	$0.043$	${0.99}_{-0.34}^{+0.38}$	$46352-69247$

Table 2. Table 2: The number of analyzed time-resolved spectra for each burst, together with the number of bins that passed the PPC as described in Section  4.4 . The fourth column shows the number of accepted spectra after additional consideration has been taken to the posterior distributions and manual inspection of the fits. The fifth column shows the number of bins which are accepted when we consider GBM data only.

GRB	Analyzed	Accepted (PPC)	Accepted (total)	Accepted (total)
		(GBM+XRT)	(GBM+XRT)	(GBM)
080928	1	0	0	1
100728A	3	1	1	1
100814A	1	1	1	1
100906A	5	4	2	4
140206A	3	0	0	0
140512A	5	4	3	5
151027A	6	5	5	5
161117A	8	6	4	6
$\Sigma$	32	21	16	23

Table 3. Table 3: Point estimates of model parameters and r d subscript 𝑟 d r_{\mathrm{d}} for all 21 fits which passed the PPC in our sample. We also include values of r d subscript 𝑟 d r_{\mathrm{d}} , which is not one of the fit parameters, but obtained from the relation r d = r ph / τ = L ​ σ T / 4 ​ π ​ τ ​ Γ 3 ​ c 3 ​ m p subscript 𝑟 d subscript 𝑟 ph 𝜏 𝐿 subscript 𝜎 T 4 𝜋 𝜏 superscript Γ 3 superscript 𝑐 3 subscript 𝑚 p r_{\mathrm{d}}=r_{\mathrm{ph}}/\tau=L\sigma_{\mathrm{T}}/4\pi\tau\Gamma^{3}c^{3}m_{\mathrm{p}} , as given in Section  2.1

GRB	Time bin
	${\epsilon }_d$	$L_{0,52}$	$\Gamma$	$r_d$
	$N_r$			(${10}^{12}$ cm)
100728A	134.2-150.5	${0.38}_{-0.01}^{+0.01}$	${15.7}_{-0.6}^{+0.6}$	${89.5}_{-1.3}^{+1.3}$	${11.0}_{-0.4}^{+0.4}$	${1.32}_{-0.07}^{+0.07}$
100814A	147.0-150.9	${0.27}_{-0.05}^{+0.05}$	${11.9}_{-1.7}^{+1.4}$	${96.9}_{-3.8}^{+1.5}$	${6.6}_{-0.4}^{+0.5}$	${0.98}_{-0.09}^{+0.09}$
100906A	96.8-100.5	${0.07}_{-0.02}^{+0.02}$	${53.6}_{-5.8}^{+5.8}$	${110}_{-8}^{+8}$	${20.6}_{-2.9}^{+2.9}$	${1.10}_{-0.09}^{+0.08}$
100906A	100.5-105.0	${0.09}_{-0.02}^{+0.02}$	${72.5}_{-7.2}^{+7.0}$	${114}_{-9}^{+8}$	${25.2}_{-3.8}^{+4.0}$	${1.00}_{-0.08}^{+0.08}$
100906A	105.0-110.0	${0.08}_{-0.02}^{+0.02}$	${74.8}_{-7.5}^{+8.0}$	${113}_{-9}^{+9}$	${26.4}_{-4.0}^{+4.2}$	${1.12}_{-0.08}^{+0.08}$
100906A	110.0-115.4	${0.02}_{-0.01}^{+0.01}$	${89.3}_{-10.0}^{+10.3}$	${113}_{-9}^{+9}$	${31.5}_{-4.6}^{+4.7}$	${1.10}_{-0.08}^{+0.08}$
140512A	108.9-114.2	${0.37}_{-0.02}^{+0.02}$	${1.7}_{-0.1}^{+0.1}$	${95.5}_{-1.2}^{+1.2}$	${1.0}_{-0.0}^{+0.0}$	${0.99}_{-0.09}^{+0.09}$
140512A	128.2-142.6	${0.40}_{-0.00}^{+0.00}$	${3.6}_{-0.1}^{+0.1}$	${95.8}_{-0.4}^{+0.4}$	${2.1}_{-0.0}^{+0.0}$	${1.19}_{-0.05}^{+0.05}$
140512A	142.6-145.6	${0.28}_{-0.05}^{+0.05}$	${1.4}_{-0.2}^{+0.2}$	${82.9}_{-5.0}^{+5.1}$	${1.3}_{-0.1}^{+0.1}$	${0.96}_{-0.09}^{+0.09}$
140512A	145.6-156.0	${0.22}_{-0.02}^{+0.02}$	${1.2}_{-0.1}^{+0.1}$	${86.2}_{-2.0}^{+2.2}$	${0.9}_{-0.0}^{+0.0}$	${1.04}_{-0.08}^{+0.08}$
151027A	97.1-105.7	${0.31}_{-0.04}^{+0.04}$	${1.6}_{-0.2}^{+0.2}$	${81.6}_{-3.7}^{+3.6}$	${1.5}_{-0.1}^{+0.1}$	${1.09}_{-0.08}^{+0.08}$
151027A	105.7-109.1	${0.39}_{-0.01}^{+0.01}$	${4.2}_{-0.2}^{+0.2}$	${94.6}_{-1.0}^{+1.0}$	${2.5}_{-0.1}^{+0.1}$	${1.05}_{-0.07}^{+0.07}$
151027A	109.1-112.1	${0.38}_{-0.01}^{+0.01}$	${5.6}_{-0.2}^{+0.2}$	${94.7}_{-1.0}^{+1.0}$	${3.3}_{-0.1}^{+0.1}$	${1.07}_{-0.07}^{+0.08}$
151027A	112.1-115.2	${0.35}_{-0.04}^{+0.04}$	${4.5}_{-0.4}^{+0.4}$	${86.5}_{-2.9}^{+2.8}$	${3.5}_{-0.2}^{+0.2}$	${1.14}_{-0.07}^{+0.08}$
151027A	115.2-120.0	${0.21}_{-0.04}^{+0.04}$	${3.2}_{-0.5}^{+0.4}$	${77.5}_{-6.1}^{+5.0}$	${3.5}_{-0.3}^{+0.3}$	${1.20}_{-0.08}^{+0.08}$
161117A	72.5-90.2	${0.07}_{-0.01}^{+0.01}$	${20.1}_{-2.1}^{+2.3}$	${153}_{-7}^{+10}$	${2.8}_{-0.2}^{+0.2}$	${0.93}_{-0.08}^{+0.08}$
161117A	90.2-92.7	${0.09}_{-0.02}^{+0.02}$	${35.7}_{-5.4}^{+5.6}$	${162}_{-13}^{+13}$	${4.2}_{-0.5}^{+0.4}$	${1.00}_{-0.09}^{+0.09}$
161117A	92.7-100.2	${0.07}_{-0.01}^{+0.01}$	${35.5}_{-2.5}^{+2.8}$	${137}_{-4}^{+4}$	${7.0}_{-0.4}^{+0.4}$	${1.08}_{-0.08}^{+0.09}$
161117A	100.2-104.1	${0.08}_{-0.01}^{+0.01}$	${52.2}_{-4.0}^{+4.6}$	${134}_{-5}^{+6}$	${10.9}_{-0.9}^{+0.7}$	${1.05}_{-0.09}^{+0.09}$
161117A	131.0-138.5	${0.09}_{-0.01}^{+0.01}$	${36.0}_{-2.2}^{+2.2}$	${97.4}_{-0.8}^{+0.8}$	${19.6}_{-1.1}^{+1.1}$	${1.01}_{-0.07}^{+0.07}$
161117A	138.5-144.2	${0.05}_{-0.01}^{+0.01}$	${26.9}_{-2.4}^{+2.5}$	${94.1}_{-1.2}^{+1.3}$	${16.2}_{-1.3}^{+1.3}$	${0.94}_{-0.09}^{+0.09}$

Equations

$$\Gamma =$$

$$400,450,500,600,700,800,900,1000$$

$$L_{0,52}=$$

$${\epsilon }_{\mathrm{d}}=$$

$$0.2,0.25,0.3,0.35,\mathrm{0.4.}$$

$$\mathrm{Pr}(\theta |y)=\frac{\mathrm{Pr}(\theta )\mathrm{Pr}(y|\theta )}{\int \mathrm{Pr}(\theta )\mathrm{Pr}(y|\theta )\mathrm{d}\theta }\propto \mathrm{Pr}(\theta )\mathrm{Pr}(y|\theta ),$$

$$\mathrm{Pr}(\log {\epsilon }_{\mathrm{d}})$$

$$\mathrm{Pr}(\log L_{0,52})$$

$$\mathrm{Pr}(\Gamma )$$

$$\mathrm{Pr}(N_{\mathrm{R}})$$

$$p(y^{\mathrm{rep}}|y)=\int p(y^{\mathrm{rep}}|\theta )p(\theta |y)\mathrm{d}\theta ,$$

$$p_{\mathrm{b}}=\mathrm{Pr}\left(T(y^{\mathrm{rep}})>T(y)|y\right),$$

Full text

gamma-ray burst: general, radiation mechanisms: thermal

††journal: ApJ††facilities: Fermi(GBM), Swift(XRT)††software: astropy (Astropy Collaboration et al., 2013), PyMultiNest (Buchner et al., 2014), Xspec (Arnaud et al., 1999), Scipy (Jones et al., 2001), Matplotlib (Hunter, 2007), Pandas (McKinney, 2010), Seaborn (Waskom et al., 2014)

1 Introduction

The prompt phase of gamma-ray bursts (GRBs) is characterized by strongly variable gamma-ray emission that typically lasts less than a minute. While almost all models agree that this emission originates from internal processes in a relativistic jet, the mechanism producing the emission is not understood. GRB spectra are often described using empirical models, particularly a smoothly broken power law known as the Band function (Band et al., 1993). This function has often been interpreted in terms of synchrotron radiation, see e.g. Tavani (1996); Briggs et al. (1999); Abdo et al. (2009); Zhang et al. (2016). However, it was argued by Preece et al. (1998) that fits with the Band function show that a large fraction of observed spectra are harder than can be accounted for by synchrotron radiation. This has been one of the reasons for considering other emission mechanisms. There are also some GRBs observed with very hard spectra, which can be well described by blackbody or multi-color blackbody models, (Ryde, 2004; Ryde et al., 2011; Ghirlanda et al., 2013; Larsson et al., 2015). Although pure blackbody emission is clearly too hard to describe most observed spectra, there are several possible broadening mechanisms which soften the low-energy slope of photospheric emission, including geometric effects (Pe’er, 2008; Lundman et al., 2013), and subphotospheric dissipation (Rees & Mészáros, 2005; Pe’er et al., 2006; Giannios, 2006; Beloborodov, 2010; Vurm et al., 2011; Chhotray & Lazzati, 2015). Photospheric emission also provides a viable explanation for the Yonetoku relation (Yonetoku et al., 2004, 2010), as shown by Ito et al. (2018). See also Parsotan et al. (2018) for further discussion on photospheric emission as an origin of the Yonetoku relation, as well as for the related Amati (Amati et al., 2002) and Golenetskii (Golenetskii et al., 1983) relations.

To evaluate different physical models for the prompt emission it is important to directly fit the models to data. Indeed, considerable work in recent years have shown that inferences about emission processes based on fits with the Band function, including the aforementioned hardness-problem, can be misleading (Burgess et al., 2014; Burgess, 2017; Burgess et al., 2018; Ahlgren et al., 2019). Examples of previous physical models fitted to data include a physical synchrotron model (fit to GRBs observed by the BATSE instrument; Lloyd & Petrosian 2000), the ICMART model (fit to GRB 080916C; Zhang & Yan 2011), and the external shock model (fit to GRB 141028A; Burgess et al. 2016). Unfortunately, fitting these models to data is generally time consuming, making broad usage difficult. However, there have recently been new development, with Burgess et al. (2018) showing successful fits to a sample of 19 GRBs using a physical synchrotron model.

Another recent example is provided in Ahlgren et al. (2019) (A19 from now on), where we tested a specific photospheric model for localized subphotospheric dissipation by internal shocks with no magnetic fields (DREAM, introduced in Ahlgren et al. 2015; A15 from now on). We analyzed time-resolved spectra of 36 GRBs observed by the Gamma-ray Burst Monitor (GBM; Meegan et al. 2009) on board the Fermi Gamma-ray Space Telescope and found that $\sim~{}30\%$ of the spectra could be well described by the model, with 10 GRBs having at least half of their spectra accepted. The model consistently failed to describe the GRBs with the highest luminosities, which may be a result of the specific dissipation scenario considered. The level of dissipation was found to be around $\sim~{}5\%$, consistent with the internal shock scenario, and the luminosity and Lorentz factor of the model were found to be correlated. The latter has independently been reported in other studies (Lü et al., 2012; Ghirlanda et al., 2018). In our study there was no correlation between the fitted Band function $\alpha$ and $\beta$ parameters and any properties of the fits using the DREAM model. The current implementation of the physical scenario is still being improved upon, and the large number of well-described spectra motivates further exploration of the model scenario.

While Fermi observes GRB prompt emission over a wider energy range than any other GRB mission, it is limited by its lower energy bound of $8$ keV. Physical models (including DREAM) often predict distinct spectral features below this energy and observations of the prompt emission down to soft X-rays therefore have the potential to be very constraining. For instance, in the case of DREAM, the curvature at low energies is a result of incomplete Comptonization of the seed photon blackbody, which may be located at energies as low as $\sim$1 keV in the observer frame. The X-ray Telescope (XRT; Burrows et al. 2005) on board the Neil Gehrels Swift Observatory observes in the $0.3-10$ keV energy band, but the observations are limited by the fact that they typically start $\sim 100$ s after a trigger from the Burst Alert Telescope (BAT). For most GRBs the prompt gamma-ray emission has already ended by this time.

Studies of early XRT observations have shown that the emission is due to a combination of late prompt emission and early afterglow emission from the interaction between the jet and circumstellar medium (O’Brien et al., 2006; Burrows et al., 2007). The former is manifested by flares in the light curve, which have properties similar to the prompt gamma-ray emission (Chincarini et al., 2010). Given these considerations, we note that it is possible to study the prompt emission from soft X-rays to gamma rays in GRBs for which the Swift/XRT light curve is dominated by flares and the prompt gamma-ray emission has a sufficiently long duration.

In previous joint analyses of XRT and GBM prompt emission, Oganesyan et al. (2017, 2018) found that many spectra exhibit a second spectral break at around a few keV. They performed fits using a doubly smoothly broken power-law, and interpret their results in terms of synchrotron radiation. Additionally, Nappo et al. (2017) performed a time resolved joint analysis of data from several instruments, including GBM and XRT, of GRB 151027A, where they detect the presence of a thermal component at low energies. These different results are interesting also from the perspective of photospheric models, as these often have a curvature at energies around a few keV.

In this work we present further investigation of the model presented in A19 by performing joint analysis to data from Fermi/GBM and Swift/XRT using Bayesian inference. The only difference to the model presented in A19 is a small expansion of the parameter space. We analyze a sample of 8 GRBs which have overlapping GBM and XRT observations and known redshifts. The goal is to use the large energy window offered by the joint observations to provide new constraints on the physical scenario and parameter space of our model. We also want to assess the impact of the XRT data in prompt emission analysis.

In Section 2 we describe the physical scenario and how it is implemented as a numerical model. We describe the data sample and analysis in Sections 3 and 4, followed by a presentation of the results in Section 5 and discussion in Section 6. In Section 7 we summarize our findings. Throughout the paper we assume standard $\Lambda$-CDM cosmology using the constants $H_{0}=67.3$, $\Omega_{\lambda}=0.685$, and $\Omega_{\mathrm{M}}=0.315$ (Planck Collaboration et al., 2014).

2 Model

We study a photospheric emission model, in which localized subphotospheric dissipation occurs due to internal shocks. Moreover, we assume no magnetization and ignore off-axis emission and other geometrical effects. The model used in this work is identical to that of A19, apart from a small expansion of the parameter space. For completeness, we here briefly describe the physical scenario and its numerical implementation. For a more detailed description of the model, including validation tests, see A19 and the references given below.

2.1 Physical scenario

The physical scenario we consider is a hot fireball (for an overview, see e.g. Pe’er 2015), with localized subphotospheric dissipation at a moderate optical depth. We assume that a central engine emits a hot plasma of baryons, electrons and photons, at a luminosity $L_{0}=L_{0,52}10^{52}$ erg s*-1*. The outflow accelerates due to the thermal pressure from the photons until it obtains its coasting bulk Lorentz factor, $\Gamma$. We assume a dissipation radius, $r_{\mathrm{d}}=r_{\mathrm{ph}}/\tau=L\sigma_{\mathrm{T}}/(4\pi\tau\Gamma^{3}c^{3}m_{\mathrm{p}})$, where $r_{\mathrm{ph}}$ is the photospheric radius and $\tau$ is the optical depth (see e.g. Pe’er et al. 2006). The internal shock assumption furthermore implies that $r_{\mathrm{d}}=\Gamma^{2}r_{0}$, where $r_{0}$ is the nozzle radius. This relation couples the photon temperature at the dissipation site to $\Gamma$, $L_{0,52}$ and $\tau$. At $r_{\mathrm{d}}$, a fraction $\varepsilon_{\mathrm{d}}\varepsilon_{\mathrm{e}}$ of the bulk kinetic energy of the outflow begins to dissipate to the electrons, and a fraction $\varepsilon_{\mathrm{d}}\varepsilon_{\mathrm{b}}$ to the magnetic fields. The dissipation is assumed to continue until $2r_{\mathrm{d}}$. In the particular scenario considered here, we assume that $\varepsilon_{\mathrm{b}}=10^{-6}$ and $\varepsilon_{\mathrm{e}}=0.9$, i.e. that magnetic fields are negligible and that almost all the energy is dissipated to the electrons.We also assume that the heated electrons are Maxwellian distributed. The remaining $\sim 0.1\varepsilon_{\mathrm{d}}$ energy is considered not dissipated. The parameterization with $\varepsilon_{\mathrm{d}}$, $\varepsilon_{\mathrm{e}}$ and $\varepsilon_{\mathrm{b}}$ is for practical reasons only and that the total efficiency in this case is given by is $\varepsilon_{\mathrm{d}}\varepsilon_{\mathrm{e}}$. Further, we assume that $\tau=35$, which means that we test a scenario where the dissipation occurs moderately deep below the photosphere. Tests have shown that this is often a good approximation since many bursts are largely insensitive to this parameter (see A19 for further discussion).

Since we assume a scenario where the dissipation occurs below the photosphere, the heated electrons will interact with the photon field in a non-equilibrium situation. In our model we consider Compton and inverse Compton scattering, pair production as well as pair annihilation. In principle we also account for synchrotron and synchrotron self-absorption. However, these effects are very small in the case of negligible magnetization.

For localized dissipation in outflows with negligible magnetization, hadronic collisions (Beloborodov, 2010) and internal shocks (Kobayashi et al., 1997; Daigne & Mochkovitch, 1998; Rees & Mészáros, 2005) have been suggested as possible dissipation mechanisms. The scenario we consider here is based on the internal shock scenario (Pe’er, 2015).

To simulate the physical scenario we use the kinetic code by Pe’er & Waxman (2005), which treats all processes described above. Our treatment does not include any spatial effects, such as geometric broadening (Pe’er, 2008; Lundman et al., 2013), or jet hydrodynamics (Lazzati et al., 2009). We also assume that the photosphere is sharp, as opposed to a fuzzy photosphere (Beloborodov, 2011; Bégué et al., 2013). While we restrict ourselves to a specific dissipation scenario, alternative scenarios are also possible, see e.g. Vurm et al. (2011).

2.2 Table model

In order to perform fits with the model we construct an XSPEC compatible table model (Arnaud et al., 1999), which consists of a grid of spectra simulated for different parameter values. Model predictions for parameters between the grid points are obtained by linear interpolation during the fitting. The simulations are costly to perform, which means that we are not able to explore all available parameter space. In A19 we presented a model in 3 dimensions with 891 grid points, consisting of the level of dissipation, $\varepsilon_{\mathrm{d}}$, the luminosity, $L_{0,52}$, and the bulk Lorentz factor, $\Gamma$.

For this study we have expanded the model with additional grid points in $\Gamma$, extending the range down to $50$ from $100$ and from $500$ up to $1000$. We added these grid points to obtain a better coverage of values of $\Gamma$ inferred from observations (Ghirlanda et al., 2018). We have also added one additional grid point in $L_{0,52}$, at $1000$. This grid point has been added to account for the possibility of very luminous bursts with narrow jets. There are no differences in the underlying code used to construct the grid. The model is now spanned by

[TABLE]

As noted in A19, high values of $\varepsilon_{\mathrm{d}}$ would yield a re-acceleration of the outflow, which is not accounted for in our code. Thus, we limit $\varepsilon_{\mathrm{d}}$ to values where the effects of such a re-acceleration are expected to be small. In A19 we demonstrated that the finite resolution of the grid introduces systematic uncertainties of $5-8\%$ in the best-fit parameter values, which is smaller than the typical statistical uncertainties. When creating a table model we also obtain one parameter for the redshift and one for the normalization. The latter is set by the redshift as ${d_{\mathrm{L}}^{2}(z=1)}/{d_{\mathrm{L}}^{2}(z_{\mathrm{obs}})}$, where $d_{\mathrm{L}}(z)$ is the luminosity distance. Both of these parameters are kept fixed in the fits.

In line with previous work, we denote this version of the table model DREAM1.3, using the naming convention introduced in A19.111This naming convention was introduced because of the plan to make the model publicly available. It is then convenient to be able to distinguish between different versions of the model used in different articles. Throughout this work, ‘model’ refers to DREAM1.3, unless otherwise stated. The underlying physical model scenario is referred to explicitly as ‘physical scenario’, to avoid confusion. The physical scenario we are testing is subphotospheric dissipation with localized dissipation at a moderate optical depth in a jet with negligible magnetization.

3 Observations

3.1 Sample selection

We examine all GRBs with a known redshift which have overlapping observations in the Swift XRT and Fermi GBM detectors, up until 2018-11-01. As mentioned in Section 2.2, the redshift is needed in order to obtain the luminosity distance for the model. We also require the overlapping interval to be at least 5 seconds long and that it is possible to bin the XRT data into at least two time bins using the method outlined in Section 4 below. Finally, we require that at least one time bin has a signal-to-noise ratio (SNR) in one GBM NaI detector of at least $3$, as described in Section 4 (the SNR in the XRT is always higher than this). Note that we do not use data from the BAT, also on board the Swift satellite, due to its narrower energy range ($15-150$ keV) compared to GBM. Contrary to the burst sample in A19 we do not perform a fluence cut in this sample selection. These criteria result in a relatively small sample of 8 GRBs, the properties of which are summarized in Table 1. In Fig. 1, we show the count rate light curves of all bursts, including both GBM and XRT data.

XRT observations typically start $\sim 100$ s after the GBM trigger (corresponding to the time it takes for Swift to slew following a BAT trigger). For most GRBs in the sample the GBM and XRT data overlap for $40-50$ s at the end of the prompt emission. However, there are two exceptions: GRB 080928 and GRB 140206A. These bursts have XRT data from the start of the GBM trigger. This is because BAT triggered on precursors at $-204$ s and $-56$ s, respectively, relative to the GBM trigger (not included in Fig. 1). Fig. 1 also shows that all the XRT light curves contain flares, as expected if the emission belongs to the prompt phase. In most cases the GBM and XRT light curves are clearly correlated. However, we note that this correlation is less prominent in GRB 100814A and GRB 100906A. This could be due to that the GBM data are particularly weak in these intervals. Additionally, in GRB 151027A, there is a clear delay in the XRT with respect to the GBM.

Comparing with the sample in A19, six bursts (GRB 100728A, GRB 100814A, GRB 100906A, GRB 140206A, GRB 140512A, and GRB 151027A) are contained in both samples. In A19 we found that only two of these bursts, GRB 140512A and GRB 151027A, had more than half of their analyzed time bins well described by our model. We discuss this further in section 6.3.

3.2 Data reduction

3.2.1 Fermi data

In the Fermi analysis we use data observed with the GBM. Specifically we use the Time-Tagged Event (TTE) data from both the NaI and BGO detectors. We include up to three NaI detectors with an angle of incidence less than $60\degree$ in the analysis, as well as the BGO detector with the lowest angle of incidence (Bhat et al., 2016). However, there is one exception where we have excluded an additional NaI from the analysis. In GRB 151027A we use only the NaI0 ($19\degree$) and NaI3 ($37\degree$) detectors, ignoring the n6 ($25\degree$) detector, which has a significantly different spectrum in several bins (example in Fig. 2). This may be due to blockage of the detector by a part of the satellite, or some other issue with the detector response (Goldstein et al., 2012).

The background is determined from a polynomial fit to the light curve, which gives a model for the background with Gaussian errors. We account for temporal evolution of the background spectrum and the changing position of the spacecraft by using rsp2 response files when available.

3.2.2 XRT data

For the spectral analysis of XRT data we use observations taken in the Windowed Timing (WT) mode, which is used when the count rate is high. We also download late time Photon Counting (PC) data for each GRB, in order to determine the intrinsic column density ($N_{\mathrm{H,intr}}$) for the absorption at low energies (see Section 4.3).

The WT data are downloaded as locally-reprocessed data from the UK Swift Science Data Centre XRT GRB repository222http://www.swift.ac.uk/xrt_live_cat/. We create time bins of the spectrum locally, as described in Section 4. These time bins are then used to specify the limits on time-sliced spectra, which we create and download from the online repository. All spectra are grouped such that each energy bin contains at least one count. This the recommended approach when fitting with the cstat fit statistic in XSPEC333see the XSPEC manual appendix https://heasarc.gsfc.nasa.gov/xanadu/xspec/manual/XSappendixStatistics.html. The background is supplied online when downloading the data. The background spectrum is constructed by sampling in an area around the burst position and it is assumed to be Poisson distributed.

All XRT spectra were checked for known calibration issues following Valan et al. (2018). This includes redistribution issues, which may cause a bump below 1 keV and/or a turn-up below $\sim$0.6 keV. Pile-up is automatically dealt with by the online tool. We find redistribution issues in GRB 100728A and GRB 151027A. We accommodate these by ignoring channels below $0.6$ keV in these bursts.

4 Data Analysis

We perform a time resolved spectral analysis using time bins defined by a Bayesian blocks binning (Scargle et al., 2013) of the XRT data. We use battblocks444Heasoft version 6.17 with default settings to bin the XRT light curve. The GBM data are binned in matching time bins. We use HEASARC’s online tool xTime555https://heasarc.gsfc.nasa.gov/cgi-bin/Tools/xTime/xTime.pl to convert between Swift mission time and Fermi mission time.

We calculate the SNR of each spectrum as described in Vianello (2018). For the GBM data we use the Poisson-Gaussian significance, whereas we for the XRT data use the Poisson-Poisson significance. The XRT data consistently have a high SNR, but as in A19 we choose to apply an SNR cut to the GBM data in order to only analyze spectra which contain a significant signal in GBM. Thus, we only analyze spectra with SNR$>3$ in the brightest NaI detector. This leaves us with 32 out of 51 time-resolved spectra to analyze.

4.1 Fitting

We set up the analysis as a Bayesian inference procedure. Bayes’ theorem states that the posterior probability is

[TABLE]

where $\theta$ are our model parameters and $y$ the observed data. $\Pr(\theta)$ is the prior, $\Pr(y|\theta)$ the likelihood, and the denominator the marginalized likelihood (also referred to as the evidence). We use PyMultiNest (Buchner et al., 2014), a python implementation of MultiNest (Feroz & Hobson, 2008; Feroz et al., 2009, 2013), to sample from the model posterior using 600 live points. We have chosen this particular number by testing the analysis with different number of live points to ensure stability.

For the Fermi data, we consider the energy range 8-1000 keV and 200 keV - 40 MeV for the NaI and BGO detectors, respectively. We also ignore the interval 30-40 keV in the NaI detectors, because of the iodine $K$-edge (Bissaldi et al., 2009). For the XRT data we consider the nominal energy range 0.3-10 keV. However, we lower the high-energy limit in cases when the signal stops below 10 keV. Additionally, the low-energy limit is modified in the presence of calibration issues, as described in Section 3.2.2 (which affects GRB 100728A and GRB 151027A).

For the Fermi data we use a likelihood for a Poisson distributed signal with Gaussian background. In XSPEC the corresponding statistic is known as pgstat (for a description of pgstat, see e.g. Burgess et al. 2018). For the XRT data we adopt the Cash statistic (Cash, 1979), for data with a Poisson signal and background. In XSPEC this statistic is referred to as cstat. Note that pgstat and cstat denote the log likelihood, $\log\left[\Pr(y|\theta)\right]$, for the respective data sets. For the joint analysis we consider the fit statistic as the sum of the pgstat and cstat statistics. We perform the analysis using PyXspec, a python implementation of HEASARC’s XSPEC 12.8.1g (Arnaud, 1996).

The GBM detectors are calibrated with a relative uncertainty on the order of 10 % in effective area (Bissaldi et al., 2009). In A19 we found that the best-fit parameter values did not change significantly whether or not we allowed for an effective area correction between the GBM detectors. There is no information on the possible difference in effective area calibration between GBM and XRT. However, we expect that the uncertainty is greater between GBM and XRT than between the GBM detectors. We thus keep the relative normalization between the GBM detectors fixed to unity, whereas we introduce a free relative normalization parameter, $N_{\mathrm{R}}=\mathrm{norm}_{\mathrm{XRT}}~{}\mathrm{norm}_{\mathrm{GBM}}^{-1}$, between the GBM and XRT data (also referred to as the cross calibration constant). We let this parameter be fit separately in each time bin.

4.2 Priors

We choose most of our priors to be uninformative. For the luminosity, $L_{0,52}$, and the level of dissipation, $\varepsilon_{\mathrm{d}}$, we choose log-uniform priors. This means that each decade in these parameters correspond to an equal prior probability. For the Lorentz factor, $\Gamma$, we choose a uniform prior.

Finally, for the effective area correction between the two instruments, $N_{\mathrm{R}}$, we choose a normal distribution prior centered around $1$ with a standard deviation of $0.1$. Thus, we have

[TABLE]

4.3 Absorption

When analyzing X-ray data below 2 keV, Galactic as well as intrinsic (extra-Galactic) absorption becomes relevant. We account for this absorption by using the multiplicative XSPEC models tbvarabs (for Galactic absorption) and ztbabs (for intrinsic absorption). We use the values of the solar abundance vector from Wilms et al. (2000) and the cross-section values listed in Verner et al. (1996). We obtain the weighted total Galactic column density, including the molecular component, $N_{\mathrm{H,gal}}$, from the Swift $N_{\mathrm{H,tot}}$ online tool666http://www.swift.ac.uk/analysis/nhtot/index.php (Willingale et al., 2013). In the case when the fraction of molecular hydrogen lies in the range of 10-30% we have replaced the tbvarabs model by tbabs, where the fraction of molecular hydrogen is fixed to 20%. This is the case for GRB 080928, GRB 151027A, and GRB 161117A.

We determine the intrinsic absorption by fitting late time PC data with a tbabsztbabspow model. At late times the XRT data are dominated by afterglow emission. These data are fainter and are thus captured in PC mode. Pure afterglow spectra are typically well described by power laws, and at late times we expect little spectral evolution (e.g. Racusin et al. 2009). Thus, this method is a good way to determine the intrinsic absorption in a way that does not introduce any degeneracy between the spectral model we wish to test and the absorption. For each burst we create a time averaged spectrum consisting of as late PC data as possible. In order to avoid spectral evolution we also require the light curve at these times to be well described by a power law without any breaks. We use the light curve fits from the online catalogue777http://www.swift.ac.uk/xrt_live_cat/ to choose these time intervals, which we present in Table 1. The fitting to determine $N_{\mathrm{H,intr}}$ is performed using XSPEC’s native Maximum-Likelihood scheme.

Both $N_{\mathrm{H,gal}}$ and $N_{\mathrm{H,intr}}$ are assumed constant for the duration of the burst and kept fixed in all fits with the DREAM model. In Table 1 we summarize the values of $N_{\mathrm{H,gal}}$ and $N_{\mathrm{H,intr}}$ used. Comparing to the $N_{\mathrm{H,intr}}$ distributions of Campana et al. (2012), we note that the value for GRB 100728A is at the high end of the distribution. Additionally, for GRB 100814A, we find values of $N_{\mathrm{H,intr}}$ consistent with 0.

4.4 Posterior predictive checks

Since we are performing a Bayesian analysis we can use posterior predictive checks (PPCs; Gelman et al. 1996; Meng 1994; Lynch & Western 2004; Gelman & Shalizi 2013) to assess the quality of our fits. We draw replicated data from the posterior predictive distribution (PPD) using the XSPEC command ‘fakeit’. We then use this new data to assess the quality of our spectral fits. The PPD is the probability of observing some replicated data, conditioned on the observed data, and can be written as

[TABLE]

where $\theta$, $y$, and $y^{\mathrm{rep}}$ are the model parameters, observed data, and replicated data, respectively. $p(\theta|y)$ is the posterior from which we sample using MultiNest. Thus it is easy to construct a posterior predictive $p$-value ($ppp$-value),

[TABLE]

which corresponds to the classical $p$-value averaged over the posterior, $p(\theta|y)$, and where $T$ is a test statistic (Rubin, 1984). We let $T$ be the fit statistic and calculate $p_{\mathrm{b}}$ for each fit based on $1000$ realizations from the PPD. We consider a fit rejected if $p_{\mathrm{b}}<0.05$.

We stress that an accepted fit only indicates that we cannot reject the fit at the given significance level. It does not mean that the model necessarily can fully describe the data. Additionally, the $ppp$-value of a rejected fit does not tell us how or where the model fails to describe the data. There are many variants of PPCs which can be used to assess the model fitness. We complement our current choice of PPC with manual inspection of posteriors and fits.

5 Results

In this section we present the results of the Bayesian analysis described in Section 4. When performing Bayesian inference we prefer to consider the posterior in its entirety. However, it is often convenient to also use point estimates, here particularly when comparing to the results of A19 and to give an overview of the results in a table. Thus, for point estimates we use the mean of the marginalized posterior for the parameter in question. The associated uncertainties correspond to the 1$\sigma$ credible interval centered around the mean, symmetrical in terms of cumulative likelihood. In Fig. 3 we present examples of corner plots from an accepted and rejected fit, respectively. In Fig. 4 we show the corresponding fits in data space. Corner plots showing the full posterior of all fits are available as online material. Additionally, online we also provide plots corresponding to those in Fig. 4 for all accepted fits.

5.1 Accepted fits and parameter estimates

A total of 32 time bins were analyzed (17 spectra were not analyzed due to SNR$<3$ in the GBM), with 21 bins being accepted under the posterior predictive checks presented in Section 4.4. In Table 2 we summarize the number of accepted bins for each burst in our sample. We also provide an additional column of spectra which are accepted after further examination of the posterior and fits. This leaves a total of $16$ accepted spectra, as presented further in Section 5.2. In Table 5.1, we present the point estimates of all fits accepted under the PPC. We note that GRB 140512A, GRB 151027A, and GRB 161117A have at least half of the analyzed time bins accepted while having more than 1 analyzed spectrum (see Fig. 1).