Non-dissipative photospheres in GRBs: Spectral appearance in the Fermi/GBM catalogue

TL;DR

This study investigates whether a significant portion of GRB spectra can be explained by non-dissipative photospheric emission, finding that over a quarter of observed spectra are consistent with this simple model, challenging previous assumptions.

Contribution

It demonstrates that a substantial fraction of GRB spectra are compatible with non-dissipative photospheric emission, using synthetic spectral modeling and empirical fits to Fermi/GBM data.

Findings

Over 25% of GRB spectra are consistent with non-dissipative photospheric emission.

Empirical models tend to soften the spectral index compared to theoretical predictions.

The fraction of compatible spectra could increase with sub-photospheric dissipation.

Abstract

A large fraction of gamma-ray burst (GRB) spectra are very hard below the peak. Indeed, the observed distribution of sub-peak power-law indices, $\alpha$, has been used as an argument for a photospheric origin of GRB spectra. Here, we investigate what fraction of GRBs have spectra that are consistent with emission from a photopshere in a non-dissipative outflow. This is the simplest possible photospheric emission scenario. We create synthetic spectra, with a range of peak energies, by folding the theoretical predictions through the detector response of the FERMI/GBM detector. These simulated spectral data are fit with typically employed empirical models. We find that the low-energy photon indices obtain values ranging $-0.4 < \alpha < 0.0$, peaking at around $-0.1$, thus covering a non-negligible fraction of observed values. These values are significantly softer than the asymptotic value of the theoretical spectrum of $\alpha \sim 0.4$. The reason for the $\alpha$-values to be much softer than expected, is the limitation of the empirical functions to capture the true curvature of the theoretical spectrum. We conclude that more than a $1/4$ of the bursts in the GBM catalogue have at least one time-resolved spectrum, whose $\alpha$-values are consistent with a non-dissipative outflow, releasing its thermal energy at the photosphere. The fraction of spectra consistent with emission from the photosphere will increase even more if dissipation of kinetic energy in the flow occurs below the photosphere.

Tables

Table 1. Table 1: The percentages of simulated NDP spectra that can explain the α 𝛼 \alpha -values in different subsamples in the the GBM catalogue.

Simulation	Sample	within FWHM (%)	(%)	$\alpha$ $>$ ${\alpha }_{\min }$ (%)
NDP coasting phase	Full catalogue	12	28	38
NDP accelerating phase	Full catalogue	0	1	3
NDP coasting phase	Cluster 1	14	34	43
	Cluster 2	2	2	9
	Cluster 3	20	43	54
	Cluster 4	0	0	0
	Cluster 5	9	30	43

Table 2. Table 2: Distributions of α 𝛼 \alpha for synthetic photospheric spectra (coasting phase and no dissipation). Simulations were made with E pk subscript 𝐸 pk E_{\rm pk} -values sampled from the distributions in the GBM catalogue, and in clusters 1 to 5 defined in Acuner & Ryde ( 2018 ) . The parameters are the medians and inter-quantile ranges (IQR).

Sample	${\alpha }_{median}$	IQR
Full GBM	-0.16	0.15
Cluster 1	-0.13	0.044
Cluster 2	-0.13	0.09
Cluster 3	-0.18	0.068
Cluster 4	-0.2	0.113
Cluster 5	-0.08	0.043

Equations

$$N_{\mathrm{E}}=K{\left(\frac{E}{E_{\mathrm{pivot}}}\right)}^{0.4}{e}^{-{\left(\frac{E}{E_{\mathrm{c}}}\right)}^{0.65}},$$

$$N_{\mathrm{E}}=K{\left(\frac{E}{E_{\mathrm{pivot}}}\right)}^{0.66}{e}^{-\left(\frac{E}{E_{\mathrm{c}}}\right)}$$

$$p(\theta \mid y)=\frac{p(\theta )p(y\mid \theta )}{p(y)}$$

$$p(y_{rep}\mid y)=\int p(y_{rep},\theta \mid y)d\theta .$$

Full text

††pagerange: Non-dissipative photospheres in GRBs: Spectral appearance in the Fermi/GBM catalogue–14††pubyear: 2019

1 Introduction

It has been argued that the observed characteristics of gamma-ray bursts (GRBs) can be explained by emission from the photosphere in a relativistic outflow. Two possible scenarios have been suggested. First, the photospheric emission could be accompanied by a non-thermal component (e.g., Mészáros et al., 2002; Ryde, 2005; Axelsson et al., 2012; Guiriec et al., 2013; Giannios & Uzdensky, 2019) in a hybrid scenario producing the range of spectral shapes observed. Prominent examples of such cases are GRB090902B (Ryde et al., 2010) and GRB190114C (Wang et al., 2019, Fermi collaboration, in prep.). Alternatively, the entire emission could be from the photosphere. In such a case, dissipation of the flow kinetic energy close to the photosphere is required in order to broaden the spectra to the observed shapes (Rees & Mészáros, 2005; Pe’er et al., 2006; Giannios, 2006; Beloborodov, 2010; Ryde et al., 2010, 2011; Vurm et al., 2013), with examples analysed in Ahlgren et al. (2015); Vianello et al. (2018); Ahlgren et al. (2019).

However, the simplest scenario for photospheric emission is a non-dissipative outflow (NDP: non-dissipative photosphere, henceforth) in which the thermal energy content of the flow is released at the photosphere, unaltered by heating in the flow (Goodman, 1986; Paczyński, 1986; Rees & Mészáros, 1994; Ryde, 2004; Ghirlanda et al., 2013; Larsson et al., 2015). In order for such emission to be detectable the flow has to become transparent close to, or below, where the flow saturates to its final outflow Lorentz factor, $\Gamma$. Otherwise, adiabatic expansion will diminish the thermal component and very little emission will be released. For non-dissipative photospheres the expected spectra are known in detail. Most importantly, geometrical effects cause the spectra to differ from a Planck spectrum. In particular, in the case of the photosphere occurring in the coasting phase, i.e., above the saturation radius, $r_{\rm s}$, the spectrum is significantly different from a Planck function (Goodman, 1986; Beloborodov, 2011; Lundman et al., 2013). Its spectral shape is much broader, while the asymptotic, low-energy spectral slope is still very hard, with $\alpha\sim 0.4$.

