Testing a model for subphotospheric dissipation in GRBs: fits to Fermi data constrain the dissipation scenario

It has been suggested that the prompt emission in gamma-ray bursts (GRBs) could be described by radiation from the photosphere in a hot fireball. Such models must be tested by directly fitting them to data. In this work we use data from the Fermi Gamma-ray Space Telescope and consider a specific photospheric model, in which kinetic energy of a low-magnetisation outflow is dissipated locally by internal shocks below the photosphere. We construct a table model with a physically motivated parameter space and fit it to time-resolved spectra of the 36 brightest Fermi GRBs with known redshift. We find that about two thirds of the examined spectra cannot be described by the model, as it typically under-predicts the observed flux. However, since the sample is strongly biased towards bright GRBs, we argue that this fraction will be significantly lowered when considering the full population. From the successful fits we find that the model can reproduce the full range of spectral slopes present in the sample. For these cases we also find that the dissipation consistently occurs at a radius of $\sim 10^{12}$ cm and that only a few percent efficiency is required. Furthermore, we find a positive correlation between the fireball luminosity and the Lorentz factor. Such a correlation has been previously reported by independent methods. We conclude that if GRB spectra are due to photospheric emission, the dissipation cannot only be the specific scenario we consider here.


INTRODUCTION
Prompt GRB emission mechanisms are still a matter of debate, and remain an unsolved problem in GRB physics. The principal approach to finding the origin of GRB prompt emission has been to fit spectra with empirical functions and use the best-fitting parameters to constrain different emission scenarios. For a considerable time the empirical Band function (Band et al. 1993), a smoothly broken power law, has been the predominant way of fitting these spectra. The Band function often provides good fits to the prompt emission spectra and the fits have commonly been interpreted in terms of synchrotron emission (Tavani 1996;Briggs et al. 1999;Abdo et al. 2009;Zhang et al. 2016). However, it has been shown that synchrotron radiation cannot consistently produce the full distribution of spectral shapes we observe when fitting with the Band function, including the E-mail: bjornah@kth.se low-energy slope (Preece et al. 1998) and the spectral width of the peak (Axelsson & Borgonovo 2015;Yu et al. 2015). In addition, the Band function does not represent a physical scenario, and in order to obtain physical information we must instead ultimately use a physically motivated model and fit to data (see also Burgess et al. 2014aBurgess et al. , 2015. Observations of GRBs with spectra close to blackbodies (BBs) (Ryde 2004;Ghirlanda et al. 2013;Larsson et al. 2015), as well as spectra described by a BB and an additional component (Ryde et al. 2010;Guiriec et al. 2011;Axelsson et al. 2012;Burgess et al. 2014b;Arimoto et al. 2016), suggest that at least some bursts have prompt emission of a photospheric origin. Conversely, for most GRBs, the simplest models of photospheric emission cannot provide an adequate description of the observed spectra. In particular, the low-energy spectral slope of most bursts is significantly softer than a Planck function.
It has been shown that more realistic photospheric models will result in spectra that are broader than a Planck function. The broadening can be a result of, e.g., subphotospheric dissipation (Rees & Mészáros 2005;Pe'er et al. 2006;Giannios 2006;Beloborodov 2010;Vurm et al. 2011;Chhotray & Lazzati 2015), or geometric effects (Pe'er 2008;Lundman et al. 2013). Subphotospheric dissipation has previously been invoked to interpret the 'top-hat' spectral shape of GRB 110920 (Iyyani et al. 2015) and the broadening of the initially very narrow spectrum of GRB 090902B (Ryde et al. 2011). Though subphotospheric dissipation has been used to describe several bursts, the different model scenarios must be examined further, particularly using larger data samples and proper techniques when fitting models to data.
In our previous work (Ahlgren et al. 2015), A15 from here on, we fitted a model for subphotospheric dissipation to prompt GRB emission data from the Fermi Gamma-ray Space Telescope. This was the first time such a model was fitted to data and the study served as an important proofof-concept. Specifically, we fitted GRB 090618, which can also be described by a Band function with typical parameters, and GRB 100724B, which has previously been fitted with a composite Band+BB model. The results showed that subphotospheric dissipation is a viable description of both bursts and justified further investigations.
In this work we build on the work presented in A15 by increasing the model parameter space and by including approximate adiabatic cooling. We have also significantly expanded our sample, now analysing 36 bursts, and have performed Monte Carlo simulations to quantify some of the model's inherent uncertainties. We also consider parameter correlations and physical implications.
The main goal of the work is to test the scenario of localised dissipation by internal shocks below the photosphere in a low-magnetisation outflow. We emphasise that the constraints found in the paper do not apply to photospheric emission in general or other dissipation scenarios. We confine our study to the case where almost all the dissipated energy goes into the electrons and explore the effects of three free parameters, as described in section 2.2. The second main objective is to present the continued development of our table model, which we plan to make publicly available.
The paper is organised as follows: section 2 describes the physical scenario and the table model implementation.
In section 3 we describe our data sample, binning and fitting procedure. Section 4 presents the results and section 5 is the discussion. Section 6 is summary and conclusions. Finally, in appendices A-E we present more technical details about the model and analysis.

Physical scenario
In this paper we consider the same physical model as in A15. It is a model based on the fireball model, specifically a hot fireball, using a kinetic code by Pe'er & Waxman (2005). We described this model in A15, but for completeness we present an overview here as well. Any deviations from the model as stated in A15 will be explicitly pointed out. Note that all quantities in this section are in the observer frame unless otherwise stated.
In the fireball model (see e.g. Pe'er 2015 for a recent review), an isotropic equivalent luminosity, L 0 = L 0,52 10 52 erg s −1 , is emitted from the central engine in the form of baryons, electrons, B-fields, and photons. The outflow starts at the initial jet radius, r 0 , where the bulk Lorentz factor is Γ ∼ 1. The photons give rise to a radiative pressure which accelerates the outflow into highly relativistic collimated jets. The outflow is accelerated until the saturation radius, r s = Γr 0 . At this point the outflow is in thermal equilibrium and coasts with constant Γ. The baryons, mainly protons, carry the bulk of the kinetic energy of the outflow. In our model we define the dissipation radius as r d = r ph /τ = Lσ T /4πτΓ 3 c 3 m p , where τ is the optical depth due to the baryonic electrons (i.e. the electrons which are associated with the protons of the outflow). At this radius we dissipate a fraction ε d of the kinetic energy of the protons, out of which the fraction ε e goes to electrons and a fraction ε b is used to generate magnetic fields.
The dissipation creates a non-equilibrium situation between the photon bath and the heated electron population, where they interact through Compton and inverse Compton scattering, as well as pair production and pair annihilation. We also treat synchrotron radiation and synchrotron selfabsorption. It is assumed that electrons are injected as a Maxwellian distribution with a dimensionless temperature, θ = kT/m e c 2 ∝ ε d ε e , in the range of 5.5 -220. Both newly injected electrons and the existing population will interact with the photons. Before the dissipation essentially all energy is carried by the baryons in the jet. After the dissipation the dissipated energy has been transferred to the photons, but the baryons still carry the majority of the energy.
Possible dissipation mechanisms in GRB jets include magnetic reconnection (Thompson 1994;Drenkhahn & Spruit 2002;Giannios & Spruit 2005;Lyutikov 2006), hadronic collisions (Beloborodov 2010) or internal shocks (Kobayashi et al. 1997;Daigne & Mochkovitch 1998;Rees & Mészáros 2005). In our current implementation we are working with the internal shocks paradigm by assuming that r d = Γ 2 r 0 . We also assume that the dissipation is a continuous process localised between r d and 2r d . Note that this is different from the scenario with continuous subphotospheric dissipation throughout the jet employed in e.g. Vurm & Beloborodov (2016).
A caveat is that large a value of ε d would remove the majority of the kinetic energy of the outflow and would inevitably decelerate it. This in turn may lead to a reacceleration phase, which is not treated in our code. We therefore avoid the corresponding part of the parameter space. Since the energy is passed to either the electrons or the magnetic fields, this translates to ε d ε e + ε d ε b being the constrained quantity. We have decided to require ε d ≤ 0.4 to minimise this effect while still keeping a large amount of parameter space to work with. Even at this level of dissipation one could expect a second acceleration phase (Bégué & Iyyani 2014), but we expect this to have a limited impact on our results.
There are additional physical effects not accounted for in the code. Examples include spatial effects, such as the aforementioned geometric broadening, the jet opening angle, jet hydrodynamics (e.g. Lazzati et al. 2009) and the effects from a fuzzy photosphere (Beloborodov 2011;Bégué et al. 2013). Finally, we do not at this stage consider neutrino emission or afterglow predictions.

Table model
The XSPEC table model implementation (Arnaud et al. 1999) allows us to build a numerical model from simulated spectra. A table model consists of a grid of model spectra, corresponding to different parameter values. Interpolation between spectra is used to obtain spectral shapes between the grid points. Using our numerical code we perform a large number of simulations for different parameter values, spanning the grid Γ = 100, 150, 200, 250, 300, 350, 400, 450, 500 L 0,52 = 0.1, 0.5, 1, 5, 10, 50, 100, 200, 300 ε d = 0.01, 0.025, 0.05, 0.75, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4 which yields a total of 891 grid points. Note that other physical parameters often discussed in the context of subphotospheric dissipation, such as r d and r 0 , are derived quantities according to the expressions defined in section 2.1. Furthermore, the model obtains a normalisation and redshift as two additional parameters when it is implemented in XSPEC. These two parameters are kept constant during the fits, with the redshift given from observations and the normalisation set by the luminosity distance. As in A15 we let ε b = 10 −6 and ε e = 0.9, in order to test the scenario where almost all the energy goes into electrons and where we have no significant contribution from synchrotron. We note that this leaves ε b + ε e < 1. However, this simply means that the remaining energy is never dissipated. The main difference to the parameter space employed in A15 is that we now have set τ = 35 instead of letting it be a free parameter. This represents a scenario with the dissipation occurring moderately deep in the outflow. There are several reasons for this choice. Firstly, testing shows that most fits prefer high values of τ, with very few spectra being described by τ ∼ a few. However, at τ 15 this parameter turns out to be relatively insensitive. Changes in τ at these values usually result in small changes of the fit statistic, predicted spectral shape in the part of the GBM energy range where the statistics are best, and small changes in other model parameters. Additionally, τ cannot always be given a meaningful physical interpretation. Due to significant pair production, it does not directly correspond to the number of scatterings in the outflow, it also does not correspond directly to r d . Finally, keeping τ fixed decreases degeneracies in the model. Compared to the model in A15, we have also expanded ε d to lower values, from 0.1 down to 0.01, and made all parameters more tightly spaced in the grid. These improvements from the previous model reduce uncertainties pertaining to interpolation in the grid (see discussion in appendix B).
With these modifications we have captured most of the physically relevant parameter space for the three free parameters of our model. The main exception is Γ, which could in principle lie outside the range 100-500, although most Lorentz factors are expected to lie within our parameter range, see e.g. Ghirlanda et al. (2018). The upper bound for ε d is set by dynamical concerns, as discussed in section 2.1, and the lower bound by efficiency requirements. To a good approximation, the fireball luminosity, L 0,52 , relates to the observed isotropic equivalent luminosity as L iso ∼ L 0,52 ε d . Hence, the maximum fireball luminosity is L 0,52 < L iso,max /ε d , for a maximal observed luminosity, L iso,max . Atteia et al. (2017) find the limit E iso < 3 × 10 54 erg for observed isotropic energies. Assuming a typical duration for long GRBs of ∼ 50 s, (Bhat et al. 2016b), we obtain the limits L 0,52 < 600 and L 0,52 100, for the minimum ε d = 0.01 and a more typical ε d = 0.05, respectively. Thus, considering variations in ε d and burst duration, the limit, L 0,52 ≤ 300, we have set is reasonable, though not absolute. In section 5.2 we discuss the implications of our choices of parameter limits.
In A15 we referred to our table model as DREAM (Dissipation with Radiative Emission as A table Model). Since the goal is to make DREAM publicly available, we adopt a naming convention to help discriminate between and refer to the different versions of the table model as we develop it. The first number relates to the physical scenario considered when generating the model spectra, using the numerical code from Pe'er & Waxman (2005). The physical scenario decides the parametrisation in the numerical code, without changing the underlying numerical scheme or the relevant physics. The second number denotes which grid has been used, including minor revisions such as post-processed treatment of adiabatic cooling. With this naming convention, the model of A15 becomes DREAM1.1, and the model herein is DREAM1.2. Throughout the paper, "model" is, unless otherwise stated, referring to the table model, DREAM1.2. "Physical scenario" is used to refer to the underlying physical model (subphotospheric dissipation by internal shocks with negligible magnetic fields). Any references to a "code" refers to the numerical code presented in Pe'er & Waxman (2005).
In Fig. 1 we exemplify some spectral shapes as functions of the free parameters. As the dependency of spectral shape on our parameters is non-linear and we have a multidimensional parameter space, it is not possible to give an exact accounting for how the spectral shape changes with any single parameter. However, there are some clear features that are typical to certain parameters. Using Fig. 1 we give a qualitative overview of the main impact of our parameters on the spectral shapes.
Generally, the spectrum will consist of a BB peak and an inverse Compton peak at higher energies. For a small amount of dissipation we expect the BB peak to dominate and to yield a single-peaked spectrum. When the dissipation increases the inverse Compton peak becomes more pronounced and we can obtain a single-peaked spectrum with the peak at high energies or, if the peaks are comparable in strength, we obtain a flat, 'top-hat' shape, or a double peaked spectrum. For low values of ε b , as we consider here, the spectrum is dominated by inverse Compton scattering and the contributions from synchrotron radiation are negligible. Another feature is the pair production peak at ∼ Γm e c 2 /(1 + z).
The parameter with the clearest effect on the spectra is the fireball luminosity, L 0,52 , which mainly sets the flux, since it is positively correlated to the number of photons in the spectrum and spans several orders of magnitude. Additionally, L 0,52 has the effect of setting the distance between the peaks in the spectrum caused by inverse Compton scattering. This is caused by the electron temperature being independent of L 0,52 , while the temperature of the seed BB is inversely dependent on L 0,52 (Pe'er & Waxman 2005).
One of the major effects of Γ is the Doppler shift, which sets the observed energies of the spectrum. Additionally, the initial BB temperature increases as a function of Γ, which sets the position and normalisation of the first peak in the spectrum. Increasing Γ also yields a higher pair multiplicity, which is why we see the spectra with increasing Γ in Fig. 1 get increasingly thermalised as Γ increases. Note that this is different from the commonly assumed solution to the compactness problem, where an increasing Γ suppresses photonphoton annihilation. Due to the assumption that r 0 = r d /Γ 2 , Γ will alter quantities such as volume and electron number, resulting in a net effect of positive correlation between Γ and pair multiplicity. Note that Γ also contributes to setting the number of photons, but with an inverse dependence. However, due to its limited range and its other effects on the spectrum, it is the luminosity which works as the main driver of the predicted flux.
The last free parameter of our model, ε d , mainly changes the comoving electron temperature, θ, (Pe'er et al. 2006). This sets the normalisation of the inverse Compton peak, with a higher ε d resulting in a stronger peak. However, the position of the this peak is not sensitive to ε d for much of the parameter space. This is because the first upscattered photons may pair produce and yield a population of cooler electrons, usually outnumbering the baryonic electrons by a factor of several. The observed inverse Compton peak in the spectrum is a result of multiple scatterings with electrons of a wide energy distribution, making its position almost independent of the initial electron temperature.

Additional scatterings and adiabatic cooling
In A15 the numerical code was run only as long as there was dissipation, i.e. between r d and 2r d . This takes care of the most important physics and relevant processes, but it does not include the last of the scatterings, between 2r d and r ph . Running the code all the way to the photosphere from an optical depth of τ = 35 is not feasible with the current set-up. However, it is sufficient to run the code a fraction of the way to the photosphere to capture most of the spectral changes occurring after the dissipation. In Fig. 2 we show examples of spectral shapes as a function of where the simulation ends. It is clear from the figure that there is little change already after 2.3r d In this work we have chosen to let the code run between r d and 2.3r d . We note that the factor of 2.3 was found empirically to be a good trade-off between computational efficiency and capturing changes in the spectral shape.
As stated in A15, the code does not in its current form include adiabatic cooling. However, we would expect adiabatic cooling to yield a shift of the entire spectrum to lower energies., see e.g. Pe'er & Waxman (2004). This effect is caused by a continuous cooling of the photons as they do work on the outflow for as long as they scatter, i.e. roughly as long as r < r ph in our scenario. In a more realistic scenario we would also expect the photons to start decoupling before r ph , mainly leading to a broadening of the spectral shape (Beloborodov 2011;Bégué et al. 2013). The photons continuously isotropise during collisions by converting internal photon energy to bulk kinetic energy of the outflow, and thus stay isotropic in the rest frame. As long as we have a thermal spectrum this results in a constant shift of the spectrum to lower energies. For non-thermal spectra, each scattering will also distort the spectral shape. We begin by considering only pure adiabatic cooling. As it is a constant shift of the spectrum, this effect may be treated by analytic post-processing. We introduce an adiabatic cooling factor a = 2(τ d ) −2/3 , where τ d is the optical depth at the end of the dissipation, with which we simply shift the spectrum. Hence the number of photons at a given energy, N(E) is shifted such that N(E) = N ad (E ad ), with E = aE ad and where E ad is the shifted energy. For more details on the adiabatic cooling factor, including its limitations, see Beloborodov (2011) and Bégué et al. (2013). The constant shift using a is an approximation to taking into account all scatterings up to the photosphere.
In Fig. 3 we show an example of an original model spectrum, as calculated from our numerical code, compared with the spectrum after having applied the adiabatic cooling factor, as well as the spectrum after a full treatment of all scatterings, adiabatic cooling and geometric effects up to the photosphere. The full treatment was performed using a numerical code from Lundman et al. (2017), which dynamically couples hydrodynamics to Monte Carlo radiative transfer inside a spherical outflow. The radiation can then be followed as it cools adiabatically, and self-consistently decouples from the flow at a fuzzy photosphere. This code can be used to study the spectral evolution after dissipation. However, the code is computationally very expensive, which prevents us from using it for all spectra. Fig. 3 illustrates that our approximation of a constant shift yields good agreement of the integrated model flux. However, we note a small discrepancy between the spectrum with the constant shift approximation and the full treatment in the low energy slope < 20 keV. This effect is a result of the full treatment taking a fuzzy photosphere and geometric effects into account. We thus conclude that the approximate treatment of the adiabatic cooling is good, and that it represents a significant improvement over the original model spectrum, but that the inclusion of geometric effects would lead to a slightly softer low-energy slope.

Data sample
Our sample contains the brightest GRBs with known redshifts observed by Fermi before 2016-06-01. The known redshift 1 helps us fix the normalisation parameter of the model spectrum via the corresponding luminosity distance, instead of leaving it as a free parameter. This is important because we want to be able to test the model's ability to correctly predict the GRB flux. We chose a fluence cut of > 10 −5 erg cm −2 , in order to allow us to perform a time-resolved analysis with good signal strength. The resulting sample contains 36 GRBs, listed in Table 1.
For all GRBs we use Time Tagged Event (TTE) data from the Gamma-ray Burst Monitor (GBM) NaI and BGO detectors (Meegan et al. 2009). We select detectors by the scheme outlined in Bhat et al. (2016a), choosing up to three NaI detectors with an angle towards the source lower than The thick lines represent the output spectrum from the code, from which we produce the table model. For illustrative purposes we also show the seed photon BB distribution, plotted in the same line style as the corresponding spectrum, but with thinner lines. The red, vertical lines show the Fermi GBM energy range, which is the fitted energy interval for the majority of GRBs. It should be noted that the plotted spectra include the adiabatic cooling factor a, introduced in section 2.3.
10 1 10 2 10 3 10 4 10 5 10 6 Energy (keV) 10 0 10 1 10 2 10 3 10 4 Figure 2. Examples of the spectral shape for different stopping radii of the numerical code. Note that the dissipation stops at 2r d for all spectra. The shape asymptotically approaches the spectrum at r ph , with most of the spectral evolution occurring early. These spectra were produced using the model parameters L 0,52 = 300, Γ = 300 and ε d = 0.4. 60 • as well as the BGO detector with the lowest angle towards the source. The energy range fitted was 8-1000 keV (NaI) and 200 keV -40 MeV (BGO).
There are also LAT-LLE data available at high energies for nine bursts in the sample (GRB 080916C, GRB 090323A, GRB 090902B, GRB 090926A, GRB 110731A, GRB 130518A, GRB 141028A, GRB 150403A and GRB 160509A). The LLE data are a type of LAT data designed for the use with bright transients to bridge the energy range of the LAT and GBM (Pelassa et al. 2010). The LLE data (pass7) were obtained using gtburst version v02-01-03p50 and Fermi Science Tools version v9r32p5. The fitted energy range is 30-100 MeV. As noted in section 2.3, it is in this energy range our model has the largest uncertainties. However, we note that the inclusion or exclusion of the LLE data has little to no effect on the results (see appendix C for a discussion on the impact of the LLE data on our analysis).  The spiky behaviour at high energies in the curve of the full treatment is numerical noise and the small peak in the original and artificially cooled spectra is the pair annihilation peak.

Time-resolved spectral analysis
Our spectral analysis is time-resolved, using Bayesian blocks (Scargle et al. 2013) as binning method 2 . The algorithm is performed on unbinned TTE data, with the false alarm rate parameter p 0 = 0.05. We perform the binning using the NaI detector with the lowest angle of incidence. All signal-tonoise ratio (SNR) presented are also calculated from this detector. The background is fitted with polynomials. We account for repointing of the spacecraft and temporal evolution in the background by using rsp2 responses when available.
Bayesian blocks with an SNR cut is the most appropriate choice of binning given our model. The model does not include any time-evolution of the physical properties of the jet, meaning that any spectral evolution observed is assumed to be caused by central engine activity. Hence we assume constant central engine activity and no spectral evolution in each time bin we analyse. A constant central engine activity may be expected to give rise to a constant Poisson rate of photons in the light curve, which is what we recover with the Bayesian blocks scheme. However, this method may yield bins with a very low SNR. Thus, an SNR cut of 4 was placed on the bins in accordance with what we found in the validation tests (see appendix A). We use an expression for the SNR as derived in Vianello (2018), for a Poisson measurement and a Gaussian background. Previous work by e.g. Yu et al. (2016) also employs an SNR requirement when performing time-resolved spectral analysis. The large difference in the level of the SNR cut between our work and theirs is a result of the different expressions used for calculating the SNR. See also Burgess (2014) for a discussion about different binning techniques, and Worpel & Schwope (2015) for general transient timing using Bayesian blocks.
In Fig. 4 we show an example of a light curve binned with Bayesian blocks, including 4 bins excluded at SNR< 4. Our burst sample, defined in section 3.1, contains two bursts which yield no bins above the SNR threshold and are thus not analysed. The bursts which did not make the SNR cut are GRB 090516 and GRB 140423A. In total 145/779 bins are removed before analysis due to low SNR, giving us 634 time bins to analyse.
The analysis is performed using PyXspec, a python implementation of HEASARC's XSPEC 12.8.1g, with pgstat statistics and with PHA (Pulse Height Analyser) data. We fit both our table model, DREAM1.2, as well the Band function, since the latter is used as a standard function in the field. We want to examine possible correlations between Band function parameters (i.e. the low energy slope, α, the peak energy, E p , and the high energy slope, β), and the model parameters of DREAM1.2.
We note that the finite resolution of the  Figure 4. Example of Bayesian blocks binning with SNR cuts applied. The burst in this example is GRB081121. The red vertical lines outline the time bins found using Bayesian blocks. Time bins 8, 10, 12 and 13 (corresponding to times 16.1 − 17.4, 17.9 − 18.8, 19.5 − 21.8 and 21.8 − 47.7 s) are excluded from the analysis due to the SNR cut at 4. These bins are shaded red. The light curve is produced with 0.3 s bins, using the NaI 11 detector, which is the same detector used to bin the data. The figure presents raw detector counts from all energy channels, with no background subtraction or off-axis effects accounted for.
the regular fitting with Maximum Likelihood. Examining the points in the parameter space where we expect the largest resolution-related uncertainties, we find the median of the uncertainties to be less than 4 − 8 per cent, depending on the parameter. L 0,52 and ε d have the largest uncertainties. This evaluation is presented in appendix B. These uncertainties are not propagated into further analyses.
We have tested for effective area corrections between the detectors and find that the results do not change significantly if we allow for this correction or not. The size of the errors on fitted parameters are also not systematically or significantly affected. We present the fit results obtained without effective area corrections.
We assume standard Λ-CDM cosmology with a Hubble constant of H 0 = 67.3, a cosmological constant, Ω Λ = 0.685 and a matter density Ω M = 0.315, (Planck Collaboration et al. 2014). All uncertainties are 1 σ unless otherwise stated. Statistical uncertainties on parameter estimates are obtained using the XSPEC error command, where model parameters are varied until the fit statistic changes by the desired level.

Goodness of fit
We wish to have a metric by which to judge if a fit should be rejected since we do not expect a fit of poor quality to give physically meaningful information. For this purpose we employ a Monte-Carlo (MC) parametric bootstrap method, using XSPEC. For each fit to real data we simulate 5000 spectra from the given best-fitting parameter values, using XSPEC's built-in function fakeit. The spectra are simulated using the response and background files of the original data. We then re-fit the model to the fake spectra and sample the fit statistic.
We perform a goodness of fit (GOF) test by comparing the pgstat value of the model fitted to the real data with the pgstat distribution obtained from the MC simulations. We choose to reject fits which have a fit statistic outside the 99.7 per cent central confidence interval of the sampled distribution (corresponding to 3σ). We performed tests with cuts corresponding to values between 1 and 5 σ, and found that it does not affect the overall shape of the distribution of parameter values in the accepted fits. This indicates that our results are qualitatively robust to the level of this cut.

RESULTS
After binning the data in Table 1 and performing the SNR> 4 cut we fit DREAM1.2 to the data. We discard all fits with a parameter on the boundary and perform our GOF test on the remaining fits. Out of 634 time bins with SNR > 4, 267 have no parameter on the boundary of the parameter space. 171 of these fits pass the GOF test and constitute our sample of accepted fits, which we analyse further. This corresponds to approximately 27 per cent of analysed spectra being accepted, with 10 bursts having at least 50 per cent accepted fits. The rejected fits and their implications on our model assumptions are discussed in section 5.2. See appendix E for a table summarising the number of accepted and rejected time bins for each burst.

Best-fitting parameters
In Fig. 5 we present the global parameter distributions for the free model parameters of our accepted fits, as well as the dissipation radius, r d . We have included Gaussian kernel density estimates in the histograms to aid with visual interpretation. The histograms span the entire parameter space, except for r d , which has allowed values in the range 3 × 10 8 − 4 × 10 15 cm. We include r d since it may give information about the dissipation scenario and/or the GRB's relation with the progenitor. Γ is almost evenly spread in the range 100 − 250, with a small peak at Γ ∼ 120. L 0,52 has the majority of its fits in the range 10 − 100.   The introduction of kernel density estimates scales the histograms so that they represent probability distributions, and the y-axes are normalised accordingly.

Correlations
It is of interest to search for correlations between the different best-fitting parameters of the model, as well as between our model parameters and the observed flux, isotropic equivalent rest frame luminosity, L iso,z , and redshift. We also examine correlations between our model parameters and the corresponding Band function parameters. Fluxes and luminosities were calculated with k-correction (Hogg et al. 2002), using a cutoff power law in the stellar rest frame 3 energy range 1 keV -10 MeV. When performing the correlation analysis we consider only fits with well constrained parameter errors, i.e. where the statistical errors on all parameters obtained during fitting are contained inside the parameter space. This leaves us with a subset of 114 bins for the correlation analyses involving E p,z and 118 for the other correlations.
To search for correlations, we use a hierarchical Bayesian model for fitting a straight line to data with errors in both the dependent and independent variables, described in Kelly (2007). We use this technique to find estimates of the linear relation parameters, as well as estimates for the population linear correlation coefficient, ρ. The errors on derived parameters have been obtained with standard error propagation.
As expected, we find a very strong correlation between the logarithm of the model luminosity, L 0,52 , and the logarithm of the observed luminosity, L iso,z . If we instead consider log(ε d L 0,52 ) − log L iso,z , the correlation becomes even tighter. This correlation serves as a sanity check. Additionally, we find a correlation between the model parameters Γ and log L 0,52 , as well as between log ε d − log E p,z , where E p,z is the Band function peak energy in the stellar rest frame. All correlations are presented in Fig. 6. The estimated linear correlation coefficients, including their 1 σ credible intervals, as well as the regression slopes are also shown in the figure.
Finally, we also find a correlation between log E p,z and Γ, which is not present across the entire sample, but instead only exists within some bursts. We have only investigated this correlation in bursts with at least 5 accepted time bins. In Fig. 7 we present the two bursts in which the correlation is significant. The correlation is considered significant if the 95 per cent highest posterior density region of the marginal posterior for the correlation coefficient does not include 0. The correlations are discussed further in section 5.4. Notably, there are no correlations with the Band function α or β. The implications of this are discussed in section 5.3.

Accepted versus rejected fits
To determine why only 171 out of 634 fits are accepted, we consider the parameter distributions of best-fitting parameters from our model, the physical quantities flux, redshift and L iso,z , as well as parameter values of the corresponding Band function fits.  Examining the parameter distributions of rejected fits, we find that 302 out of 463 rejected fits have L 0,52 = 300. The remaining 161 spectra, which have not reached the upper bound in L 0,52 , exhibit a greater variation in their parameter values. Additionally, out of these 302 spectra, 260 also have ε d = 0.4. Furthermore, manual inspection of spectra shows that the fits that hit the upper boundary in L 0,52 and ε d are significantly worse than those that do not.
When analysing correlations between physical quantities and whether a fit is successful or not, we find that bursts with higher redshifts have a lower fraction of successful fits, and similarly for the flux. Even tighter is the correlation between a high L iso,z and a high fraction of rejected fits. The latter relation is shown in Fig. 8. These results, together with the pegged values of L 0,52 in most rejected fits, demonstrate that it is the model failing to reach the required normalisation, i.e. the predicted flux, which is the cause for the majority of the rejected fits. There is no discernible pattern in the best-fitting Band α or β parameters for rejected fits. This is an important result which we discuss further in section 5.3. However, the spectra with the 26 highest value of E p are not successfully fitted with DREAM1.2, and there is a decrease in the fraction of accepted spectra from E p 300 keV. This suggests that DREAM1.2 works better to describe spectra with lower values of E p . In Fig. 9 we show examples of typical accepted and rejected fits applying DREAM1.2 to our sample. This clearly illustrates the normalisation problem for the rejected fits.

Time-evolution within bursts
For a physical model the time evolution may serve as a consistency check (since we expect physical parameters to behave in some systematic way) as well as a way to assess the model. Many of the analysed bursts do not have a sufficient number of bins with acceptable fits for a study of the time evolution. 11 bursts in our sample have more than 4 surviving bins. Only these are considered for the analysis of the temporal evolution of best-fitting parameters. See table E1 for a summary of the number of accepted time bins for each burst.
For these bursts, there are some tentative systematic parameter evolutions. As expected, we find that L 0,52 approximately follows the light curves. As mentioned in A15, we observe a decreasing value of Γ for some bursts, whereas there is no apparent pattern in others. The dissipation radius, r d , generally increases during a burst. The level of dissipation, ε d , appears to have no consistent temporal behaviour.
In Fig. 10, we provide an example of the time evolution of L 0,52 , Γ, ε d and r d for the accepted fits of GRB 150821A, plotted together with the flux light curve. L 0,52 correlates with the light curve, as seen in most bursts. In this particular burst Γ also tends to follow the light curve, and ε d appears to anti-correlate with it. Finally, r d remains almost constant throughout the burst duration. We provide equivalent plots for all bursts with more than 4 accepted time bins in appendix D. We discuss the implications of the temporal evolutions in section 5.5.
In the flux light curve of Fig. 10, we note that there are a few very thin time intervals with unusually large or small fluxes. This is a consequence of our Bayesian blocks binning and these fluxes are not reliable. However, due to their very short durations, these bins have not passed the SNR cut and are thus not part of the analysis.

DISCUSSION
We begin by relating to the discussion of A15 and point out some of the key differences between the two works. Secondly we consider the rejected fits, to see where the model fails and discuss how this relates to the underlying assumptions. We also discuss the relation to the Band function, correlations as well as other physical models and previous works.

Comparison with A15
As shown in section 2.1, there are several differences between the model presented in this work and that in A15. Although the models are produced using the same numerical code, the improvements we have made since A15 have changed the model predictions. The single largest improvement is the introduction of adiabatic cooling (see section 2.3). One of the conclusions in A15 was that τ is an important parameter which sets the distance between the two peaks in the spectrum and decides if the spectrum appears single or double peaked. In this work we have instead chosen to keep τ fixed and argue that it has low impact for most fits. There are two reasons for these differences. Most importantly, the introduction of the adiabatic cooling has increased the degeneracy between τ and L 0,52 . As a result, the model generally prefers higher values of τ and becomes increasingly insensitive to high values of τ. Secondly, low values of τ predominately lead to double peaked spectra, but we find that our sample mainly consist of single peaked spectra. This is not surprising given the small number of GRBs with strong double peaks reported in the literature (Guiriec et al. 2011;Axelsson et al. 2012). We found that most spectra in our sample prefer τ 15, allowing us to keep it fixed without changing the number of accepted fits significantly. Individual GRBs may still prefer different values of τ and it may be warranted to keep τ as a free parameter in some studies.

Implications and analysis of rejected fits
The fact that 463 time bins are not well described by our model shows that the simplest version of the internal shock scenario below the photosphere with negligible magnetisation employed here does not fully account for GRB prompt emission. As shown in section 4.3, the main reason for fits For GRB 091127 data from NaI 6,7 and 9, as well as BGO 1, are shown in black, red, blue and green, respectively. For GRB 130518A data from NaI 3, 6 and 7, as well as BGO 0 and the LLE data, are shown in black, red, blue, green and magenta, respectively. The data have been re-binned to errors of 3 σ for visual clarity. The residuals are produced using the XSPEC, 'delchi' command (for the statistic used here this does not correspond to contributions to the statistic). The best-fitting parameters are ε d = 0.07 +0.02 −0.01 , L 0,52 = 41 +4 −4 , Γ = 188 +9 −11 , for GRB 091127, and ε d = 0.40, L 0,52 = 300, Γ = 150, for GRB 130518A. Note that the parameter values found for GRB 130518A represent a model which does not fit the data.
being rejected is that the model cannot reproduce the observed flux. Below we discuss how this result is affected by our model assumptions and the limitations on the parameter space.
We start by considering these results in relation to the luminosity distribution of GRBs. In Fig. 11 we compare the maximal L iso,z of an accepted fit (7.5 · 10 52 erg s −1 ), to the L peak distribution from Yonetoku et al. (2010), where L peak is the value of L iso,z for the 1 second peak spectrum of each burst. 27/101 bursts in this sample have an observed L peak > L iso,z,max . Thus, to the degree that the sample in Yonetoku et al. (2010) is representative of the GRB population, the current version of our model is incapable of describing the brightest ∼ 30 per cent of GRBs, only considering its low predicted observed luminosity. We note that this is only a first-order estimate. It does not imply that the remaining 70 per cent of all GRBs can be described by the model, only that this fraction is not immediately ruled out by the luminosity.
One possibility to accommodate the problem of insufficient predicted flux would be to expand the parameter space, e.g. to higher values of L 0,52 . As noted in section 2.2, the boundary of L 0,52 = 300 that we have imposed is not absolute, but nearing the limit of what is physically reasonable. To investigate this we constructed a model including the parameter value L 0,52 = 1000. This model does yield additional accepted fits (219 in total), but is not sufficient to fully solve the issue, with many spectra still having higher fluxes than the model can account for. The highest L iso,z we can produce with DREAM1.2 is max(L 0,52 ) · max(ε d ) · ε e · a ≈ 4 · 10 53 , where a is the adiabatic cooling factor presented in section 2.3. As shown in Fig. 8, the fraction of accepted fits goes down significantly well below this value. This is because the spectral shape changes and the model becomes less flexible when L 0,52 and ε d are forced to high values, and the fits are instead discarded in the GOF test. Larger values of ε d are possible in principle, but we are not able to test this with our current set-up without ignoring the possibility of a re-acceleration phase. Additionally, as shown in Fig. 5, most fits prefer low values of ε d , with only 18 per cent of fits having ε d > 0.2, suggesting that higher values of ε d may not be required.
We note that the distribution of Γ presented in Fig. 5 peaks at the lower end of the parameter range. However, there are very few failed fits that actually hit the the lower boundary, making it unlikely that lower values of Γ will resolve the flux issue. The current range of this parameter is also in good agreement with other observations (see section 2.2). Changing the parameters that were kept fixed in this study may also help resolve the issue. Most interesting are ε b and τ. In particular, a higher magnetisation will increase the number of photons provided by synchrotron emission. We have performed fits with ε b = 0.1 and found that it does not result in a significantly larger number of successful fits. A study including other values of ε b needs to be conducted before any conclusions regarding the magnetisations can be drawn. In the case of τ, we note that lower values lead to more photons and less adiabatic cooling, which results in higher fluxes. However, we have investigated this by using a model where τ is allowed to vary in the range 1 − 500 and found that it does not significantly change the number of accepted fits. It is possible that dissipation in the optically thin region would provide additional accepted fits, but this is not a part of the subphotospheric dissipation scenario that we are testing.
Another option is to consider an alternative dissipation scenario. This would allow us to increase the predicted flux while still considering subphotospheric dissipation in an outflow with low magnetisation. We therefore examine our model assumptions, particularly the internal shock scenario, which implicitly sets the nozzle radius as r 0 = r d /Γ 2 . This in turn sets the temperature of the initial BB, which, given a specific luminosity, sets the number of photons in the outflow. One way to accommodate the flux issue may therefore be to remove the assumption r 0 = r d /Γ 2 , instead making r 0 a free model parameter. This will be investigate in a future work. Alternatively, by invoking additional dissipation deeper in the outflow, in the Wien or Planck zone, it would be possible to increase the number of seed black body photons (Beloborodov 2013).
Another possible origin for rejected fits is that there are additional emission components present. In our model we have assumed a one-zone emission and no synchrotron emission. There are cases of GRBs where, e.g., an additional power-law component is used to fit the data, such as Arb. norm Figure 11. L peak distribution from Yonetoku et al. (2010), showing the k-corrected isotropic equivalent luminosity in the 1 second peak interval spectrum, L peak , from 101 bursts. The red line indicates the maximum L iso,z for which we have an accepted fit with DREAM1.2 (7.5 · 10 52 erg s −1 ). The blue line shows a Gaussian kernel density estimate of the distribution.
in GRB 090902B, (Ryde et al. 2010). This particular burst is in our sample, but only contains rejected fits, due to the model under-predicting the flux. This is discussed further in section 5.5.

Relation to the Band function
As noted in section 4.2, there is no relation between the Band function α and β parameters and whether the DREAM model provides an acceptable fit to the spectra. A plausible explanation for this is that there is spectral complexity not captured by the Band function. In Fig. 12, we show the distribution of best-fitting Band function parameters for the distribution of accepted spectra, as well as for all analysed spectra. We note that the accepted fits covers essentially the same ranges as the full distribution.
By possibly ignoring features in the spectrum, incomplete, or even misleading, spectral information might be retrieved when fitting with the Band function, leading to erroneous conclusions. For example, as pointed out in A15, the possible existence of two breaks in the spectrum could imply that we need a more nuanced notion of E p , meaning that a single peak energy might not be a relevant characteristic for all spectra. To illustrate these issues we show an example in Fig. 13 of fits to the spectrum of GRB 150821A in the time interval 21.1 − 47.3 s using DREAM1.2 and the Band function, plotted together. Particularly noteworthy is that α = −1.2 is not what is usually associated with the lowenergy slope of photospheric models. The key point is that it does not need to be, because, as shown in Fig. 13, values of α may be averages of more complex spectral slopes.
Comparing values of the Band function α with the asymptotic low-energy slope of the thermal radiation is only meaningful if a) the asymptotic behaviour of the BB is contained in the fitted energy range and b) if the data actually represents a power law. Neither of these conditions are necessarily true. We show here, as well as in A15, that a) the fitted energy range is above the asymptotic slope of the seed BB and b) that the data may be equally well described by a non-power law behaviour at low energies.
We stress that Fig. 13 does not show that the presented spectrum must be curved at low energies. However, the forward folding and measurement uncertainties make it hard to assess any model by comparing it to fits with another model. We emphasise that physical models must instead be assessed based on direct comparison to data, as well as their physical predictions. These concerns, as well as the shortcomings of the Band function for describing the underlying data, are also discussed in detail in, e.g., Burgess et al. (2014aBurgess et al. ( , 2015 and Burgess (2017).
Given that photospheric emission is usually associated with very hard values of α, it is an important result that the model can describe spectra with α as soft as -2 (see Fig. 12).
Our finding that the model can describe the full range of Band parameters is with the approximation of no geometric broadening. The effect of geometric broadening is in general a softening of the low-energy spectral slope (cf. Fig. 3).

Correlations
The correlations between the model parameters of DREAM1.2 and the Band function have no physical interpretation. Instead, the ε d − E p,z and Γ − E p,z correlations can be understood in terms of the position of spectral peaks. In DREAM1.2, most spectra are slightly double peaked, which means that E p,z corresponds to some average of the two peaks. For high ε d , E p,z will be found closer to the highenergy peak (cf. section 2.2 and 5.3) and thus we see a positive correlation between these parameters. As can be seen in Fig. 1, Γ yields more complex changes to the spectral shape. This leads to, not surprisingly, that the Γ − E p,z correlation is less prevalent in the sample than the ε d − E p,z correlation.
Regarding the log L 0,52 − Γ correlation, we note that an L iso − Γ correlation has previously been reported by Lü et al. (2012) and Ghirlanda et al. (2018). Lü et al. (2012) argue that the positive correlation is consistent with the idea of a hyper-accreting black hole central engine, where the jet is powered by a neutrino-driven dominated accretion flow. In these studies L iso is the time-integrated luminosity. Lü et al. (2012) calculate the Lorentz factor used to find the correlation with the afterglow onset method (Sari & Piran 1999), whereas Ghirlanda et al. (2018) compare several different methods. The correlation we present here is self-consistent and we find it independently of previous works, which used completely different methods to calculate Γ. Additionally, our model also allows for an alternative origin of the correlation, as being a result of the radiative transfer process in the outflow. This interpretation is independent of assumptions about the central engine, other than the assumptions from the fireball model.
In appendix A we demonstrate that the log L 0,52 − Γ relation is not caused by degeneracies in the model, or by intrinsic correlations. Although Γ and L 0,52 are related in the model through their mutual effect on r d , there are no assumptions on the value of r d , and the two parameters are not a priori functionally correlated or dependent. However, as shown in appendix A, when excluding spectra with low SNR we remove spectra with simultaneously high Γ and low L 0,52 . This indicates that the correlation may be partially due to selection effects. The impact of selection effects on GRB parameter correlations have been studied before. Heussaff et al. (2013) consider the E peak − E iso correlation (E peak and E iso being the peak energy of the prompt spectrum and the isotropic equivalent energy, respectively), and find that there is a bias against detecting GRBs with large E peak and low E iso due to their lower SNR. Ghirlanda et al. (2018) also consider the impact of selection bias and find that the L 0,52 − Γ correlation they observe is not a result of selection effects. Similarly, it seems like the selection effect we find in appendix A cannot fully account for the observed log L 0,52 − Γ correlation. However, a more detailed statistical analysis where selection effects are specifically modelled for will be carried out in future work.

Implications for the physical properties of jets and progenitors
Even though most fits are rejected (463 out of 634 fits), we have shown that the problem occurs predominately for bright bursts and that our sample is biased towards bright bursts. The majority of all bursts have luminosities which are consistent with the model (see section 5.2). Additionally, there are 6 bursts in our sample where more than 60 per cent of the analysed time bins are accepted (see table E1). This justifies considering the implications of the best-fitting parameters of the accepted fits. From Fig. 5 we see that r d is generally distributed around 10 12 cm. All values of r d from successful fits are above the ∼ 10 11 cm radius typically associated with Wolf-Rayetlike (WR) stars, usually assumed to be the progenitors of long GRBs. As pointed out in section 4.1, values of r d found from the fits inhibit only a fraction of the allowed parameter space, with values as low as 3 × 10 8 cm being allowed. Hence, the relatively narrow interval in which we find r d indicates that the dissipation is not caused by interaction with the progenitor star itself. These radii could instead be associated with internal or external shocks, as well as interactions of the outflow with a baryonic shell lifted from the WR star, as found by  and .
Furthermore, we find that most values of ε d lie around a few percent, which is consistent with what is typically assumed for internal shocks (Mochkovitch et al. 1995;Kobayashi et al. 1997). Although we use an internal shock scenario for our dissipation process, there is nothing in our numerical code which forces the efficiency to this level, which makes the result intriguing. Regarding the nozzle radius, r 0 , we find, from the relation r 0 = r d /Γ 2 , that r 0 typically takes on values r 0 ∼ 10 7 − 10 8 cm. This result lies in proximity of the usually assumed value of r 0 ∼ 10 7 cm, which is obtained when a variability time-scale of milliseconds is assumed for the central engine. Similar results (10 6.5 r 0 10 9 ) have been found also by Pe'er et al. (2015). Note that, as for ε d , there is nothing in the fitting procedure or code which forces these values. As mentioned in A15, the time evolution of Γ can be interpreted to correlate with central engine activity and could in principle be used to constrain accretion scenarios. Having Γ decrease during a burst implies a central engine luminosity that decreases faster than the mass loading, or even an increasing mass loading.
When studying the temporal evolution of our model, we note that most GRBs have accepted and rejected fits interspersed, cf. Fig. 10 and corresponding figures in appendix D. We do no not interpret this as subphotospheric dissipation being turned on and off throughout the GRB. This instead points to the fact that our model is a limited implementation of the physical scenario considered, having only three free parameters. Although additional free parameters or expanded parameter space are not expected to be able to fully solve the issue of the under-predicted flux (see section 5.2), it can likely improve fits for bursts with moderate luminosity.

Other physical models and previous work
In this section we consider previous work on subphotospheric dissipation, as well as other physical models, particularly those which have previously been fitted to bursts in our sample. GRB 090902B is one of the best examples of a burst with photospheric emission and it has been the focus of several papers. Ryde et al. (2010) describe GRB 090902B by a two-component scenario by modelling the photospheric emission as a multicolour BB and adding an additional power-law component for the non-thermal emission. Both the thermal and non-thermal emission are modelled empirically. Pe'er et al. (2012) also describe GRB 090902B by a two-component scenario. They use an empirical multicolour BB to model the photospheric emission component. Using the numerical code of Pe'er & Waxman (2005), this multicolour BB is then injected as the seed photon population for optically thin dissipation to produce the non-thermal power law component.
In contrast, in this work we tried fitting GRB 090902B with a single component model from the photosphere. We were unable to produce any successful fits using this model, due to insufficient predicted flux. We also tested fitting this burst with a two component scenario, where the thermal component was produced by DREAM1.2, and the nonthermal component was represented by a power law. Also in this case the thermal component obtained an insufficient flux. The conclusion from our fits to GRB090902B is therefore that the dissipation scenario we assume in DREAM1.2 is not the correct one for this burst. However, we point out that a photospheric origin of the main component in GRB090902B is still probable, but a different dissipation scenario is needed. For instance, GRB090902B was modelled by Vurm & Beloborodov (2016), assuming a completely different scenario for subphotospheric dissipation. One of the key differences between the model used by Vurm & Beloborodov (2016) and our model is that they consider a jet which starts to dissipate at r d ∼ 3 × 10 10 cm, which is below their saturation radius. This is in contrast to our dissipation radius which, by construction, is always above the saturation radius. Allowing photon production from thermal and non-thermal mechanisms deeper in the outflow yields a higher photon production, particularly around the Wien peak, which makes it possible to obtain a higher predicted flux without resorting to higher fireball luminosities.
Another burst in our sample which has been fitted with physical models is GRB 080916C. Zhang & Pe'er (2009) suggested a multi-component description of this burst by showing that there is an additional energy source above the photosphere. This burst was also described using the ICMART model of Zhang & Yan (2011). In our analysis, GRB 080916C could not be fitted due to insufficient predicted flux, and a comparison of model predictions cannot be made.
There are other physical models which have been fitted to or compared to data, such as the external shock model of Burgess et al. (2016) for GRB 141028A, and the synchrotron model of Zhang et al. (2016) for GRB 130606B. Out of these, only GRB 141028A is in our sample. However, only 2 of the 15 time-resolved spectra are accepted, while the others are rejected, mainly due to the high flux. It is clear that the scenario of localised subphotospheric dissipation cannot describe this GRB as a whole, and a comparison with Burgess et al. (2016) is not possible. Studies of larger overlapping samples of GRBs using physical models are needed in order to make more thorough comparisons with our work.
The fact that all spectral slopes present in the Band function parameter space can be fitted with our model is a clear indication of our model's strength and its relevance. The fact that the current version of the model clearly does not suffice to describe all observed GRBs is an important conclusion of this work. This, together with the indications that there are possible alterations of the dissipation scenario which may resolve some of the identified key issues, warrant further studies.

SUMMARY AND CONCLUSION
We have considered a model for subphotospheric dissipation as the origin of GRB prompt emission and fitted it to Fermi GRB data. The dissipation is localised and we assume internal shocks as the dissipation mechanism. Additionally, the model does not take into account geometric effects, a fuzzy photosphere or jet hydrodynamics. We consider the scenario where there are no significant magnetic fields presents. A table model, DREAM1.2, was created from simulations using the numerical code of Pe'er & Waxman (2005). Using DREAM1.2 we performed a timeresolved analysis of 36 bursts, using Bayesian blocks as a binning method. We analysed a total of 634 time-resolved spectra. Out of these we find that 171 spectra are well described by our model, passing the GOF test and having no parameters on the boundaries of the parameter space. This corresponds to an acceptance-rate of about 27 per cent, with 10 bursts having at least 50 per cent accepted fits.
From the successful fits we conclude that: • The model can fit all spectral slopes found in our sample in terms of Band function parameters. The spectra are dominated by Comptonisation and have negligible contributions from synchrotron radiation.
• Most values of ε d are in good agreement with the level of efficiency expected from the internal shock scenario at ε d ∼ 1 − 10 per cent. It should be noted that there is nothing in the underlying numerical code which ties ε d to these values, despite assuming internal shocks.
• We find a correlation between the model parameters L 0,52 and Γ. This correlation is a possible equivalence of the L iso − Γ correlation found by Lü et al. (2012) and Ghirlanda et al. (2018). Simulations show that the correlation is not a model artefact. However, we find that it may be partially due to selection effects.
• The fact that there is very little correlation between our model parameters and the Band parameters emphasises the importance of using physical models to analyse GRB data.
From the rejected fits we conclude that: • With 73 per cent of spectra not being successfully described using DREAM1.2 it is clear that the scenario of photospheric emission with localised subphotospheric dissipation by internal shocks and without magnetic fields, does not suffice to describe the majority of bursts in the sample.
• The main reason for fits being rejected is that the model under-predicts the observed flux. The inclusion of additional parameter space of the model, or the inclusion of significant magnetic fields may yield additional accepted fits. However, we do not expect this to be able to completely resolve the flux issue. Alternative scenarios of subphotospheric dissipation or additional emission components could instead be invoked to describe the brightest GRBs.
• The analysed sample is biased toward bright GRBs. Comparing the maximum L iso,z of time bins successfully described by DREAM1.2 to the L peak,iso distribution of Yonetoku et al. (2010) yields the estimate that ∼ 1/3 of the GRB population are too luminous to be described using DREAM1.2.
We conclude that our model cannot describe the full population of GRBs, but that it appears to be capable of describing a significant part of the population. Thus, if all GRB prompt emission is of photospheric origin, the scenario we have considered in this work cannot be the only dissipation scenario. In order to further constrain the dissipation scenario, a possible next step is to include Swift data at lower energies, outside the Fermi GBM energy range. Additionally, the parameter space of the model can be explored further, including the effect of significant synchrotron emission or optically thin dissipation.  Figure A1. Relative error as a function of the SNR. The SNR is spaced logarithmically in 50 steps in the range 0.5 − 40. It is the average of the relative error in each SNR interval which is plotted. The horizontal black line denotes 0 relative error. The colours and marker symbols indicate the different parameters as given in the legend. 0.5−40. As expected, the relative error is 0 at large SNR, but increases when the SNR is low. The fact that the relative error becomes very large for L 0,52 and ε d at SNR 2, whereas it remain relatively low for Γ, is due to the fact that Γ can only vary over half an order of magnitude, whereas L 0,52 and ε d vary over several orders of magnitude. We note that there is a significant decrease in the relative errors between SNR 2 − 4 and that the errors approach 0 asymptotically for a larger signal. We stress that the exact values of the errors presented in Fig. A1 will vary depending on the redshift, as well as on the background and response considered.
Most spectra in the real data sample with very low SNR are discarded due to fits ending up on the parameter boundaries, while others fail the GOF test presented in section 3.3. However, even after these two cuts, there would remain fits to spectra with SNR so low that the fits do not provide meaningful information. In order to remove these spectra we apply an SNR cut of SNR> 4 in the analysis. This cut is a trade-off between ensuring reliable fits and having enough spectra to analyse. As seen in Fig. A1, the average error due to the imperfect signal at SNR=4 and at a redshift of 1 is, 0.08, 0.07 and 0.23 for L 0,52 , Γ and ε d , respectively.
Here we consider the correlations presented in section 4.2. In Fig. A2 we show the resulting parameter values from fits to 2000 faked spectra where z = 0.1 was assumed. L 0,52 and ε d have been drawn from uniform logarithmic distributions over the entire range of values covered by our model, and Γ was similarly drawn from a uniform linear distribution spanning all available parameter values. It is clear from this figure that there are no intrinsic model dependencies which forces the log L 0,52 − Γ correlation found in section 4.2. The correlation is thus not a model artefact.
However, the simulated data presented in this figure represent data of a far better quality than what is available in our real data sample. Hence, we also consider the same plot using the 2000 previously simulated spectra assuming z = 1. We make a cut at SNR> 4 to remove spectra with a poor signal, as we do in the analysis of real data. The results are presented in Fig. A3. The empty regions in this plot show that we expect some parts of the parameter space to be devoid of fits at certain redshifts. This is a result of the model's physical nature and is not surprising. However, it suggests that the nature of the L 0,52 − Γ correlation may be partially due to selection bias. To investigate this further, we simulate 200 spectra for each data point in the L 0,52 − Γ correlation in Fig. 6, at the redshift of the respective data point. This gives us a simulated data set with as SNR cutoff representative of the real data. In Fig. A4 we show the best-fitting L 0,52 and Γ values from these simulated spectra, plotted together with the L 0,52 −Γ correlation from Fig. 6. We see that the data points in the observed correlation generally do not lie along the boundary to the empty region. This indicates that the log L 0,52 −Γ correlation found in section 4.2 is unlikely to be caused primarily by selection bias.

APPENDIX B: EVALUATION OF GRID RESOLUTION
When fitting with a tabulated model we naturally introduce uncertainties due to the finite resolution of the grid. Specifically, the interpolated spectra between grid points will differ slightly from spectra simulated at the same parameter values. In Fig. B1 we provide an example of an interpolated spectrum plotted together with the corresponding spectrum simulated using our numerical code. The spectra in the figure are generated for the parameter values ε d = 0.125, L 0,52 = 75 and Γ = 375, to show what we may expect from a spectrum that has been maximally interpolated in every parameter. We would naturally expect smaller discrepancies the closer we are to a grid point.
To quantify these uncertainties we simulate a subgrid, consisting of the holes in the parameter grid used to construct DREAM1.2. Each grid point in the subgrid was chosen to have maximal separation between its grid points and the grid points of DREAM1.2. We thus obtain the subgrid The scatter plots show the pairwise relation between all model parameters. The data has been pruned so that only spectra with an SNR>4 are present in the plot. Note that the empty regions are indicative of some possible selection biases in the real data. However, the figure also shows that the log L 0,52 − Γ correlation is not an artefact of the model.  Figure A4. Best-fitting L 0,52 and Γ values (orange points) from fits to spectra simulated from redshifts found in the data of the log L 0,52 − Γ correlation presented in Fig. 6. Only spectra with SNR>4 are present in the plot. The blue points, as well as the black line and grey confidence band, are taken from Fig. 6. spanned by the parameters L 0,52 =0. 25, 0.75, 2.5, 7.5, 25, 75, 150, 250 Γ =125, 175, 225, 275, 325, 375, 425 We simulate spectra from all sets of parameter values of the subgrid and fit DREAM1.2 to them. The spectra were simulated without noise and background in order to isolate the effects of the limited resolution. We construct the relative error of a parameter as presented in appendix A. In Fig. D1 we present the resulting relative error distributions for our three model parameters. The median values and the twosided 1 σ errors for the different parameters are 0.08 +0.08 −0.04 , 0.05 +0.05 −0.03 and 0.08 +0.12 −0.05 for L 0,52 , Γ and ε d , respectively. We stress that these distributions represent the points in the model grid where the errors pertaining to evaluating the model between grid points are the largest. These uncertainties are smaller than the statistical errors on the best-fitting parameters, which are found to be 0.12 +0.16 −0.06 , 0.09 +0.43 −0.06 and 0.20 +0.21 −0.10 for L 0,52 , Γ and ε d , respectively. We conclude that the uncertainty in parameter estimates are not dominated by effects of the limited resolution and that the grid of our model is adequately spaced.
Finally, XSPEC provides an option of whether to use linear or logarithmic interpolation in the table model. When examining the effects of interpolation we do so using one model with linear interpolation and one model with logarithmic interpolation and compare the results. These tests indicate that we obtain a a slightly better precision using linear interpolation, although the discrepancy is not prevalent in the entire parameter space. We have thus chosen to use the linear interpolation for our model.

APPENDIX C: IMPACT OF LAT-LLE DATA
Here we discuss the impact of LAT-LLE data on our analysis. This is relevant since DREAM1.2 has its largest uncertainties at high energies (see section 2.3), and we do not want these uncertainties to bias our fits. To test this we performed all analysis also without the LLE data. The results are summarised in Fig. C1. We find that fitting without the LLE data yields that only 3 spectra change their best-fitting parameter values by more than 1 per cent. Additionally, no additional spectra are accepted or rejected in the GOF test described in section 3.3.
The low impact of the LLE data in this study is mainly a result of two effects. Firstly, when performing the binning described in section 3.2, the SNR of the LLE data are often very low. For comparison, the median SNR of the LLE spectra is 2.3, whereas it is 20.5 for the corresponding NaI spectra. When the difference in SNR is large and there is no strong tension with the model at high energies, the fits are predominately driven by the larger number of counts at low energies.
The second reason why the LLE data do not significantly affect the results is that the vast majority of spectra with bright LLE data are already rejected due to their high flux. Specifically, there are only 21 analysed spectra with LLE data where ε d and L 0,52 do not hit the boundary of the parameter space. Out of these, only 6 spectra have SNR> 1 in the LLE data. Thus, it is not surprising that the general impact of the LLE data is negligible in this study. Therefore, we conclude that the relatively larger model uncertainties at high energies have negligible to no effect on the results.

APPENDIX D: FLUX LIGHT CURVES AND PARAMETER EVOLUTIONS
In Fig. D1 we provide plots equivalent to Fig. 10, for all bursts in our sample with at least 4 accepted time bins.

APPENDIX E: SUMMARY OF FITTED SPECTRA FOR EACH BURST
This section contains table E1, which summarises the number of time bins created by our Bayesian blocks binning. It also shows the number of bins that remain after the GOF test described in 3.3.  Figure C1. Best-fitting parameter values of all 230 spectra that have LLE data, where SNR> 4. The y-axis and x-axis show parameter values from the analysis with and without LLE, respectively. 211 out of the 230 fits hit the boundary in either ε d or L 0,52 , and are thus rejected (note that it is hard to see most points since they are almost all plotted on top of each other). Orange and blue points represent fits corresponding to accepted and rejected fits, respectively. The red line shows the 1:1 correspondence between parameter values from the respective fits.