The Fraction of Gamma-ray Bursts with an Observed Photospheric Emission Episode

There is no complete description of the emission physics during the prompt phase in gamma-ray bursts. Spectral analyses, however, indicate that many spectra are narrower than what is expected for non-thermal emission models. Here, we reanalyse the sample of 37 bursts in \citet{Yu2019}, by fitting the narrowest time-resolved spectrum in each burst. We perform model comparison between a photospheric and a synchrotron emission model based on Bayesian evidence. We choose to compare the shape of the narrowest expected spectra: emission from the photosphere in a non-dissipative flow and slow-cooled synchrotron emission from a narrow electron distribution. We find that the photospheric spectral shape is preferred by $54 \pm 8 \%$ of the spectra (20/37), while $38 \pm 8 \%$ of the spectra (14/37) prefer the synchrotron spectral shape; three spectra are inconclusive. We hence conclude that GRB spectra are indeed very narrow and that more than half of the bursts have a photospheric emission episode. We also find that a third of all analysed spectra, not only prefer, but are also compatible with a non-dissipative photosphere, confirming previous similar findings. Furthermore, we notice that the spectra, that prefer the photospheric model, all have a low-energy power-law indices $\alpha>-0.5$. This means that $\alpha$ is a good estimator of which model is preferred by the data. Finally, we argue that the spectra which statistically prefer the synchrotron model, could equally well be caused by subphotospheric dissipation. If that is the case, photospheric emission during the early, prompt phase would be even more dominant.


INTRODUCTION
Gamma-ray bursts (GRBs) are bright flashes of gamma-rays emitted from a highly relativistic jet that is ejected during the formation of a black hole in a transient collapse (see, e.g., Mészáros 2019). Within the fireball model of GRBs fluctuations in the jet lead to internal dissipation, which provides energy that can be radiated away. The observed emission depends on the strength and localisation of such dissipation. However, the smallest radius from which observable emission can originate is from the photosphere, at r = r ph , where the jet becomes optically thin. The strength of the photospheric emission depends on the relative amount of energy that is dissipated below the photosphere and on the adiabatic energy losses. Numerical models of subphotopsheric dissipation show that intense and broad spectra can be produced (e.g., Rees & Mészáros 2005;Pe'er et al. 2005;Giannios 2006;Beloborodov 2011).

Sample
Our choice of sample is motivated by our desire to study the intrinsic emission in the early phase of gamma-ray bursts and get unambiguous attributes of the emission process. This task is, however, not straight forward since the GRB emission typically evolves during its duration, with the intensity and spectral peak energy changing (e.g., Golenetskii et al. 1983;Kargatis et al. 1994), but equally important, also changing in the spectral shape and width (e.g., Wheaton et al. 1973;Crider et al. 1997). These facts have two consequences. First, to observe and study the intrinsic emission, it is vital to study as narrow time interval of the light curve as possible, in order to avoid the observed spectrum to be smeared by spectral evolution. At the same time, though, the signal-to-noise ratio needs to be sufficiently high. Second, a large fraction of the timebins will have spectra that are so broad that, as mentioned in the introduction, they are ambiguous since many emission models can produce them.
Therefore, the best chance of identifying the intrinsic spectrum is to restrict the investigation to individual pulses in the light curve, which naturally avoids overlapping emission episodes and rapidly varying spectra. Moreover, in order to assess the compatibility of the data with different emission models, the narrowest (intrinsic) spectrum in each burst is of particular interest.
We consequently choose to study the 37 pulses in the catalogue of Yu et al. (2019). These pulses were identified as individually connected emission episodes, which can offer sufficiently narrow time-intervals, while maintaining a high statistical significance. This selection ensures that the spectra represent the intrinsic emission mechanism, by avoiding integration of temporally varying emission. Yu et al. (2019) showed that these time-resolved spectra are all well fit with an exponentially cut-off power law function, with a power law index, α, which is defined as where N E and K CPL is the photon flux and normalisation (cm −2 s −1 keV −1 ), E c and E piv = 100 keV is the cutoff energy and pivot energy. We examine the spectral fits to all of the time intervals in Yu et al. (2019) in order to identify the time intervals which have the largest value of α, one from each pulse. This selection thus creates a sample of 37 spectra, which all are the narrowest spectra in each individual pulse, and which, as far as possible, represent the intrinsic emission. This sample is presented in Table 1. For the analysis, we used the standard Fermi/GBM analysis proceedure (e.g. Goldstein et al. 2012;Yu et al. 2016). This includes, for instance, using an energy range of 8 keV to ∼ 850 keV for the NaI detectors, and ∼ 250 keV to 40 MeV for the BGO detectors, and using three NaI detectors and one BGO detector, with viewing angles of less than 60 degrees from the location of the burst.

Emission Models
Both synchrotron emission and emission from the photosphere could be expected during the prompt emission phase in GRBs (e.g., Mészáros 2019). These emission processes have in common that they can produce a variety of spectra, with varying widths. However, specific to each model is the theoretical limits to how narrow the spectra are permitted to be. By comparing the narrowest observed spectra with these two limiting spectral shapes, one can compare the models abilities to explain these narrow spectra. The primary purpose of this study is therefore to make a model comparison between the narrowest expected photospheric and synchrotron spectra.
The narrowest synchrotron emission spectrum is expected from mono-energetic electrons, which are in a given Bfield, and that are in the slow cooling regime. However, since any plausible synchrotron-emitting source does not have mono-energetic electrons, we choose to use the emission from a set of electrons with a distribution in energies instead. The theoretical prediction for Fermi acceleration of the electrons provides a power-law index of the electron distribution of either p = 2 or in the range of p = 2.2-2.4 in the relativistic case (Sironi et al. 2015). Likewise, the GRB afterglow emission is well described by synchrotron emission from electrons with p ∼ 2 (Wijers & Galama 1999). We know, though, that the theory of particle acceleration is incomplete and we cannot exclude other values at this points. We therefore choose a larger value of p = 3.8, which makes the synchrotron spectrum narrower. This value is motivated mainly by the fact that it is the typical value found when slow-cooling synchtrotron spectra are fitted the data (e.g., Tavani 1996;Burgess et al. 2019a;Ronchi et al. 2019).
Moreover, a plausible synchrotron-emitting source is not either expected to have a single value of B in the emitting region. However, for simplicity, we still retain a single value of B, which is an approximation that avoids an additional broadening of the spectrum and that is typically used in most spectral fits of synchrotron emission. Finally, we assume an isotropic distribution of the electrons' pitch angle, which is the angle between the electron momenta and the magnetic field. In the specific case of an anisotropic pitch angle distribution, the small-pitch-angle synchrotron emission allows for harder spectral slopes below the peak (Epstein & Petrosian 1973;Lloyd & Petrosian 1999;Medvedev 2000;Lloyd-Ronning & Petrosian 2002). However, it is not clear that such emission is applicable to GRBs (Baring & Braby 2004;Yang & Zhang 2018). In the analysis below we use the emissivity function of synchrotron radiation in a random magnetic field for a power-law distribution of electron energies (Aharonian et al. 2010), implemented in the Python package Naima (Zabalza 2015). The emission spectrum used here is shown by the green curve in Figure 1. A power-law with photon index of −2/3 is shown by the blue line, which corresponds to the asymptotic power-law slope for synchrotron emission. This model is referred to as the slow cooling synchrotron model (SCS) below.
Likewise, for the photospheric emission the narrowest spectrum is expected in the absence of any significant energy dissipation in the jet that would alter the photon distribution 1 . The spectral shape is given by numerical simulations Normalized photon energy, E Normalized SCS NDP α = -2 / 3 Figure 1. Theoretical spectral shapes that are used to fit to the data, plotted as νFν -spectra. Slow cooled synchrotron emission with electron index p = 3.8 (green curve) and non-dissipative photosphere emission (red curve). The grey power-law has a power law photon index α = −2/3, corresponding to the asymptotic slope of the synchrotron spectrum. (Beloborodov 2010;Bégué et al. 2013;Ito et al. 2013;Lundman et al. 2013) and shown by the red curve in Figure  1 (see also, Fig. 1 in Ryde et al. (2019)). An analytical approximation to the spectrum is given in Equation (1) in Acuner et al. (2019). This model is called the non-dissipative photosphere model (NDP) below.
We point out that our primary aim is not to determine the values of the models' parameters, but rather to compare the statistical strength of support for the spectral shapes. However, the peak energy and the normalisations of the models can be translated into physical parameters (Tavani 1996;Pe'er et al. 2015). The ranges of these parameters can serve as additional indications to the validity of the models by comparing them to expected values from theory (Beniamini & Piran 2013;Ghisellini et al. 2019). However, to do this, further assumptions are needed about the values of other unknown parameters, such as the bulk Lorentz factor, dissipation efficiencies, and the redshift.
Although the primary goal is to compare the narrowest expected physical emission spectra, we expand the study by including fits with the the empirical cutoff power-law (CPL; Eq. 1). CPL is not a physical model, however, it has a flexibility to vary the low-energy power-law slope, α, which indicates the shape of the spectrum below the peak. As such the CPL spectrum can characterize expansions of the two emission models, when the strong assumptions that were made for the two limiting spectral shapes presented above are relaxed. At the same time, the CPL is a reasonable approximation to some particular physical models such as fast-cooled synchrotron emission (see, e.g., Burgess et al. 2015) and emission from a photosphere in a photon-dominated flow (Ryde et al. 2017;Acuner et al. 2019).