The most common way to assess different emission models in GRBs is by studying the sub-peak spectral shape. This is characterised by the photon index, $\alpha$, of the Band function or, similarly, of the cut-off powerlaw function (Band et al., 1993). The values of $\alpha$ can vary significantly depending on the emission process that is taking place. In a optically-thick, thermal scenario, the Rayleigh-Jeans limit has $\alpha=1$, the Wien limit has $\alpha=2$, and the non-dissipative photosphere in the coasting phase has the expected $\alpha$ of 0.4. Non-thermal emission is always in the regime $\alpha$ $\textless$ 0. This limiting value can be reached by synchrotron emission in many specific cases, such as from electrons with a small pitch angle distribution (Lloyd & Petrosian, 2000), jitter radiation (Medvedev, 2000), with attenuation by scattering (Dermer & Böttcher, 2000), and thermal synchrotron emission (Petrosian, 1981). Likewise the limiting values of $\alpha=0$ can be reached by synchrotron self-Compton emission from monoenergetic electrons scattering off a self-absorbed seed photons field (Stern & Poutanen, 2004). More generally, instantaneous or fast synchrotron cooling can explain $\alpha$ $\textless$ -3/2 and slow synchrotron cooling is expected to occur at $\alpha$ $\textless$ -2/3 (e.g., Tavani, 1996).

However, before comparing predictions of physical emission models to the observed characteristics of bursts, such as the distibution of $\alpha$-values, limitations of the detector and the analysis methods must also be taken into account. Such limitations could cause spuriously deviation from the expected $\alpha$ (e.g., Preece et al., 1998; Lloyd & Petrosian, 2000). One such limitation is due to the limited band-width of the detector which prevents the full spectrum to be detected. Another limitation is due to the typically empirical models that are used to fit the data; the Band function and the cut-off powerlaw function (see, e.g., Yu et al., 2016). If these empirical models do not match the true spectral shape, such as its curvature, then the fitted parameters might not be readily interpreted.

There are two possible routes to address these limitations. One way is to fit physical models directly to the data, which eliminates the need for fits to empirical functions (e.g. Lloyd & Petrosian, 2000; Zhang & Yan, 2011; Ahlgren et al., 2015; Vianello et al., 2018; Burgess et al., 2018; Ahlgren et al., 2019; Oganesyan et al., 2019). However, such analysis is computationally very demanding and, moreover, with the present, limited understanding of GRBs it is not clear what models should be used. Alternatively, one can assume a physical model and study what the response of the detector would be to such a model. This can be done by first producing synthetic data for a particular physical model, by taking into account the limitations of the detector (energy range, photon detection characteristics, etc.). The synthetic data can then be fitted with empirical models, accounting for the limitations of the typically adopted analysis methods. Such a procedure would allow to identify which parameter space of the empirical fit-model that the physical model corresponds to (see e.g., Burgess et al., 2015). Observed bursts which have spectral properties that coincide with the determined parameter space, can therefore be claimed to be compatible with that particular model. The advantage of this strategy is that the spectral analyses in all previous GRB catalogues can be directly assessed, since these have been produced using empirical models.

In this paper, we follow the second strategy and study a very specific and restrictive model, namely the non-dissipative photospheric emission (NDP), by producing synthetic GBM observations and by investigating the results of fitted empirical functions.

2 Properties of non-dissipative photosphere emission

The simplest photospheric emission spectrum that can be expected from a GRB is created in an undisturbed outflow, without any significant energy dissipation in the flow. The reason for this is that any dissipation of the flow kinetic energy will energise the electron population. Below the photosphere, these electrons subsequently might be able to distort the thermal photon spectrum into broader and more complicated shapes (e.g., Pe’er et al., 2006; Giannios, 2006; Vurm et al., 2011; Ahlgren et al., 2015; Chhotray & Lazzati, 2018). Numerical simulations of a jet emerging from a progenitor star do indicate that dissipation is expected to some degree (e.g., Ito et al., 2013; De Colle et al., 2017). However, the scenario without dissipation is still of interest to study, since such spectra define the narrowest and thereby the most extreme spectral shapes expected in the photospheric scenario. Moreover, if such spectra are identified in real cases, a strong limit on the degree of dissipation can be put.

But even in the absence of dissipation, there are still factors that cause the spectrum to differ from the comoving and local original thermal shape. One such factor is the fact that the photons’ last scattering off electrons in the flow occur at significantly different radii (e.g., Pe’er, 2008; Beloborodov, 2011). This could lead to a broader spectrum, if the temperature of the flow varies with radius. Another factor is the angular distribution of the photons in the lab frame, that is affected by the radial expansion; their distribution becomes more and more anisotropic the closer they are to the photosphere (Beloborodov, 2010). Finally, differences in the magnitude of the Doppler boosts for emission at high latitudes will also cause a broadening (Abramowicz et al., 1991; Pe’er, 2008).

These factors cause the spectrum from a photosphere occurring above the saturation radius (where the Lorentz factor saturates) in a non-dissipative flow to be broader than a Planck spectrum (Goodman, 1986; Beloborodov, 2011; Lundman et al., 2013). There is no simple analytical expression for this spectrum, and it has to be derived numerically. However, the shape of its photon flux spectrum can be approximated by a powerlaw with a stretched exponential cut-off:

[TABLE]

where $N_{\rm E}$ is the photon flux (unit), $E_{\rm pivot}$ is the pivot energy and $E_{\rm c}$ is the cut-off energy. In Figure 1 the approximation given by Equations (1) is shown as the light blue dashed line, which is overlayed on the numerically calculated spectrum in Lundman et al. (2013) (dark blue, solid line; see Fig. 1 in Ryde et al. 2017). In the figure the black line is the Planck function, aligned to the same spectral peak. The approximation in Eq. (1) is useful since it can easily be implemented in spectral analysis tools, such as RMfit111http://fermi.gsfc.nasa.gov/ssc/data/analysis/rmfit/, XSPEC (Arnaud, 1996), and 3ML (Vianello et al., 2017).