Bayesian Model Comparison
We now have a sample of the narrowest spectrum in each burst which will be fitted to the narrowest allowed theoretical spectra from the two models. The fits to these two models can then be compared statistically. In this work, we choose the Bayesian approach for model comparison.
Bayes' theorem gives the probability of the model (M ) being true given a set of observed data (D) as where P (M ) is the prior probability of the model, P (D) is effectively a normalising factor that ensures the posterior probabilities of the different models sum to unity and P (D | M ) is the marginal likelihood or Bayesian evidence (Jaynes & Bretthorst 2003). The marginal likelihood of the nth model, denoted Z n = P (D|M n ) for brevity, is given by where θ n are the parameters of the nth model, P (θ n |M n ) is the prior distribution of the models parameters, and P (D|θ n , M n ) is the likelihood under the nth model given specific parameter values. In other words, the Bayesian evidence Z is the average of the prior-weighted likelihood of a model over the parameter space. The calculation of the evidence thus requires an integration over the whole parameter space and is computationally demanding, which is one reason why it is not commonly used. However, recent computational developments make this numerically accessible (Skilling 2004;Feroz & Hobson 2008). The Bayesian evidence Z is central to the task of model comparison. When Bayesian evidences are obtained for two competing models, the ratio of the respective evidences, Z 2 /Z 1 , summarizes the evidence given by the data in favor of model M 1 or model M 2 , or equivalently in natural logarithm as the Bayes factor Using this notation, the ratio of the posterior model probability between two models is given as, When there is no a priori reason to believe that one model is more probable than the other, our prior ratio for the two models becomes unity and the remaining relation at the right hand side equation gives the Bayes' factor (see, e.g., Jaynes & Bretthorst 2003). Finally, the posterior model probability when the prior ratio is unity is which gives the probability as a sigmoid function of the Bayes factor (eq. 4). As a guide, a difference in log-evidence greater than 2 will have a model probability that corresponds to a belief of larger than 90% in the model compared to the other model tested. Between the two models that are compared, we therefore state that a model is preferred if the log-evidence difference |ln(Z 2 /Z 1 )| 2. The exact value of this threshold is subjective, but similar criteria are typically employed (Jeffreys 1961(Jeffreys , 1957Kass & Raftery 1995).
In the analysis performed below, we calculate the marginal likelihood or the Bayesian evidence Z (Eq. 3) for each spectrum fitted and report the values of ln Z n for the competing models, their difference in log-evidence; ln Z 2 /Z 1 (Eq. 4), as well as, the posterior model probability P (Eq. 6).
If this model comparison shows that the non-dissipative photospheric model is statistically preferred over the slowcooled synchrotron model, then the conclusion can be drawn that the observed spectrum is so narrow that synchrotron model (slow cooled, fast cooled, or marginally fast cooled) is not the best description of the data. On the other hand, if such a model comparison shows that the slow-cooled synchrotron model is statistically preferred over the nondissipative photospheric model, then the spectrum is broad enough to be compatible with any synchrotron emission model. However, it cannot be concluded that the emission is not from a photosphere, since the possibility of dissipation below the photosphere can produce a broad spectrum (Pe'er et al. 2005;Beloborodov 2011;Ahlgren et al. 2015). To be able to draw any conclusion in such a case, a different model comparison must be performed, between a full synchrotron model and a full subphotospheric dissipation model.

Analysis and Numerical Method
To sample from the posterior distributions and evaluate evidence integrals, we use Multinest which has been shown to successfully and efficiently sample from complicated multimodal posteriors while calculating the evidences in its default run (Feroz & Hobson 2008;Feroz et al. 2011). Multinest utilizes the nested sampling (NS) algorithm (Skilling 2004). It enables efficient evaluation of the Bayesian evidence where the posterior distribution is given as a by product of the performed run. Therefore, NS requires no additional computations to estimate the marginal likelihood which would be the case for other MC algorithms. Below, we describe in a simplistic way how the algorithm functions.
Multinest improves over the standard NS algorithms by utilizing an ellipsoidal rejection sampling scheme which provides better geometrical flexibility for posterior exploration and the ability to identify distinct modes. Due to the nature of the NS method, Multinest runs are generally successful in converging on the highest likelihood region. However in some cases, this region might be left outside of the selected elliptical regions during the initial estimations and due to the gradual shrinkage of the initial region the true high likelihood space might be missed, causing misleading posterior distributions. For this reason, we use corner plots to assess the validity of the posteriors. If the posterior distribution does not reflect the expected results in any way, the fits are run multiple times to ensure convergence as well as the live points given to the algorithm are increased so that there is a greater chance to sample the highest likelihood region without eliminating it.
As the evidence Z, given in Eq.
(3), is in general computationally expensive to generate, we also calculate the Akaike Information Criterion (AIC) in order to test whether it can be used on this problem. The AIC is mathematically simpler measure that could be used as an approximation to evaluate the strength of support for a model. The AIC is defined as where L max = P (D | θ max , M ) is the maximum likelihood, θ max is the maximum likelihood parameter values, and k is the number of parameters of the model. Then, the below relation can be written for NDP and SCS models 2 , where θ NDP and θ SCS are the estimated maximum likelihood parameters from the NDP and SCS model fits to the data. We further utilize posterior predictive checks (PPC) to check for the viability of the model fits. If a model describes the data well, then the replicated data from the fit should look similar to the observed data that is used to for modelling. Any possible systematic differences in these two distributions hint to the aspects in which the preferred model cannot describe the observed data accurately. It should be noted that although PPCs can and do reveal major discrepancies, the related model can still be used for some purposes. In our case, we showcase the model fits with poor PPC results as a further tool for understanding how models can fail to be preferred by data both through evidences and the incapability of reproducing the observational quantities tied to the data (Gelman 2004). We use the two physical models for model comparison without the assumption of any one of them being the standalone processes that might have single-handedly generated our data. Hence, the output of the Bayesian model comparison (performed in §3.1) is not to determine if a certain model is "correct" but to see its potential failings and expand the model to encompass these aspects of the data for the future analyses (preformed in §3.2). In this regard, PPCs help to pinpoint the exact source of any problematic issues present in the analysis (Gabry et al. 2017). In our analysis, we do not omit the fits with bad PPC results from the presented results. However, this indicates the preferred model is not fully capable of describing the data and hence needs to be modified and/or expanded. This is further discussed in Section 3.2.

Selection of Parameter Priors for the Analysis
Priors for each of the two emission models (NDP, SCS) as well as for the empirical CPL model are estimated via the empirical Bayes method. In this procedure, which is a practical step towards more advanced hierarchical modelling, we first derive the approximate priors from the catalogue values and physical arguments. These are used as priors in the Multinest fits. The posterior distributions of all the Multinest fits for each model then are turned into priors for our actual runs which are the shrunk versions of the initial universal priors assigned. The priors are calibrated this way so that we can assure that we have good sampling and convergence with the given number of live points in a reasonable computational time. By fixing the parameter priors for each model to be the same throughout the sample, we further assure that proper evidences are obtained, which can be affected if different prior intervals are used for different bursts for the same model. We also compare the prior intervals estimated by Multinest to those estimated by emcee (MCMC) to make sure that we have properly explored the parameter spaces for all models.
Most priors are distributed as log uniform. The priors for the non-dissipative photosphere spectrum are chosen to be where E NDP is the break energy and K NDP is its normalisation, and U(x min , x max ) is a uniform and log U(x min , x max ) is a log-uniform distribution of the variable x, between x min and x max . For the slow-cooled synchrotron spectrum, we use the correpsonding priors Finally, the empirical cutoff powerlaw model, which is defined in Equation (1), has the following priors: We note that the prior ranges differ somewhat, simply due to the differences in spectral shape, and their definitions. For instance, the parameters E NDP , E SCS and E CPL are the break energies of the respective functions and hence are not identical with each other nor with the peak energy, E pk .

Sensitivity to Priors
The AIC value (Eq. 7) evaluates the best fit likelihood and is, unlike the evidence (Eq. 3), insensitive to the parameter priors assuming the true likelihood maximum is within the parameter range encompassed by the priors. Therefore, we will compare these two measures to see if they differ, in order to assess whether it is reasonable to approximate the full evidence calculation with the AIC. The AIC values are shown in Table A.2 in the Appendix.
In Figure 2, we compare the corresponding differences in log-evidence and AIC for the fits to the NDP and the SCS models. A positive difference in log-evidence indicates that the NDP model is favoured, while the opposite is true for the AIC. The figure shows that the AIC differences are strongly correlated with the difference in log-evidences. It is therefore evident that, in this case, the AIC measure captures everything important in the fit. The correlation also indicates that our results are not critically sensitive to the prior choices, since the AIC is insensitive to the parameter priors. We conclude that, at least for the priors we have adopted and with the data and model set described, the AIC difference is a reasonable approximation to the difference in log-evidences.
The robustness of the priors have also been tested by assigning different sized prior ranges to the fits and checking the fits and evidence values. We use the smallest possible prior ranges that include the greatest likelihood region in order to facilitate the numerical estimation of the Bayesian evidence. We also make an effort of approximately matching the prior ranges across parameters for the three different models used (see §2.5).

RESULTS
In order to compare the two models (NDP and SCS), we calculate the Bayesian evidence (Z, or marginal likelihood, Eq. 3) for each of the models, assuming the prior distributions specified in § 2.5. In Table 1, we show the results of the fits to the 37 spectra in our sample. The GRB number ID are followed by the log-evidence values (ln Z) of the two models, as well as their differences, ln(Z NDP /Z SCS ) (Eq. 4), and corresponding posterior model probability, P (Eq. 6). For ln(Z NDP /Z SCS ) > 0, P is the probability for the NDP model to be the better model, while for ln(Z NDP /Z SCS ) < 0, P is probability for the SCS model to be the better model. Consequently, P is always larger than 0.5. In Table A.2, in the Appendix, we include more information of the analysed spectra.

Model Preference
We will now discuss the differences in log-evidence, ln(Z NDP /Z SCS ), shown in Table 1 Table 1. CPL low-energy index (α), log-evidences for NDP and SCS, the Bayes factor ln(ZNDP/ZSCS), as well as the posterior model probability P from eq. (6) are presented. For ln(ZNDP/ZSCS) > 0, P is the posterior model probability for the NDP model, while for ln(ZNDP/ZSCS) < 0, P is the posterior model probability for the SCS model. Consequently, P is always larger than 0.5.
have | ln(Z NDP /Z SCS )| 2 and that therefore are excluded in the discussion 3 , since they are indecisive. Table 1 shows consequently that for 20/37 (54%) of the analysed spectra, the photosphere model is preferred, since ln(Z NDP /Z SCS ) > 2. Similarly, for 14/37 (38%) of the spectra, the synchrotron model is preferred, since ln(Z NDP /Z SCS ) < −2.
We note that for most of the bursts the difference is log-evidence, ln(Z NDP /Z SCS ), is large. This means that we can claim that one of the models is preferred over the other with great confidence, but not that either model is necessarily a good fit to the data. Table 1 shows that the bursts with the two largest difference in evidences are GRB100707, overwhelmingly in favour of the photosphere model, and GRB110721A, overwhelmingly in favour of the synchrotron model. Apart from these burst, there are four bursts with ln(Z NDP /Z SCS ) > 100 and one with ln(Z NDP /Z SCS ) < −100 (see Table 1). In the Appendix, we present some of the analysis products, such as the PPC plots and the spectral fits, for two bursts (GRB150314 and GRB150902) to illustrate the analysis performed and what good and bad fits typically look like. We proceed to investigate whether the model preference depends on the α-value found from fits to the traditionally used empirical model, the cut-off power law function (CPL, Band et al. 1993). This spectrum is characterized by the low-energy power law index, α (see Eq. 1). Therefore, in our analysis we have also refit the spectra to the CPL function. For our spectra, such fits are already presented in Yu et al. (2019) and Ryde et al. (2019), who used MCMC for posterior sampling. Here, instead, we use Multinest as a sampler ( §2.4) in order to be able to be able to calculate the CPL fit evidences (used in §3.2). We find similar fit results to within errors for all the analysed spectra.
In Figure 3, we plot the difference in log-evidence between NDP and SCS versus the corresponding value of α. In the plot the line of α = −2/3 is indicated, which corresponds to the asymptotic power-law photon index for the SCS model. Figure 3 shows that there is a strong correlation between ln(Z NDP /Z SCS ) and α. When the NDP is the preferred model, then the α-value typically is −0.5, while for the cases where α is −0.6, then the SCS is systematically the preferred model, and finally the three indecisive cases lie between −0.7 α −0.5. This means that the α-parameter can be used as an approximate determination of the model preference of the data. Eluding to the 'line-of-death' of synchrotron emission (Preece et al. 1998), we can thus define a 'line of non-compatability' of synchrotron emission, lying at α = −0.5, instead.

Expanding Beyond the Two Models
In the analysis above we have compared the strength of support for two specific models, NDP and SCS and determined which model is the preferred one. However, the actual property of the burst emission could be different from any of the tested models. synchrotron emission (see, e.g., Tavani 1996;Lloyd & Petrosian 2000;Burgess et al. 2019a). A way to investigate this, is to compare the fits of the physical models to fits to the empirical cutoff powerlaw (CPL) model, eq. (1). The flexibility of this function is given by the power-law slope, α, which can capture all of the model expansions mentioned above. The log-evidence values found from fitting our spectra with the CPL model are given in Table A.2. We now focus on the 22 burst spectra which have a preference for the NDP over SCS. The bursts are shown in Figure 3 above the dashed line since they have ln(Z NDP /Z SCS ) > 0. In Figure 4 we plot the difference in log-evidence between the NDP model and the CPL model as a function of the CPL parameter α. The figure shows that there are six bursts with positive difference in ln Z NDP − ln Z CPL (Table A.2). This means that the physical NDP model is preferred over the empirical CPL spectrum, even though the latter has more flexibility by having one more free fit parameter. Most of these cases occur in an interval around α −0.2. Moreover, Figure 4 shows that for α −0.3 the CPL is the preferred model for most cases. The same is for the spectra with α 0. We interpret this as the spectra are not fully correctly represented by the NDP, in particular, the low-energy slope is not correct.
A natural explanation for this is that the α −0.3-cases are better produced by subphotospheric dissipation, which naturally produces a soft low-energy shoulder of the spectrum (see, e.g. Pe'er et al. 2005;Ahlgren et al. 2019). Of the spectra that have α ≤ −0.3 there are five cases for which the two requirements hold, namely ln(Z NDP /Z SCS ) > 2 as well as ln(Z NDP /Z CPL ) < −2. These five spectra can thus be claimed to be due to subphotospheric dissipation 4 .
Likewise, it can also be argued that the cases with α 0 are better fit by a spectrum that is even narrower than the NDP (see, e.g., Ryde et al. 2017;Acuner et al. 2019). The only possibility for this to occur within the fireball model is if the photosphere takes place in the acceleration phase, before the flow saturates (Beloborodov 2011;Ryde et al. 2017;Chhotray & Lazzati 2018). Since dissipation is not expected below the saturation radius, such spectra are naturally non-dissipative. Therefore, together with the spectra that prefer NDP (with positive difference in ln Z NDP − ln Z CPL ) we therefore can argue that all the 13 bursts with α > −0.3 are consistent with being non-dissipative. We note that the fraction of non-dissipative spectra (13/37 = 35%) that we find here is similar, but somewhat larger, than the fraction of spectra consistent with non-dissipative photospheres found by Acuner et al. (2019). They used a different approach by analysing synthetic spectral data, which were produced by simulating data from the theoretical NDP spectrum and then convolved through the GBM response. Acuner et al. (2019) showed that when such synthetic data, having a range of spectral peak energies, are fitted with a CPL function it results in −0.  Yu et al. (2016). However, it is worth pointing out that these two different methods, one based on comparing simulated spectral shapes to the catalogue shapes and the other based on model comparison, give similar results. This consolidated the conclusion that a non-negligible fraction of GRB spectra are consistent with photospheric emission in a non-dissipative flow.

NDP preferred CPL preferred
In summary, our statistical analysis and consideration concludes that, not only 20/37 (54%) of the analysed spectra are compatible with a photospheric origin, but also that 13/37 (35% ; Fig. 4) are actually compatible with photospheric emission in a non-dissipative jet.

DISCUSSION
We are addressing the fundamental mechanism of the emission during the prompt phase in GRBs. In order to do this we study the narrowest spectra in 37, well-detected, GBM bursts. We discriminate between the radiation processes by comparing the spectral shapes of a non-dissipative photosphere (NDP) spectrum and a slow-cooled synchrotron (SCS) spectrum.
In the choice between NDP and SCS, we find that 54 ± 8 % of the spectra prefer the photospheric model while 38 ± 8 % prefer the sychrotron model 5 . It is important to bear in mind that the statistical answer we arrive at is the spectra's preference between these two models only. This means that we did not, necessarily, arrive at the best model for the spectrum, which could be a different model. Indeed, we showed, in §3.2, that spectra with −0.6 α −0.3 still do prefer the NDP model over the SCS model, but there must be some additional broadening mechanism to explain the spectra fully. Therefore, we argued that a subphotospheric dissipation model is an even better description of the data, compared to NDP and SCS, for these cases.

Subphotospheric Dissipation when α −0.6
The fact that a broadening mechanism is required for some of the spectra with α −0.6 leads to the natural consequence that a fraction of the rest of the sample (spectra preferring the SCS over NDP) also could be from the photosphere. The reason is that if a photospheric broadening mechanism is needed for some bursts, there is no reason to believe that it is limited to α −0.6. On the contrary, values down to α ∼ −1.1 are theoretically expected (Thompson & Gill 2014;Vurm & Beloborodov 2016). There are only three spectra 6 in Table 1 that, given the uncertainties, are consistent with having α < −1.1. Since these three spectra are not naturally formed by subphotospheric dissipation, it is difficult to argue for anything but prompt synchrotron emission for them. This deliberation points to a natural possibility that a significant majority of the analysed spectra are indeed photospheric and only a small fraction are due to synchrotron spectra: 35% (13/37) are from non-dissipative photospheres, 57% (21/37) are from dissipative photospheres, and 8% (3/37) are from synchrotron emission.

Similar Spectra: Synchrotron and Photosphere
The deliberation in the previous section argued that many of the spectra preferring the SCS over NDP actually could be photospheric. Indeed, the spectral shape of synchrotron emission and subphotopsheric dissipation spectrum can be similar in shape. To illustrate this point, we will now create synthetic spectra from a specific scenario within the subphotospheric dissipation model and then fit them with a synchrotron model. Synchrotron spectra are, by nature, very broad, emitting over a large energy range. However, spectra formed by subphotospheric dissipation can also be broad. They are, though, inherently marked by a low-energy break, relating to the seed photons for the unsaturated Comptonization up to the peak energy. The position of the low-energy break depends on the exact prescription of dissipation, but can be several decades in energy below the main peak (Pe'er et al. . Further, such a spectrum will be marked by a high-energy cutoff at around ∼ 50 − 100 MeV due to intrinsic pair opacity (Pe'er et al. 2005;Beloborodov 2011). In any case, over the limited energy range of the GBM detector the subphotospheric spectrum can resemble a Band spectrum, or even a synchrotron spectrum.
Our synthetic spectra are produced using the model of Pe'er et al. (2005). One such spectrum is shown by the (synthetic) data points in Figure 5. The data are generated from a model (green curve), with a small amount of subphotospheric dissipation in a high luminosity jet. The model parameters are typical for other observed bursts that are fitted by the model (Ahlgren et al. 2015;Ahlgren et al. 2019). The data are generated using XSPEC (Arnaud et al. 1999) and have a signal-to-noise ratio of ∼ 20. As shown by the figure, the fit to the synchrotron model (grey curve) gives a good representation of the data. The AIC value for the synchrotron fit is 273.2, while the AIC value for a fit to the photospheric model itself gives 273.4, which are not significantly different. The data at high energies are in this case not enough to reliably discriminate between the two models, as is often the case. However, data at energies 10 keV could in this instance help distinguish between the models. By inspection of this fit alone it is not possible to tell that the data have been generated from a scenario of subphotospheric dissipation.

Compatible Model versus Preferred Model
Turning back to the analysed spectra, which were the hardest spectra in each burst, we found that 38% of the spectra strongly prefer a synchrotron spectrum over the NDP. Indeed, many works have identified GRB spectra that are compatible with synchrotron emission (e.g., Tavani 1996;Lloyd & Petrosian 2000;Burgess et al. 2011;Oganesyan et al. 2017;Burgess et al. 2019a).
Again, from a statistical perspective, all these previous findings do not prove the synchrotron model, but they show that the data is compatible with a certain spectral model. The spectral fit in Figure 5 illustrates this inherent ambiguity of the GBM data to clearly distinguish between the two physical models. A further example is GRB170114A which 6 GRBs 100122, 101126, 130614. Figure 5. Synthetic data from a subphotospheric dissipation model (green curve) for photospheric emission. The grey curve shows the synchrotron spectral fit in the GBM energy interval, which gives a good fit. The fit is performed in XSPEC excludes data in the region of 30 − 40 keV around the K-edge. Red data points are synthetic data from the NaI detector and blue data points are from the BGO detector.
provides a good fit to a slow-cooling synchrotron model (Burgess et al. 2019b). In our analysis above, we found a clear preference for NDP for this burst (Table 1; ln(Z NDP /Z SCS ) = 14). This means that there is a physical model that is statistically preferred by the data over the SCS model. The conclusion of this fact is that the emission in GRB170114A is not synchrotron but it is from the photosphere. On the other hand, including the cut-off powerlaw function (CPL) in the model comparison (see §3.2), the CPL becomes the preferred model for GRB170114A. This is mainly due to the fact that α ∼ −0.5, which is very soft for photospheric emission. The conclusion we draw from this analysis is that the best physical model for (largest α-timebin in) GRB170114A is a dissipative photospheric emission.

Afterglow Overlapping and Superseding the Prompt Emission
The spectra we analysed in this paper are from pulses in the light curve, which are mainly in the vicinity of the GBM trigger, within ∼ 10 s (see Table 1 and Yu et al. 2019). Only eight of the spectra occurred more than 10 s after the trigger. Since the afterglow onset timescale is on the order of a few tens of seconds, the afterglow component could be expected in observations much beyond this time scale.
This point is further underlined by the fact that many burst are observed to have a second emission component which is consistent with being synchrotron emission and that coexists with the prompt emission. In particular, observations of high-energy, LAT photons indicate that they are typically delayed relative to the GBM trigger, but emerge while the GBM emission is still active (e.g., Ackermann et al. 2010;Giuliani et al. 2014;Ajello et al. 2019b).
These components have been associated to the onset of the afterglow emission from an external blastwave, involving the reverse shock and forward shocks (e.g., Kumar & Barniol Duran 2009;Ghirlanda et al. 2010;Panaitescu 2017). Alternatively, such extra components have been successfully explained to be due to inverse Compton cooling of the heated e ± enriched plasma, created behind the forward shock in the blastwave, due to the exposure of the prompt keV-MeV emission Vurm et al. 2014) and that involves two phases, an inverse Compton phase and a synchrotron-self-Compton phase. Recently, a very clear example of such emission is given in GRB190114C Ajello et al. 2019a). Here a synchrotron component decidedly emerges during the prompt phase and this component is unmistakably connected to the synchrotron component of the early afterglow emission of the blastwave (Fraija et al. 2019;Derishev & Piran 2019). In the case of GRB190114C the synchrotron component is initially observed to fluctuate until it finally settles into the temporal power-law of the long term afterglow. The afterglow component therefore coexists with the prompt emission in many bursts.
In other cases, the observed photospheric emission is superseded by a synchrotron component. Examples are GRB150330A and GRB160625B, which have an initial, short photospheric episode, while their main emission episode occurs a few hundred seconds later, which are dominated by synchrotron emission Li 2019). Other similar cases have the components separated by a few 10s of seconds instead (GRB140329B and GRB160325A; Li 2019; Sharma et al. 2020). This type of bursts illustrate the fact that if the initial thermal pulse had been missed for some reason, they would have been detected with a purely synchrotron-like prompt emission.
It was, in fact, early argued (Rees & Mészáros 1992) that prompt emission might be absent and that all observed emission could be synchrotron emission from the external shock, i.e., the afterglow (see further discussion in Panaitescu & Mészáros 1998). An example of such a case is GRB141028A, where all the observed emission, both the "prompt" phase as well as the afterglow phases, are interpreted to be from one and the same process, which is consistent with an external shock ). In the analysis above, this particular burst has ln(Z NDP /Z SCS ) = −23 (Table 1), indicating a strong preference for the synchrotron model, thus supporting such an interpretation. A similar suggestion was made for the interpretation of GRB110721A (Iyyani et al. 2016), which had ln(Z NDP /Z SCS ) < −200 above.
The possibility of observing only the external shock emission is in line with the fact that when synchrotron spectra are found to be compatible with the prompt emission, a general finding from the fits is that the emission radius, the bulk Lorentz factor, and the minimum electron energy, have to be large, in order to avoid fast cooling of the electrons (Nakar & Piran 2017;Iyyani et al. 2016;Burgess et al. 2016Burgess et al. , 2018Oganesyan et al. 2019;Ravasio et al. 2019). In particular, the emission radii that are found ( 3 × 10 16 cm) are close to typical values for the deceleration radius 7 . These facts again raises the possibility that a detected "prompt" synchrotron emission could actually from the external shocks.
In order to remove the obvious ambiguity and unclarity in GRB emission nomenclature, given by the use of the "prompt" and "afterglow" emission phases, it might therefore be clearer to instead use the dichotomy between the photosphere (with various degrees of dissipation; Rees & Mészáros (2005)) and the blastwave (Synchrotron self-Compton emission; Rees & Mészáros (1992) or inverse Compton cooling; Beloborodov et al. (2014)).