If, on the other hand, the photosphere occurs well below $r_{\rm s}$, then these factors will have a smaller effect and the narrowest allowed spectrum will be emitted. An analytical equation for this spectrum is given by Eq. (2) in Ryde et al. (2017), but can again be approximated by a simple analytical function; a cut-off powerlaw:

[TABLE]

The approximation in Equation (2) is shown in Figure 1 as the black, dashed line which is overlayed on the analytical expression in Ryde et al. (2017) (red, solid line).

3 Simulation and analysis of spectra

All the analysis for this study has been carried out within the The Multi-Mission Maximum Likelihood 3ML framework (Vianello et al., 2017). As a first step, we have generated synthetic spectra observed by the Gamma-ray Burst Monitor (GBM) onboard the Fermi Gamma-ray Space Telescope. GBM consists of 12 Sodium-Iodide (NaI) with $E$ $\epsilon$ ($8-900$ keV) and two Bismuth-Germenate (BGO) with $E$ $\epsilon$ ($205$ keV - $40$ MeV) detectors, totaling to 14 detectors.

To create these observed spectra, the raw model should be folded by the detector response, which is taken from the instrument’s standard response files that are available for every observation. We have chosen GRB120711115 for this task which is a generic long burst with a non-thermal looking spectrum and a large $E_{\rm pk}$. The simulated data consists of two NaI (NaI 2 and NaI A) and one BGO (BGO 0) detector that have the smallest source viewing angles in this particular observation, following the procedure of the Fermi GBM time-integrated and time-resolved catalogues (Goldstein et al., 2012; Gruber et al., 2014; Yu et al., 2016). The source and background intervals are obtained from the GBM catalogue. The source time interval is chosen as 62 s to 106 s after the trigger time (the fluence time). The background post-source interval is between 146-191 s while the pre-source interval is from -45 to -4 s (pre-trigger). The background selections are modeled with a zeroth order polynomial in time which is determined by a likelihood test by fitting the total count rate first. Following this, the polynomial is fitted to all energy channels and integrated over time to estimate the count rate from the background in each channel and their respective errors. It should be noted that the burst choice for this analysis does not effect the results since the fitted spectra are simulated from a theoretical model and only the response files are made use of from the chosen data files.

The synthetic spectra are then fitted within the Bayesian inference framework. All spectra are modeled with a powerlaw with an exponential cut-off. This is an empirical model that has been extensively used in the community, e.g. in the Fermi GBM catalogues (Goldstein et al., 2012; Gruber et al., 2014; Yu et al., 2016; Yu et al., 2019). For posterior simulations we have picked informative priors that specify realistic parameter intervals that can be detected with GBM and are reasonable for capturing the shape of our seed function. The normalization ($K$) is assigned a log uniform prior with $K$ $\sim$ ($10^{-11}$, $10$) cm*-2* keV*-1* s*-1*. The low energy index ($\alpha$) and the cut-off energy ($E_{C}$) are assumed to have uniform priors of $\alpha$ $\sim$ ($-3$ , $2$) and $E_{C}$ $\sim$ ($1$ , $10000$) keV respectively. The posterior sampling was performed via the emcee (an affine invariant Markov Chain Monte Carlo (MCMC) ensemble sampler) implementation in 3ML (Goodman & Weare, 2010). The model viability has been tested through posterior predictive checks and the convergence of MCMC simulations have been checked for, for which the details can be found in Appendices A and B.

We have implemented both Maximum Likelihood Estimate (MLE) and Bayesian methods in the simulation process to be able to double check. We find perfect agreement between the frequentist Maximum Likelihood (MLE) and Bayesian framework results.

4 Empirical characterisation of NDP photon spectra

Typically GRB spectra are described by empirical models. If these empirical models do not capture the true curvature of the theoretical emission spectrum, then the empirical model parameters might attain spurious values. This is particularly the case when the spectra are studied over a limited energy range. This has, for instance, been shown for synchrotron emission spectra (Lloyd & Petrosian, 2000). They showed that by using an empirical model, such as the Band function, to describe the synchrotron emission spectrum, the determined low-energy photon-index, $\alpha$, will depend on the position of the spectral peak, relative to the energy of the low-energy detector limit. Indeed, the fitted values of $E_{\rm pk}$ and $\alpha$ will give a well-defined correlation, even though the underlying spectral shape is unaltered.

On the other hand, if the empirical model used is able to properly capture the curvature of the theoretical emission spectrum, then the determined parameters will be valid, even outside of a limited energy range. For instance, Burgess et al. (2015) showed that if the true spectrum is a Band function, then the correct parameters are retrieved, independent of the position of the spectral peak. Likewise, Ryde et al. (2019) showed that this is also the case for a photospheric spectrum arising during the acceleration phase, when it is fitted by a cut-off powerlaw function. The determined $\alpha$-value is independent of the position of the peak relative to the low-energy detector limit. In such cases the curvature imprinted in the spectrum will be enough to find the correct parameter values, even though its asymptotic value has not been reached within the observed energy range.

However, since the true, theoretical emission model is not yet known, such spurious effects must be taken into consideration when interpreting data and their correlations. One way to do this is to analyse fits to the analytical (theoretical) spectra using empirical functions (§4.1) and another way is to analyse fits to synthetic spectra (§4.2), that are based on theoretical emission models. Below, we perform such analyses on the spectrum from a non-dissipative photosphere.

4.1 Curvature of a non-dissipative photosphere spectrum

We will now investigate how well the empirical functions can mathematically characterise the theoretical spectra from a photosphere occurring during the coasting phase in a non-dissipative flow, e.g., given by Eq. (1).

We compared Eq. (1) to a cut-off powerlaw function, over different energy ranges. First, we do this by fitting a cut-off powerlaw function, without limiting the energy range. The value of the fitted powerlaw index $\alpha$ is shown as the black line in Figure 2. The value is consistently $\alpha=0.32$. We note that this value is slightly different from the asymptotic powerlaw slope which has a photon index of $0.4$ (Beloborodov, 2010), giving a first indication that the cut-off powerlaw has a limited ability to perfectly reproduce the spectral shape.