Spectral Evolution During Pulses
One particularity in our analysis is that we only studied one spectrum from each burst, namely the narrowest spectrum. This was motivated by the fact that the spectral shape varies significantly within pulses (Wheaton et al. 1973;Crider et al. 1997;Ryde et al. 2019) and therefore the narrowest spectrum is the most constraining for any emission mechanism. Moreover, if the emission mechanism is the same throughout the pulse, then the narrowest spectrum is the most informative for the pulse emission as a whole. We note that these spectra typically occur close to the pulse peak , which ensures high signal rates.
One advantage of such a choice, compared to data samples that include many timebins for each of the studied bursts, is the avoidance of the larger fraction of broad spectra (α < −2/3), simply caused by the spectral broadening (Wheaton et al. 1973), and that therefore makes model determination ambiguous. A second advantage is that such a choice will not be biased towards bursts with many time bins, which complicates the interpretation of the parameter distributions .
Within photospheric emission models the observed broadening of the spectra during a pulse is typically taken as an indication of a varying strength of subphotospheric dissipation (e.g., Ryde et al. 2011;Ahlgren et al. 2015). Furthermore, the narrowest spectrum is expected as the photosphere occurs close to the saturation radius, since then dissipation is less likely and at the same time the photospheric emission is the brightest . Therefore, we argue that if a pulse has one spectrum that is convincingly photospheric, then the whole pulse is photospheric.