Second, we limit the spectral range, over which the fit is performed, to be above $E_{\rm low}=8$ keV. This corresponds to the low-energy threshold of the GBM detector (Meegan et al., 2009). We find that the determined $\alpha$-value will depend on the position of the peak energy, $E_{\rm pk}$. The dark blue line in Figure 2 shows this relation between $\alpha$ and $E_{\rm pk}$. The value of $\alpha$ is consistently above 0 for $E_{\rm pk}$ larger than 100 keV.

A similar study was then done using the Band function as the empirical model. The corresponding result is shown in the red curve in Fig. 2. The results are very similar to the fits with the cut-off powerlaw function, indicating that the following investigation is valid for both functions.

We then tried a range of different low-energy limits, $E_{\rm low}$. A motivation for this is that the effective area of GBM does not reach its maximal value until 20 keV and therefore photon detections above 20 keV are more significant than below. The resulting correlations are shown by the series of light blue lines, using the range $E_{\rm low}=[10,12,14,16,18,20]$ keV. In the case of $E_{\rm low}=20$ keV, we note that a range of $\alpha>-0.6$ can be reached.

In summary, by just narrowing the analysed energy range, a positive correlation between fitted $E_{\rm pk}$ and $\alpha$ is expected, as an energy-window bias effect. Significant deviation occurs for $E_{\rm pk}$ below 100 keV. In particular, for spectra with $E_{\rm pk}$ larger than a few hundred keV, $\alpha$ is still very hard, being positive.

4.2 Synthetic Fermi/GBM spectra from a photosphere

While we above studied the functional spectral shape, we will now turn to simulating GBM data from the theoretical models. This produces synthetic photon data that can be studied with the typical data analysis methods and tools. These data will include all the effects of the dispersion of the instrument and the effective area variations according to the determined detector responses, which is not included in §4.1.

We start by producing a series of synthetic spectra which all have a similar signal-to-noise ratio, SNR $\sim 20$, but with different peak energies. We choose peak energies in the range 40–2000 keV, equally spaced in logarithmic scale. This is the range over which $E_{\rm pk}$ is typically observed (Yu et al., 2016).

These spectra are then fitted with a cut-off powerlaw. The results of the fits are shown as blue data points in Figure 3. In the figure we also plot the relations for $E_{\rm low}=8$ keV and $20$ keV from §4.1 for comparison. Two features can be noted. First, a positive correlation between the fitted $\alpha$-value and $E_{\rm pk}$-values, is in accordance to the window effect studied above, in §4.1. Second, there is an obvious offset between the mathematical fits over a limited energy band (light-blue lines) and the $\alpha$-values found from folding the spectra through the detector response (blue data points). Instead of the expected value of $\alpha\sim 0.3$ at large $E_{\rm pk}$-values, the synthetic spectra have fitted values of $\alpha\sim-0.1$.

We then produce a new series of synthetic spectra but now with a SNR $\sim 300$, in order to investigate whether the data significance affects the result. The $\alpha$-values from these fits are depicted by the green data points in Figure 3. These values are only slightly more negative and the conclusions drawn are the same as for the case with SNR $\sim 20$.

In order to investigate the reason for the off-set found in Figure 3, we plot in the upper panel in Figure 4 both the theoretical input spectrum, from which the synthetic data were produced (purple line), and the best-fit to these data by a cutoff powerlaw function (blue line). In this particular simulation $E_{\rm pk}=500$ keV. It is clear from this figure that the best-fit value of $\alpha$ is much softer than the actual spectral slope below the peak, elucidating the off-set seen in Figure 3. However, the fit captures the spectral shape relatively well above 50 keV. In the lower panel in Figure 4 the residuals of the ratio of the input spectrum (data) and the best-fit spectrum (model) is shown. The fit leaves residuals with a particular, wavy structure. The residuals are the smallest (within 10 %) between 35–2000 keV.

In summary, the cutoff powerlaw function does not properly describe the curvature of the spectrum from a non-dissipative photosphere. Therefore, as shown in Figure 3, the determined values of $\alpha$, from such a photosphere observed by GBM, lie around $\alpha\sim-0.1$, significantly differing from the asymptotic powerlaw slope of 0.4. The window effect (positive correlation between $E_{\rm pk}$ and $\alpha$) is mainly affecting the determined values of $\alpha$ below 100 keV. Therefore, one should expect $-0.4<\alpha<0$ from a non-dissipative photosphere observed by GBM.

5 Catalogue simulations

We will now redo the analysis done in the previous section, but we will instead randomly sample the $E_{\rm pk}$-values from the observed distributions from the GBM catalogue. This will give us the expected distribution of $\alpha$, in the case all observed spectra were to stem from non-dissipative photospheres.222The energy peak of the coasting NDP corresponding to peak flux $E_{\rm pk}$ values determined with the cut-off powerlaw function is first sampled randomly from a lognormal distribution with 8 $\leq$ $E_{\rm pk}$ $\leq$ 2000 keV. For each of the 100 spectra that are created, the normalization of the function is adjusted so that SNR $\approx$ 30 (also $\approx$ 300 to examine the high SNR condition).

We will sample from five different distributions. First, we will sample from the full $E_{\rm pk}$-distribution in Yu et al. (2016) (§5.1). Second, we will sample from the $E_{\rm pk}$-distributions belonging to the separate burst clusters that are defined in Acuner & Ryde (2018) (§5.2). We will focus on their clusters 1, 3, and 5, which all are characterized by having hard spectra. Finally, we will sample from the $E_{\rm pk}$-distribution associated to the hardest spectra ($\alpha_{\rm max}$) in Yu et al. (2016) (§5.3). The graphical results will be presented as Kernel Density Estimation (KDE) curves estimated with normal kernels.333See https://docs.scipy.org/doc/scipy/reference/generated

/scipy.stats.gaussian_kde.html

5.1 Synthetic $\alpha$-distribution by sampling from the GBM catalogue $E_{\rm pk}$-values

We sample from the $E_{\rm pk}$-distribution found from fits with the cut-off powerlaw function, made in the Fermi/GBM catalogue. The selected sample consists of 2228 bursts which have been detected by GBM until 19 January 2019. The median of the distribution is $\sim$ 260 keV, with minimum and maximum values of approximately 10 keV and 10 MeV, respectively.

We have randomly sampled from the above mentioned distribution to assign energy peaks to 100 NDP spectra simulations. These are then fitted with the cut-off powerlaw function. The resulting low-energy index ($\alpha$) distribution is shown by the blue line in Figure 5. The median value of the distribution is -0.16 and the inter-quantile range (IQR) is 0.15. This distribution should be compared to the asymptotic value of the theoretical NDP spectrum which is $\alpha=0.4$ (blue dashed line in Fig. 5). There is a significant offset between the two $\alpha$ estimations.

Next, we redo the analysis by creating synthetic data from a photosphere occurring when the flow is still undergoing acceleration (Eq. 2). This spectrum is much narrower and it can be approximated by an exponential cut-off powerlaw function. As mentioned above, a consequence of this is that the use of the cutoff powerlaw function in the fitting procedure will not give rise to a window bias effect. Thus, an artificial correlation between $\alpha$ and $E_{\rm pk}$ is not expected. The result from using the distribution of $E_{\rm pk}$ from the full GBM catalogue is shown by the red curve in Figure 5. The median value of the distribution is $\alpha_{med}=0.65$.

As a comparison the distribution of $\alpha$-values from the fits to the peakflux spectra in the GBM catalogue is shown by the grey line in Figure 5. For this distribution, we have applied the restrictions that the measured $\alpha$-values should be smaller than 3 and the $\alpha$-errors should be smaller than 1. This was done in order to remove obvious, and most likely, erroneous fits.

It should be noted that the peakflux spectra are produced by integrating the emission over approximately 1 second around the flux peak (e.g., Gruber et al., 2014). This means that they are not necessarily time-resolved, in the sense that they could still contain significant variability (see, e.g., Golkhou & Butler (2014), Golkhou et al. (2015), and Acuner & Ryde (2018)). Bursts with significant spectral evolution, on a time scale much shorter than the integration time scale for these spectra, will thus not reveal the instantaneous emission spectrum. As a consequence, the measured $\alpha$-values from the peakflux timebin might be significantly softer than that given by the emission process. In any case, the comparison of the $\alpha$-distribution for the coasting NDP spectra (blue curve) with the catalogue distribution (grey curve) shows that around 12% to 28% (see Table 1) of the catalogue peak flux spectra can be attributed to a non-dissipative photosphere in the coasting phase. This is determined by identifying the percentage of bursts in the GBM distribution that is inside the FWHM of the NDP $\alpha$-distribution for the former and inside the minimum and maximum values for the latter. The fraction of bursts that have $\alpha$-values larger than the minimum value is 38%. In contrast, only a small fraction of bursts from the GBM samples can be explained by photospheric spectra occurring during the acceleration phase (red curve).

5.2 Synthetic $\alpha$-distribution by sampling $E_{\rm pk}$-values from burst clusters

It is often of value to carry out an Exploratory Data Analysis (EDA) to detect any structure that might be innate in the data without assuming any models. This procedure can also reveal unexpected aspects of the data that have not been captured by existing models. Acuner & Ryde (2018) have carried out such an analysis applying an unsupervised clustering method to the Fermi GBM data which included the Band fit parameters from the GBM catalogue ($\alpha$, $\beta$, $E_{\rm pk}$) alongside the fluence and $T_{90}$. They identify five clusters and argue that 1/3 of the bursts could be explained by synchrotron emission (their clusters 2 and 4), whereas 2/3 are explained by a photospheric origin (clusters 1, 3, and 5). The latter class of bursts, which have hard $\alpha$-values, is of particular interest for this study.

We therefore redid the simulations and analysis made in §5.1, but now using the $E_{\rm pk}$ distributions from the five individual clusters from Acuner & Ryde (2018), instead. The basic statistical properties of all the simulated clusters are given in Table 2 and the distributions of clusters 1, 3, and 5 (photosphere clusters in Acuner & Ryde, 2018) are shown in Figure 6. The main difference from the full catalogue (simulated) distribution (§5.1; blue line in Fig. 5) is that the IQR are significantly smaller for clusters 1, 3, and 5. The main source of dispersion in the full catalogue (simulated) distribution of $\alpha$ is from the bursts in clusters 2 and 4. The only significant shift in median value is for cluster 5, which has a harder $\alpha$-distribution. This is due to the input $E_{\rm pk}$ distribution of the cluster 5 simulations. Compared to the other clusters, cluster 5 has the $E_{\rm pk}$ distribution with the highest $E_{\rm pk}$ mean (see Acuner & Ryde, 2018).

These distributions can now be compared to the distributions in the original, measured GBM clusters in Acuner & Ryde (2018), for which the percentages are presented in their Table 1. Furthermore, in Appendix C, we present all five simulated distributions with a comparison to the original data distributions from Acuner & Ryde (2018). The fraction of bursts in the five clusters that can be attributed to a NDP in the coasting phase are also shown in Table 1.

This comparison indicates that it is mainly clusters 1, 3 and 5 that show some overlap. The overlap for clusters 1 and 5 are still similar to that of the full catalogue, though. However, cluster 3 has the largest overlap; only less than half of the bursts have $\alpha$-values that are softer than expected from a NDP. A discussion of a possible cause for this is given in §6.2.

On the other hand, clusters 2 and 4 have very small overlap. This is not surprising since these two clusters are the ones with softest low-energy indices with median values of $\alpha$ are -0.77 and -1.43, respectively (Acuner & Ryde, 2018). This result is consistent with the arguments made in Acuner & Ryde (2018) that these bursts could be due to synchrotron emission. However, dissipation below the photosphere could also produce broad spectra. Therefore, physical model comparison should be done in order to assess the underlying physics for these two clusters.

5.3 Comparison to the maximal, time-resolved $\alpha$-value in each burst

In the previous sections we have compared the predicted range of $\alpha$-values from NDP to the observed distributions of $\alpha$ from the GBM catalogue. As mentioned above, these spectra are from the timebin at the flux peak for each burst, integrating over $\sim 1$ s (Gruber et al., 2014). Hence, this selection yields spectra that are not necessarily time-resolved enough, since the light curve variability is not considered. The five clusters discussed in § 5.2 all have different variability properties (Table 2 in Acuner & Ryde, 2018). In particular, we note that the median variability time for cluster 3 and 4 are around 1 second while for the other clusters it is $<0.5$ seconds. Therefore, the analysed timebins in clusters 3 and 4 might be time-resolved enough to exhibit the instantaneous emission spectrum.