CONCLUSIONS
We have investigated the narrowness of GRB prompt emission spectra, in order to infer constraints on possible emission mechanisms. Our sample consists of one spectrum in each of 37 GRB pulses, all selected by having the maximal α-value during the pulse, i.e., the narrowest spectrum. The spectra all occur close to the pulse peak and mostly within 10 s of the trigger time. We have shown that more than half of these spectra statistically prefer the non-dissiptive photosphere spectrum (NDP) over the slow-cooled synchrotron spectrum, thereby indicating narrow spectra and a photospheric origin. Making the natural assumption that the whole pulse due to the same emission mechanism, this conclusion translates to the whole pulse, not only for the narrowest spectrum.
We also confirm earlier claims in Acuner et al. (2019) that a non-negligible fraction of all spectra (1/3), not only statistically prefer, but are also compatible with photospheric emission in a non-dissipative jet. This result is important for two reasons. First, it shows that a noticable fraction of GRB spectra are formed in a region which is not substantially affected by dissipation. This could be due to the photosphere occurring either close to the saturation radius or in a laminar region of the flow. Second, in the photospheric scenario only a small fraction of spectra are expected to be non-dissipative, since only little disturbances in the flow is needed to change the spectral shape (see, e.g. Pe'er et al. 2006;Ryde et al. 2011Ryde et al. , 2017. Therefore, identifying a non-negligible fraction of NDP spectra means that the main part of the remaining spectra should also be photospheric, albeit suffering alterations due to disspation. Indeed, among the spectra in our sample there are cases having relatively soft spectra (with −0.6 α −0.3) for which the data still prefer the NDP model. We have shown that there must be some additional broadening mechanism to explain these spectra, presumably subphotospheric dissipation. Arguing further, the fact that a broadening mechanism is required for these spectra leads to the natural consequence that a fraction of the rest of the sample, which do not prefer the NDP, also could be from the photosphere. Consequently, there is a natural possibility that a significant majority of the analysed spectra are simply photospheric and only a small fraction are from synchrotron emission: While a third are from non-dissipative photospheres, slightly more than half would thus be from dissipative photospheres, leaving only very few to be from synchrotron emission. Such a possibility, however, cannot be firmly established based on the spectral shape alone, as illustrated in Figure 5. Other indications, such as the high degree of polarisation that is observed during the prompt phase in some bursts, is a strong indicator for a non-thermal origin of these spectra (e.g., Fraija et al. 2017;Sharma et al. 2019). Moreover, any sample of the prompt emission is most likely contaminated with afterglow emission, which naturally provides synchrotron spectra from the external shock being observed during the "prompt" phase (e.g., Burgess et al. 2016).
For the model comparison, we have made use of the Bayesian evidences. We have, however, shown that, at least for the spectra in our sample and for the parameter priors we employ, the simper AIC measure gives comparable results. Furthermore, we have shown that the traditionally used "line-of-death" for synchrotron emission (Preece et al. 1998) can be used as a rough discriminator between photospheric and synchrotron spectra. However, the limit is not exact, and α −0.5 is safer to use.
Our conclusion is therefore that the initial emission episode in most GRBs is photospheric, be it with or without spectral broadening mechanisms, such as subphotospheric dissipation. However, the actual fraction of bursts with episodes dominated by synchrotron emission during the "prompt" phase is difficult to assess and still needs to be determined, using other observables, such as gamma-ray polarisation. The large fraction of photospheric spectra that we find in this paper, is in line with the current paradigm of GRB emission with an efficient photosphere combined with an external blastwave component (see review by Mészáros (2019)). However, the properties and cause of the subphotospheric dissipation are not firmly established. This can, though, be investigated by further analysis of observations and numerical simulations of different jet compositions and dynamics. By generating data from a selected model or directly from the observed data, we obtain the count spectrum which is the sum of burst flux convolved with the GBM response and the background. These are used to create the Detector Reponse Matrix (DRM). Since this matrix is not square, the equation linking counts to energy spectrum cannot be solved and hence the method of forward folding is required. The model flux vector is folded through the response matrix which gives the model count spectrum. Then, the model count spectrum is compared to the observed counts and a new model flux vector is calculated. Once this iteration leads to satisfactory results, the best fit parameters are obtained. 8 Since the GBM data needs the forward folding method to obtain resulting spectra, we compare the generated counts and asses in which ways they deviate from the actual spectra which presents a handy tool for determining good and bad fits from the actual analysis in this paper. As shown in Figures 6 to 9, an ideal fit gives a near perfect one to one relation in the PPC plots and a randomly distributed residual plot whereas a bad fit shows predicted cumulative counts that are outside the 68 % and 95 % quantiles for a significant part of the data which looks very different from the one to one relation. Furthermore, the spectrum of such a fit would show systematic differences in the residuals plot.
A.1. Analysis of the fits from synthetic spectra In this section we perform the full statistical analysis for sets of synthetic data, produced from the two theoretical models, NDP and SCS. In such a way we can have clear examples of what good (simulated data fitted by the model it was generated from) and bad fits (simulated data fitted with another model) look like in PPC plots. Additionally, this provides a check on the method used for the analysis, identifying possible coding errors and any faulty selection of prior ranges.
The PPC plots for the analysis on synthetic data are presented in Figs 6 and 8, which summarises the abovementioned results.

A.2. Further Results of the Analysis on the Sample
In Table A.2 we present more details of the analysis of the 37 spectra in our sample. For every spectrum we present the statistical significance, S, of the timebin (as defined in Vianello 2018), and its time-interval (cf. to Yu et al. (2019)). Furthermore, we present the statistical measures AIC (Akaike 1974), BIC (Kass & Wasserman 1995), and DIC (Gelman et al. 2014) for the nondissipative photosphere and the slow cooled synchrotron spectral fits.
All the fits were initially performed with a maximum likelihood estimation. We used these point estimates as initial starting values for our fits within the Bayesian framework. We find that both sets of fits give consistent estimations within the error bars for all bursts.  Figure 6. Posterior Predictive Check (PPC) count plots from a well fitted spectrum, with three GBM detectors shown in separate panels. The red line demonstrates the one-to-one relation. In the lower right-hand panel the count spectrum and the best fit model are shown, together with the residuals. The data are from a synthetic NPD spectrum with a SNR=100 and it was fitted with the NPD spectral model.   Figure 8. Same as Fig. 6, but for a badly fitted spectrum. The synthetic data (NDP spectrum, SNR=100) are now fitted with the SCS spectral model, instead. Deviations from the one-to-one lines are apparent.

A.3. Detailed Analysis Results of Two Example Spectra
As part of the analysis the fits were judged based on a number of characteristics. First, the spectral fit for the various models is displayed on the observed counts and shown together with the residuals. A wavy structure of the residuals indicate a bad fit.
Second, we present the PPC plots, which again gives an indication of how well the model fits the data. Third, the Bayesian analysis produces corner plots of the posterior distributions of the parameters. These fits are considered good if a Gaussian like distribution is obtained for the posteriors.
To illustrate this proceedure, we present the detailed analysis on GRB150314 (Fig.10) and GRB150902 (Fig. 12) as examples. We again show the same plots as in Figure 6, but now these figures are on the real data from GBM. For GRB150314 we see a perfect one to one relation for the PPC plots and a random distribution of the residuals in the spectral fit which both indicate a good fit. For GRB150902 we see PPC plots that do not match with the one to one PPC line and the spectrum residuals are showing systematic differences which indicate a bad fit.
For both of these fits we have Gaussian like posterior distributions which shows the fits have converged properly. However, the ability of the two different models to describe the data is very different as seen from the PPC plots and the residuals in the spectra.  Figure 10. Same as Fig. 6, but for a well fitted observed GBM spectrum. GRB150314205, with SNR=20, is given as an example. The data are fitted with a NDP model.   Figure 12. Same as Fig. 6, but for a badly fitted observed GBM spectrum. GRB150902A, with SNR=93, is given as an example. The data are fitted with a SCS model.