On the other hand, time-resolved spectral catalogues have been presented, but they are for significantly smaller samples (e.g, Kaneko et al., 2006; Yu et al., 2016; Yu et al., 2019). Furthermore, the parameter distribution in these catalogues contain values from all the time-resolved spectra. As a consequence, these distributions contain varying number of spectra from each individual burst. This will give rise to a bias towards long and strong bursts, which have many time bins and therefore will dominate the $\alpha$-distribution. In addition, since spectral evolution is expected during individual bursts, the distributions also reflect the evolution in $\alpha$, even though there is no spectral evolution within timebins. The latter point is important since there are bursts which initially are unquestionably photospheric, but later evolve into a broader spectrum that can be fitted by non-thermal emission models (see, e.g. Ryde et al., 2011; Guiriec et al., 2015; Yu et al., 2019; Zhang et al., 2018; Ryde et al., 2019). Indeed, since emission processes typically are easiest distinguishable by the hardest $\alpha$-value that they allow, and if one assumes a single emission mechanism throughout a bursts, then $\alpha_{\rm max}$ is the most indicative parameter for the burst emission. In order to constrain the emission process during bursts, it should, therefore, be more informative to study the distribution of the maximal value of $\alpha$ in every individual burst instead.

To illustrate this, we plot in Figure 7 the distribution for $\alpha_{\rm max}$ for the 81 bursts in the GBM time-resolved catalogue (Yu et al., 2016). The time resolved distribution has an $\alpha$ median of -0.84, whereas the median for $\alpha_{\rm max}$ is -0.53. As can be seen from the figure, the distribution of $\alpha$ is significantly shifted to harder values. In particular, a majority of the spectra are now beyond the line-of-death for synchrotron emission, $\alpha=-2/3$ (green dashed line; Preece et al., 1998)). The fraction changes from 29% to 65%. The $\alpha_{\rm max}$-distribution in Figure 7 is similar to the corresponding distribution for the sample of 38 single pulses in Yu et al. (2019); see their Fig. 3, which reaffirms the result. The significant change in the $\alpha$-distribution thus calls for a separate investigation of the $\alpha_{\rm max}$-timebins.

In the upper panel in Figure 8, we plot the distribution of $\alpha_{\rm max}$ and their corresponding $E_{\rm pk}$-values, from the 81 bursts in the GBM time-resolved catalogue (Yu et al., 2016). The expected values for non-dissipative photospheres is shown by the fits to synthetic bursts (blue and green lines from Figure 3.) We find that around 26% of the observed points (21/81) are consistent with the synthetic burst spectra, since they cover the same region in the $\alpha$–$E_{\rm pk}$–plane: They lie above or within $1\sigma$ below the green line. In the lower panel in Figure 8, the corresponding relation is plotted for the 38 pulses in Yu et al. (2019). Here, we find that around 29% of the observed points (11/38) are consistent.

While some overlap exists between the two samples (Yu et al. (2016) and Yu et al. (2019)), we note that they are different in how the bursts were selected. Bursts in Yu et al. (2019) are selected to be single pulses, while Yu et al. (2016) also allowed multi-pulse, and complex bursts. Furthermore, Yu et al. (2019) uses the Bayesian Blocks method (Scargle et al., 2013) for creating the timebins, to which they later apply an significance threshold (Vianello, 2018), whereas Yu et al. (2016) only utilizes a signal-to-noise ratio threshold to determine the timebins.

Therefore, we argue that the clearest signature of the emission mechanism will be provided by the analysis of the Yu et al. (2019)-sample, depicted in Figure 8. This analysis indicates that around 29% of the bursts have at least one timebin with a spectrum that is consistent with a non-dissipative photosphere.

6 Discussion

We have studied the observed properties of photospheric emission in the specific case where the photon distributions are unaltered by energy dissipation in the flow, the so called non-dissipative photospheres. If such spectra are observed by a $\gamma$-ray instrument with a limited energy-band and are fitted by the typically employed empirical functions, then these spectra will yield a distribution of $\alpha$-values, lying in the range $-0.4<\alpha<0.0$. This is in contrast to the generally expected value of $\alpha=0.4$ for non-dissipative photospheres.

6.1 Explaining the observed $\alpha$-distributions

A consequence of this is that a non-negligible fraction of the $\alpha$-distribution from bursts observed by, for instance, the GBM, are compatible with being from non-dissipative photospheres. Earlier, only a handful of bursts have been claimed to be from NDP (Ryde, 2004; Ghirlanda et al., 2013; Larsson et al., 2015). Ryde et al. (2017) analysed two of these bursts, which exhibit extremely narrow spectra, namely GRB100507 and GRB101219. The spectral evolution were indeed consistent with a photosphere transitioning from the acceleration to the coasting phase of the flow, a possibility also discussed in Chhotray & Lazzati (2018). Here, our results indicate that many more bursts could be interpreted as stemming from such photospheres, just based on their $\alpha$-values (Fig. 5).

This is remarkable, since non-dissipative photospheric emission is only expected in very special circumstances, such as very smooth flows. The reason is that during the coasting phase the jet kinetic energy dominates the radiation energy. Therefore, dissipating only a small fraction of the kinetic energy can easily yield an energy density that is comparable to the energy density of the photon field prior to the dissipation. In such a situation the spectrum from the photosphere is expected to be modified from the shape in Eq. 1 (e.g., Rees & Mészáros, 2005; Pe’er et al., 2006; Giannios, 2006; Vurm et al., 2011; Ahlgren et al., 2015). Above the spectral peak, the spectrum can become harder due to Comptonization, while below the peak, low-energy, synchrotron photons, that are produced at the dissipation site, can make the spectrum softer. The change in the spectrum below the peak will cause smaller $\alpha$-values. The absence of such broadening thus places strong constraints on the amount of dissipation and photon production that is allowed.

It is thus a natural expectation that a part of the kinetic energy of the flow below the photosphere is allowed to dissipate in many bursts. A consequence of this is that an even larger fraction of the $\alpha$-value distribution could be consistent with emission from the photosphere. For instance, bursts that have $\alpha_{\rm max}<-0.4$ (limiting value for NDP) will then be compatible with a photospheric scenario as well.

This could also explain the distributions shown in Appendix C, in particular the limited overlaps between the observed and simulated $\alpha$-values (see Table 1). The fraction of bursts overlapping are from bursts not affected by dissipation, while the others are affected. The characteristics of the observed photospheric emission, would then depend on the individual properties of the flow yielding a varying amount of dissipation.

On the other hand, the complementary fraction of spectra in the $\alpha$-distribution ($\alpha<-0.4$) could be due to a wholly different emission process, such as synchrotron emission (see, e.g., Beniamini & Piran, 2013; Acuner & Ryde, 2018). In such a case one could imagine that during an initial phase of a burst, the emission from the photosphere dominates upon which synchrotron emission supervenes from a different part of the flow, either from internal or external shocks (e.g., Ghirlanda et al., 2004; Zhang et al., 2018; Li, 2018). However, the evolution of spectral parameters are typically smooth, which questions a two-zone emission scenario in which the two zones supplant in dominance of the emission (e.g., Ryde et al., 2019).

6.2 High-energy spectral shape and bursts in cluster 3

Low-energy indices are widely used to interpret radiation mechanisms from GRBs because of their distinguishability between various models. In this study, we concentrated on the $\alpha$-values which determine the seed radiative process that defines the lower energy tail of the spectra.

On the other hand, significant information also lies in the spectral shape above the peak energy, which therefore should not be neglected whilst inferring the emission physics. For instance, a high-energy power-laws can be the result of thermal Comptonization (Rybicki & Lightman, 1986), inverse Compton scattering (Stern & Poutanen, 2004), or thermal synchrotron emission (e.g, Petrosian, 1981).

However, spectral analysis in large samples of GRBs suggest that the best empirical model is the cut-off powerlaw function and that the data do not require a high-energy power-law (Goldstein et al., 2012; Gruber et al., 2014; Yu et al., 2016). In particular, Yu et al. (2019) showed that if the analysed timebins are chosen to avoid spectral evolution, and therefore can be considered to exhibit the instantaneous emission spectrum, then the Band function $\beta$ becomes softer (more negative) and thus more compatible with an exponential cut-off. In such cases, the $\alpha$-value of the spectra is the decisive parameter and the conclusions on the fraction of bursts consistent with NDP spectra based on the $\alpha$-distribution are valid.

On the other hand, a particular feature of the spectra in one of the clusters in Acuner & Ryde (2018) (cluster 3) are very hard values of the high-energy, powerlaw index $\beta\sim-2$. The presence of such a powerlaw indicates that energy dissipation has occurred and that the Compton $y$-parameter is $y\sim 1$. The high-energy, power-law emission could either be due to optically-thin Compton scattering, which would give rise to a wavy high-energy spectrum (as indicated in Acuner & Ryde (2018)) or optically thick Comptonisation, producing a pure power-law. In either case, such a feature also indicates that the flow has dissipated an amount of energy which is approximately of the same amount as the thermal photon energy.

Interestingly, we found above (§5.2) that, comparing clusters 1, 3, and 5 (photosphere clusters in Acuner & Ryde, 2018), cluster 3 differs in that it has the largest overlap between the measured and the expected NDP $\alpha$-values (see Table 1 and Appendix C). However, even though a large fraction of the bursts in cluster 3 have $\alpha$-values that correspond to the expectation of NDP spectra, they are obviously not all non-dissipative, due to the presence of a high energy tail in the spectra.

This could therefore, in part, explain the fact that cluster 3 has a larger overlap between the measured and the expected NDP $\alpha$-values, compared to clusters 1 and 5. The fraction of spectra affected by dissipation might be similar between the three clusters, however, the effects of the dissipation could be different. In clusters 1 and 5, the dissipation mainly affects the low-energy part ($\alpha$) of the spectra: The dissipation produces low-energy (synchrotron) photons which are energised and populate the spectrum below the peak. This results in that a majority of all spectra are not compatible with the NDP. In contrast, the effect of the dissipation in the cluster 3-bursts could be different. For instance, for some of the burst, the low-energy photon production could be suppressed, leaving the $\alpha$-slope is unchanged. At the same time, the energy dissipation occurring in these bursts would mainly cause changes in the high-energy part of the spectrum. Thus, the actual fraction of bursts in cluster 3 without dissipation, i.e., truely non-dissipative, should be lower than what is given in Table 1.

As a conclusion, only focusing on $\alpha$-values, therefore, does not necessarily reveal the actual presence of dissipation. An investigation of the low-energy index should therefore be followed by an investigation of the whole spectrum, in order to address the effects of Comptonisation. This requires different physical models to be fit to data and to be compared through a statistically reasonable methods.

6.3 Interpretation of $\alpha$ and extrapolation to the X-ray spectrum

Figure 4 shows that the fitting process of minimizing the fit statistics of the cut-off powerlaw function, results in a fit that only correctly determines the spectral curvature above $\sim 35$ keV. However, the fit is not able to correctly reproduce the actual powerlaw slope at low energies. The $\alpha$ parameter attains a value which optimizes the match to the actual spectral curvature around the peak and therefore it does not necessarily fit the spectral shape below $\sim 35$ keV.

This can in part be understood by the significance of the counts per channel for the synthetic photosphere spectra. The significance of the data will be the largest at a few 100 keV. This is partly due to the decrease in effective area towards lower energies and partly due to the subpeak powerlaw slope of the photon spectrum (Eq. 1) having a slope steeper than $N_{\rm E}\propto E^{0}$, that is, an increasing photon number with energy. For typical spectra the significance can vary by a factor of four. Therefore, more weight in the fit is given to the count data around a few 100 keV, even baring in mind that dispersion of incoming photons occur in the detector to produce the count spectrum (Meegan et al., 2009). As a consequence the data far below the peak will have less influence on the over all fit. Hence, the fitted $\alpha$-value does not necessarily represent the asymptotic slope of the incoming spectrum.

In such a case, $\alpha$ should be interpreted as a parameter that is used to fit the spectral curvature around the peak, and not the asymptotic powerlaw slope of the incoming spectrum, which is typically done. Hence, the determined value of $\alpha$ should, for instance, not be used to extrapolate the GBM data to lower energies, such as, the X-ray or optical bands (Ghirlanda et al., 2007; Sakamoto et al., 2009). The residuals from the spectrum in Figure 4 can become large. For instance, at 1 keV, their ratio is 0.2. This means that if (i) the true spectrum is a NDP, and (ii) if an extrapolation were to be made into the X-ray energy range and below, based on the $\alpha$-value alone, the spectral flux will be overestimated by nearly an order of magnitude. Note that the parameter $\alpha$ is not necessarily badly determined statistically. It is the interpretation of it that might be wrong, simply due to the fitted model not being the accurate model, which, in this case, is assumed to be the NDP.

6.4 Determination of $E_{pk}$

In the catalogue simulations made above, we have sampled the $E_{\rm pk}$-values that we use for the NDP spectra, from the $E_{\rm pk}$-distribution found from the cut-off powerlaw fits. This assumes that the $E_{\rm pk}$ is well determined by the cut-off powerlaw. On the other hand, the fitted $\alpha$-values can be quite off from the expected value. This is shown in Figure 5, where the determined $\alpha$-values are shown by the blue line, which is significantly off-set from the asymptotic power law slope ($\alpha=0.4$) of the NDP spectra, which is shown by the blue, dashed line. However, Figure 9 shows how well the $E_{\rm pk}$ is determined by showing the ratio of $E_{\rm pk}$ and $E_{\rm pk,input}$ versus $E_{\rm pk,input}$. The ratio is consistently close to 1.0, and there is no strong trend. This, therefore, shows that there is no disadvantage of sampling from the fitted $E_{\rm pk}$ from using the cut-off powerlaw function.

7 Conclusion

We have shown that GRB spectra that have $\alpha>-0.4$ have low-energy slopes that are compatible with a photospheric emission from a non-dissipative outflow. This corresponds to a non-negligible fraction of all GRB spectra (Fig. 5). The reason one should not expect a hard $\alpha$-value from photospheric emission, such as $\alpha=1$ (Rayleigh-Jeans slope) or $\alpha=0.4$ (including geometrical effects) is that the empirical models used in spectral analysis (cut-off powerlaw and Band) do not match the curvature of the incoming spectrum. A consequence of this is that the determined value of $\alpha$ will be consistently softer than the asymptotic slope of the theoretical spectrum ($\alpha\sim 0.4$). The fitted value is instead closer to $\alpha\sim-0.1$.

We find that more than a quarter of all bursts have at least one time-resolved spectrum, whose $\alpha$-values are consistent with a non-dissipative outflow photosphere. This is the most restrictive and simplest photospheric scenario. However, since it is natural to expect a degree of dissipation of kinetic energy in the flow below the photosphere, the fraction of spectra consistent with emission from the photosphere should be even higher.

The conclusion is that (i) a non-negligible fraction of the observed $\alpha$ distribution can be explained by a non-dissipative photospheric model. (ii) It is not advisable to extrapolate the determined value of $\alpha$ to lower energy bands, such as the X-ray and the optical bands. (iii) The interpretation of $\alpha$-values to assess emission models can be misleading and must be done with caution.

Acknowledgements

We thank Dr. Lundman for enlightening discussions. We acknowledge support from the Swedish National Space Agency and the Swedish Research Council (Vetenskapsrådet). FR is supported by the Göran Gustafsson Foundation for Research in Natural Sciences and Medicine.

Appendix A Bayesian formalism

Bayesian inference applies the Bayes’ theorem to infer and update parameter estimations, probabilities and distributions regarding a model after the experimental data is obtained. Bayes’ theorem is given as

[TABLE]

where $p(\theta\mid y)$ is the posterior distribution of the estimated parameter $\theta$ given the data $y$.

A model’s ability of describing the data at hand can be put to the test by assessing the predictions it makes about the process that has created the observed data (y). If $y_{rep}$ is the replicated data set from the model’s fit to the data, then the distribution of $y_{rep}$ conditional on y is called the $posterior$ $predictive$ $distribution$ and is shown as (Gelman et al., 2014a),

[TABLE]

A model is capable of describing the data gathered from a certain process if the future unknown observables ($\tilde{y}$) are successfully described by $y_{rep}$. In other words, the observed data $y$ and the replicated data $y_{rep}$ should be coming from the same generative process. In this paper, checks of posterior predictive distribution have been used to assess the quality of the fits done (see Appendix B).

Appendix B Checks of model fit viability

B.1 Corner plots and convergence diagnostics

For all fits the corner plots displaying the posterior distributions of the resulting parameters have been visually checked. Posteriors which show clear bi or multi-modalities have been refit with refined initializations of the parameters to get unimodal normal-like posterior distributions for all simulations.

The convergence of the MCMC chains have also been assessed through trace plots to make sure that the true stationary distribution has been reached.

B.2 Posterior Predictive Checks and QQ plots

The viability of the spectral fits for the various simulations in this work has been assessed via posterior predictive checks (PPC). This includes the replication of the posterior distributions from the fitted spectra which are then compared to the observed data. A mismatch of the observed and replicated data would then indicate in what ways the model fails to represent the data, pointing to possible aspects of the model that should be modified or extended.

A compact way of comparing the observed and replicated data is via plotting the quantile-quantile (QQ) plots. Any significant deviation from the one to one line would then indicate the potential shortcomings of the current model to describe the data at hand (Gelman et al., 2014b).

Following the method of Burgess et al. (2018), we have generated 500 realizations from each posterior given by the spectral fits. We then simulated the counts form these spectra and compared them to the observed data counts in QQ plots with 95% and 68 % quantiles superimposed. A deviation from 95% quantile for a significant portion of the plot would indicate an unviable fit. We find no significant deviations in the presented simulation fits. We also compare the QQ plots of the fits to the cutoff powerlaw function to that of the real model used for the simulations. The former plot should not have great deviations from the latter for the fit to be acceptable. We find that this criterion is also satisfied in our fits.

Appendix C Original and simulated clusters comparisons