Observational characteristics of radiation-mediated shocks in photospheric gamma-ray burst emission

Emission from the photosphere in gamma-ray burst (GRB) jets can be substantially affected by subphotospheric energy dissipation, which is typically caused by radiation-mediated shocks (RMSs). We study the observational characteristics of such emission, in particular the spectral signatures. Relevant shock initial conditions are estimated using a simple internal collision framework, which then serve as inputs for an RMS model that generates synthetic photospheric spectra. Within this framework, we find that if the free fireball acceleration starts at $r_0 \sim 10^{10}~$cm, in agreement with hydrodynamical simulations, then the typical spectrum consists of a broad, soft power-law segment with a cutoff at high energies and a hardening in X-rays. The synthetic spectra are generally well fitted with a standard cutoff power-law (CPL) function, as the hardening in X-rays is commonly outside the observable energy range of current detectors. The CPL-fits yield values for the low-energy index, $\alpha$, and the peak energy, $E_{\rm peak}$, that are centered around $\sim -0.8$ and $\sim 220~$keV, respectively, similar to typical observed values. We also identify a non-negligible parameter region for what we call ``optically shallow shocks'': shocks that do not accumulate enough scatterings to reach a steady-state spectrum before decoupling and thereby produce more complex spectra. These occur for optical depths $\tau \lesssim 55 \, u_u^{-2}$, where $u_u = \gamma_u\beta_u$ is the dimensionless specific momentum of the upstream as measured in the shock rest frame.


INTRODUCTION
In a gamma-ray burst (GRB) jet, the photon mean free path is much longer than the kinetic scales of the particles in the plasma.This means that in the regions where the jet is optically thick, sufficiently fast shocks will be mediated by radiation (e.g., Blandford & Payne 1981a,b).The photon spectrum in the downstream of a radiation-mediated shock (RMS) is different compared to the synchrotron spectrum emitted by highenergy particles accelerated in collisionless shocks (see Levinson & Nakar 2020, for a recent review on RMSs).In an RMS, the photons themselves dissipate the incoming particle kinetic energy, which, in the case of an at most mildly relativistic shock with upstream dimensionless specific momentum of u u ≡ β u γ u ≲ a few, results in an extended power-law photon spectrum with a cutoff at high energy.Here, β u is the velocity in units c of the incoming plasma as measured in the shock rest frame and γ u = (1 − β 2 u ) −1/2 .Such spectra resemble those observed in the prompt phase of GRBs (Lundman et al. 2018;Ito et al. 2018;Samuelsson et al. 2022).GRB prompt spectra are commonly fitted with a phenomenological cutoff power-law (CPL) function (e.g., Yu et al. 2016Yu et al. , 2019;;Poolakkil et al. 2021).As the name suggests, the CPL function consists of a power-law segment with an exponential cutoff: N E ∝ E α e −E/Ec , where N E is the specific photon number flux, E is the photon energy, α is the spectral index, and E c is the cutoff energy.The α-distribution found in time-resolved analysis of GRB prompt spectra with the CPL function is a smooth distribution centered at α ∼ −0.8 (e.g., Yu et al. 2016).
The soft (i.e., low) values of α commonly found in catalogues of GRB observations are difficult to obtain in non-dissipative photospheric models, which generate hard observed spectra with α ≳ −0.5 and higher (Deng & Zhang 2014;Acuner et al. 2019).However, dissipation below the photosphere can greatly alter the shape of the photon spectrum, which may not have time to relax to a thermal equilibrium before decoupling from the plasma (Rees & Mészáros 2005).Energy dissipation via RMSs is a natural expectation in the optically thick regions of a GRB jet.Firstly, when the jet is drilling through the star or ejected neutron star material, a high-pressure cocoon forms that keeps the jet collimated.This leads to collimation shocks extending up to, and in some cases even above, the edge of the collimating material (Lazzati et al. 2009;López-Cámara et al. 2013;Gottlieb et al. 2018Gottlieb et al. , 2019Gottlieb et al. , 2021)).Secondly, irregularities within the jet can lead to internal shocks occurring below the photosphere.The irregularities can form either due to mixing between the jet and the surrounding material or stem from a variable central engine (Bromberg et al. 2011;Gottlieb et al. 2019).
As a natural cause for subphotospheric dissipation, RMSs have been thoroughly discussed in the context of GRBs (Eichler 1994;Levinson & Bromberg 2008;Katz et al. 2010;Budnik et al. 2010;Bromberg et al. 2011;Levinson 2012;Keren & Levinson 2014;Beloborodov 2017;Ito et al. 2018;Lundman et al. 2018;Lundman & Beloborodov 2019;Levinson & Nakar 2020;Lundman & Beloborodov 2021;Samuelsson et al. 2022).Unfortunately, the vast difference of interaction scales between the photons and the plasma particles and the high number of photons per charged particle render simulations of RMSs that properly couple the radiation and the plasma extremely time consuming.Therefore, there existed a need for an approximate, faster method that could be used when fitting the model to data.In Samuelsson et al. (2022), we developed such an approximation called the Kompaneets RMS approximation (KRA).The KRA allows us to accurately reproduce spectra from full-scale radiation hydrodynamic simulations using four orders of magnitude less computational time.This massive time reduction made it possible to probe a large enough parameter space such that we could perform a spectral fit using an RMS model to prompt GRB data for the first time.
In this paper, we use the KRA to investigate what typical observational signatures are expected from RMSs in GRB jets.Specifically, we focus on the observed spectrum.To estimate relevant RMS initial conditions in the context of GRBs, we employ a simple internal collision framework (the advantages and disadvantages of this approach are discussed in Section 4).Within this framework, we use the KRA to generate synthetic photospheric spectra, whose typical spectral characteristics can be studied.Furthermore, we compare the synthetic spectra to catalogue distributions of real observations by forward-folding the spectra through the response matrix of the Gamma-ray Burst Monitor (GBM, Meegan et al. 2009) onboard the Fermi Gamma-ray Space Telescope and subsequently fitting them with a CPL function.
The paper is organized as follows.In Section 2, we give a short summary of the KRA.In Section 3, we show that a large parameter range results in shocks where the radiation do not have time to reach steady state, which we call "optically shallow shocks".In Section 4, we describe the methodology: we briefly outline the internal collision framework used in 4.1, the details of which are given in Appendix B, show how to go from the internal collision parameters to the KRA parameters in 4.2, and describe how we obtain a CPL fit from the KRA parameters in 4.3.In Section 5, we present our results.In Section 6, we list the underlying assumptions of the model and in Section 7, we discuss our findings.Finally, we summarize and conclude in Section 8. We will employ the notation Q x = Q/10 x throughout the text.

THE KOMPANEETS RMS APPROXIMATION
In this section, we give a brief summary of the KRA for completeness.We refer the interested reader to Samuelsson et al. (2022) for more details.
The KRA is an approximate but accurate method to simulate RMSs using only a tiny fraction of the computational time required by self-consistent radiation hydrodynamic simulations.One run with Komrad, the simulation code that implements the KRA, takes ∼ 2 minutes to run on a modern laptop.A corresponding run using the full-scale radiation hydrodynamic simulation radshock (Lundman et al. 2018) takes several days on tens of cores.KRA is valid for mildly relativistic or slower shocks (u u ≲ 3) and in regions where radiation pressure dominates over magnetic pressure.In an internal collision framework, the former condition translates into Γ f /Γ s ≲ 36, where Γ f and Γ s are the bulk Lorentz factors of the faster and slower blob, respectively (see Appendix B in Samuelsson et al. 2022).As appropriate in the context of GRBs, the shock is assumed to be photon rich, meaning that the photon number in the downstream is dominated by advection of photons from the upstream (Bromberg et al. 2011).The KRA can model shocks that have finished their dissipation in the optically thick regions.Specifically, we let the RMS propagate for a doubling of the radius, as this corresponds to the shock crossing time of the causally connected region.That means that we do not model shock breakout radiation, which results in more complicated dynamics (e.g., Lundman & Beloborodov 2021).
As photons traverse an RMS, they experience a converging flow.The photons scatter in the velocity gradient and, on average, they gain energy in each scattering.This leads to a bulk Comptonization of the photon spectrum.As both the average relative energy gain per scattering, ⟨∆ϵ/ϵ⟩, and the probability of escape from the RMS are energy independent, a power-law spectrum forms.Here, ϵ is the photon energy and ∆ϵ is the en-ergy change in a scattering, both given in units of electron rest mass energy, m e c 2 .The bulk Comptonization is a Fermi-process, analogous to, e.g., the acceleration of cosmic rays in diffusive shock acceleration.In a planar parallel geometry and with a thermal upstream photon distribution, the evolution of the photon spectrum in the RMS is fully characterized by three parameters: the temperature of the photons far upstream, θ u , the incoming velocity of the upstream as measured in the shock rest frame, β u , and the number density of photons per baryon (here taken to be protons) in the upstream, ñ ≡ n γ /n p .In a photon rich shock, the photon production rate is negligible and ñ is roughly constant across the RMS transition region.
The principle process described above is remarkably similar to that of thermal Comptonization of photons, which occurs if one continuously injects low-energy photons into a region of hot (nonrelativistic) electrons, and let the escape probability be energy independent (e.g., Rybicki & Lightman 1979).Thermal Comptonization is described by the Kompaneets equation, which has the major advantage of being numerically very fast to solve.This realization is the basis of the analogous and approximate description of RMSs introduced in the KRA.In the KRA, we split the, in reality, continuous RMS transition region into three discrete zones: the upstream, the RMS, and the downstream zone.In each zone, the photon distribution is evolved with the Kompaneets equation and photons are advected from the upstream to the downstream, via the RMS, using source terms.Dissipation in the RMS zone is modeled by prescribing the electrons with a high, effective electron temperature, θ r .
Here and throughout the paper, the subscripts u, r, and d when used for the KRA parameters refer to quantities measured in the upstream, RMS, and downstream zones, respectively.When the subscripts u or d are used for the physical RMS parameters, they described quantities measured in the far upstream or downstream, respectively.Furthermore, all temperatures are proper (comoving with respect to each region) and given in units of m e c 2 .
Apart from the effective temperature of the shock, there are two additional parameters that describe the evolution of the photon distribution in the RMS zone in the KRA.The second parameter is R, defined as the ratio of θ r to the upstream temperature, θ u,K , i.e., R = θ r /θ u,K , where the subscript K stands for Komrad.The distinction is necessary as in fact θ u,K ̸ = θ u .The upstream temperatures are different because in a real RMS, adiabatic compression from the upstream to the downstream results in an achromatic energy increase, which is not captured by Komrad.Hence, Komrad re-quires a higher upstream temperature in order to match the low-energy cutoff in the final spectrum (see equation (1) below).
The final free parameter is the Compton y-parameter of the RMS zone, y r .The Compton y-parameter is a measure of the average photon energy gain.In the KRA, it determines the hardness of the power-law segment in the RMS spectrum.A value of y r = 1 indicates a flat νF ν spectrum, where ν is the photon frequency and F ν is the specific flux, and higher values of y r correspond to harder slopes.Specifically, for a given value of θ r , the parameter y r dictates the escape probability of photons from the RMS zone to the downstream zone.A larger value of y r means a lower escape probability and, thus, photons diffuse in the RMS zone for longer, gaining more energy on average, which results in a harder power-law segment.
The effective temperature of the hot electrons in the RMS zone, the temperature ratio between the RMS zone and the cold upstream, and the escape probability from the shock can be chosen such that the photon spectra obtained are remarkably similar to those generated by full-scale radiation hydrodynamic simulations (see Figures 2 and 3 in Samuelsson et al. 2022).To account for the adiabatic compression, we use (Blandford & Payne 1981b) where u d is the specific momentum of the downstream plasma in the shock rest frame.By equating the average relative energy gain per scattering, ⟨∆ϵ/ϵ⟩, in both models, we found in Samuelsson et al. (2022) that where εu (ε d ) indicates the average photon energy in the far upstream (downstream).Here, ξ is a constant and the only free parameter in the conversion.In Samuelsson et al. (2022), we found empirically that ξ = 55 gave good agreement across the whole parameter space.Equation (1) assures that both systems get the same low-energy cutoff.Equation (2) determines the highenergy cutoff of the spectrum, since the maximum energy in the shock, ϵ max , is roughly equal to ⟨∆ϵ/ϵ⟩.The third and final relation is given by equating the average photon energy in the downstream, εd , between the two models.In the case of a real RMS, εd is determined by the shock jump conditions.In Komrad, εd is uniquely determined by the value of y r (given fixed values for θ r and R).Unfortunately, there is no analytic way to relate them to one another.Hence, εd is found in Komrad by numerical integration of the downstream photon distribution.
The shape of the spectrum in the RMS zone is completely determined by the three parameters discussed above.In the downstream zone, the photons continue to interact with the electrons.The downstream electrons have a temperature θ C , where θ C is the Compton temperature defined as the temperature with which there is no energy transfer between the photon and electron population.Thus, there is no energy gain in the downstream for the photon population as a whole.The Compton temperature θ C is not a free parameter in the model: similarly to εd , it is uniquely determined by the parameters R, y r , and θ r .
For the downstream spectrum, the optical depth to the observer at which the shock is initiated, τ i , is an additional parameter that affects the level of thermalization that occurs in the downstream before decoupling.In Samuelsson et al. (2022), we found that the relevant parameter determining the degree of thermalization is actually the product τ i θ r ≡ τ θ, as increasing either one increases the level of thermalization (τ i is a measure of the number of scatterings and θ r is a measure of the energy gain per scattering).1 Indeed, the shape of the spectrum at the photosphere is degenerate as long as τ θ, R, and y r are constant.Note that it is the shape of the photon spectrum that is relevant, since the bulk Lorentz factor and the luminosity of the GRB can shift the spectrum in frequency and flux.If the Lorentz factor is known, as in the case with the internal collision framework below, then the spectra are no longer degenerate and τ i and θ r can be decoupled.

OPTICALLY SHALLOW SHOCKS
Typically, when RMSs are studied they are assumed to have already reached steady state (Levinson & Nakar 2020;Samuelsson et al. 2022).Here, we show that this is not necessarily fulfilled for a wide range of optical depths.In such cases, effects of the shock formation need to be taken into account.
The optical depth required for the radiation inside an RMS to reach steady state, τ ss , can be approximated as the number of scatterings a cold upstream photon (after adiabatic compression) needs in order to reach the maximum energy in the RMS, ϵ max .The average energy increase of a photon as a function of optical depth is dϵ ∼ ⟨∆ϵ⟩ dτ , where ⟨∆ϵ⟩ is the average energy gained in each scattering.Dividing by ϵ and integrating, one gets where εu,K is the average energy of the upstream photons after adiabatic compression.Inside an RMS, ⟨∆ϵ/ϵ⟩ is given by the right-hand-side of equation ( 2).Let us parameterize ln( Here, C is a proportionality constant, which for the results presented in Section 5 is small, ∼ 1.5.However, C could potentially become large in very weak shocks were εd ≳ εu .
With this parameterization, we get As ξ ≈ 55, the equation above indicates that shocks with upstream velocities u u ≲ 1 require optical depths of the order τ ∼ 100 to fully form.Before the shock has reached steady state, traces of the shock formation2 may still be imprinted in the downstream spectrum.However, in this work, u u is commonly larger than unity and steady state is reached quicker, with τ ss ∼ 15.
In terms of the KRA parameters, equation (3) reads where we used ϵ max ∼ ⟨∆ϵ/ϵ⟩ = 4θ r , εu,K = 3θ u,K , as valid for a Wien distribution, and R = θ r /θ u,K .In Komrad, the RMS exists for a doubling of the radius, as this is the time it takes the shock to traverse the causally connected region.This corresponds to a halving of the optical depth.An estimate of the condition required for the RMS to reach steady state is then where the 4/3 in the logarithmic term has been dropped.In Samuelsson et al. (2022), we considered shocks that occurred deep enough such that the RMS had time to reach steady state.Hence, we did not elaborate on the details regarding the shock formation.However, motivated by the fact that a wide range of optical depths can potentially result in optically shallow shocks, we have extended the KRA to more accurately capture these initial phases.The details of the extension is outlined in Appendix A. In short, we have made two changes: 1) the RMS zone now contains an initial thermal population of photons with average energy εd , and 2) the effective electron temperature in the RMS zone is allowed to vary to ensure the energy in the downstream remains constant, as imposed by the shock jump conditions.The latter results in a lower effective electron temperature initially, before it asymptotes to its steady state value θ r as given in equation (2) after around τ ss number of scatterings (see Figure 6).

METHODOLOGY
The aim of this work is to find expected observational characteristics of RMSs in GRB spectra using the KRA.However, relevant ranges for the KRA parameters are not clear a priori and the resulting spectra can appear widely different depending on what initial conditions are used.To estimate what shock initial conditions are expected in the context of GRBs, in this work we employ a simple internal collision framework.This choice is not unique, for instance, one could instead have chosen parameter values based upon the output of numerical simulations of GRB jets.Yet, the internal collision framework is commonly used to characterize typical shock properties and often serves as a benchmark in the literature (e.g., Rees & Mészáros 1994;Kobayashi et al. 1997;Daigne & Mochkovitch 1998, 2000;Lazzati et al. 1999).It has the advantage that all relevant shock properties are analytically obtained from a few higher order parameters.Furthermore, it accounts for possible underlying correlations among the parameters, for instance, the optical depth of a collision should be correlated with the upstream temperature as they are both decreasing with increased distance to the central engine.However, the framework is highly idealised and nature is likely to be more complicated.We stress, therefore, that the intended use of the internal collision framework in this work is solely as a means to obtain relevant RMS initial conditions.
In the next subsection, we give a brief outline of the internal collision framework (all equations used can be found in Appendix B).In subsection 4.2, we show how we go from the higher order internal collision parameters to the KRA parameters.A KRA spectrum is generated from the parameters, forward-folded through the Fermi response, and fitted with a CPL function as explained in subsection 4.3.

Subphotospheric internal collision framework
In this work, we use a simple internal collision framework to estimate the relevant parameter space for RMSs in GRBs.A detailed description of the internal collision framework used is given in Appendix B. Most of the equations in Appendix B can be found in relevant reviews on GRBs such as Mészáros (2006); Pe'er (2015).In this subsection, we give a brief overview.
We consider a variable central engine that launches a slower and a faster fireball with terminal Lorentz factors Γ s and Γ f ≡ ψΓ s , respectively, where ψ > 1 is a free parameter.The two fireballs are launched with a time separation δt, which we parametrize as δt = χr 0 /c.Here, χ is a constant and r 0 is the radius where the free fireball acceleration begins.This parametrization implies that we assume the time separation to be proportional to the light crossing time as the base of the jet.Throughout this work, we use χ = 1, however, in principle χ could be both smaller or larger (Rees & Mészáros 2005).The initial Lorentz factor of the fireballs at the base of the jet, Γ 0 , is also a free parameter, although its value has little effect on the final result.
The speed difference between the two fireballs eventually leads to an internal collision at some radius r coll .Given that ψ ≫ 1, the collision radius is given by r coll = 2Γ 2 s χr 0 .The collision will be subphotospheric as long as r coll is smaller than the photospheric radius of the slower fireball, r ph,s .As we are only interested in subphotospheric collisions, we use the optical depth for the slower fireball at which the collision occurs, τ s , as a free parameter.For a conical outflow, one has τ s = r ph,s /r coll .When the collision occurs, a reverse and a forward RMS are launched into the two fireballs and the shocked plasma is accumulated in a common downstream in between the two fireballs.Somewhat surprising is that the photospheric radius for the common downstream is slightly larger compared to that of the slower fireball due to shock compression of the plasma, even though the downstream has an increased velocity (see equation (B14)).
When the RMS gets closer to the photosphere, photons start to leak out of the shock and the decreasing radiation pressure leads to the formation of two collisionless shocks (Pe'er 2008;Beloborodov 2011;Lundman & Beloborodov 2021).As Komrad cannot model these higher order effects, we restrict ourselves to shocks where τ s > 5 (see Beloborodov 2011).With this choice, we find that the corresponding optical depth for the downstream, τ i , is always larger than 10.After the collision, the energized radiation is contained in the downstream, which travels at a bulk Lorentz factor Γ = Γ s Γ f = Γ s √ ψ.A schematic of the internal shock framework is shown in Figure 1.To minimize the number of free parameter, we fix the values of r 0 and Γ 0 and study only two A schematic of the subphotospheric internal collision framework.Two fireballs are accelerated from the base in a conical outflow, starting at radius r0 and with an initial Lorentz factor Γ0. The slower of the two fireballs has a bulk Lorentz factor Γs and the faster fireball has a Lorentz factor Γ f = ψΓs.The collision occurs at an optical depth (as measured for the slower fireball) τs, with the requirement τs > 5 to avoid the dynamical changes that occur in the RMS at small optical depths.
different jet base configurations (see Section 5 for motivation and further discussion).We title these two scenarios the early acceleration scenario with r 0 = 10 7 cm and Γ 0 = 1.5 and the delayed acceleration scenario with r 0 = 10 10 cm and Γ 0 = 4 .

From the internal collision parameters to the KRA parameters
To arrive at the results presented in Section 5, we start by randomly assigning values to log(Γ s ), log(τ s ), and ψ using normal distributions as given in Table 1.Each parameter has a minimum requirement as Γ s > 5, τ s > 5, and ψ > 2, which collectively removes ∼ 1% of generated parameter triplets.For each valid triplet, we calculate θ u and ñ using equations (B9) and (B10).3With the velocity transformations γ u = (1 − β s β fs )Γ s Γ fs and γ d = (1 − ββ fs )ΓΓ fs , where β fs is the dimensionless velocity of the forward shock as measured in the lab frame and Γ fs = (1 − β 2 fs ) −1/2 , we find u u , u d , and εd using the shock jump conditions (e.g., Beloborodov 2017;Samuelsson et al. 2022).We only include runs where u u < 3 in the final results, as above this limit, relativistic effects such as pair production, Klein-Nishina cross section, and anisotropy start to have an effect (Samuelsson et al. 2022).This removes ≲ 4% of remaining spectra.
With u u , u d , and εd determined, one can obtain the KRA parameters τ θ ≡ τ i θ r and R using equations (1), (2), and (B14).In the two table models that are used to generate the synthetic spectra (one for the early and one for the delayed acceleration scenario, see next subsection), we set bounds as 1 < τ θ < 50 to increase the resolution of τ θ around the most relevant parameter region.Less than 3% of remaining spectra have values of τ θ outside of these limits, and are thus removed.Bounds are similarly set for R and y r in the table model, which are different in the the early and delayed acceleration scenarios, but these have a minor effect.All in all, ∼ 93% of initially generated parameter triplets survive all cuts.
Unfortunately, there is no exact analytical expression that relates y r to εd .Preferably, one would iteratively run Komrad (without adiabatic cooling as this is not accounted for in the jump conditions) and integrate the downstream spectrum until one finds the unique value of y r that, for the current value of R and θ r , gives an average downstream energy of εd (note that εd is independent of τ i when there is no adiabatic cooling).However, this is computationally too costly.Instead, we create a grid with 27 000 spectra (30 values each for θ r , R, and y r , neglecting adiabatic cooling) over the relevant parameter space and save εd .We then find the spectrum in the grid that has the closest values to the current θ r , R, and εd , and use the value of y r from that simulation as the true value.We use different grids for the early and the delayed acceleration scenarios to maximize resolution.This prescription results in a maximum error in y r of 0.09 and 0.05 for the early and delayed acceleration scenario, respectively.We ensure that a recovered value of y r is never on a grid border.

Synthetic spectrum and fit
Having obtained the KRA parameters, we generate a synthetic Komrad spectrum using the Multi-Mission Maximum Likelihood Framework (Vianello et al. 2015, 3ML).In 3ML, we have constructed two table models consisting of 1 000 Komrad spectra in the early acceleration scenario (10 values each for τ θ, R, and y r ) and 5 586 spectra in the delayed acceleration scenario (21, 14, and 19 values for τ θ, R, y r , respectively).The two table models have different ranges of R and y r to maximize resolution.With the given KRA parameter values, 3ML interpolates between the values in the table model to generate the synthetic spectrum.All spectra in the table models have been broadened to account for high latitude emission and different radial contributions to the photospheric emission, as discussed in Appendix C. The broadened spectrum is Doppler boosted into the observer frame using the correct Doppler factor obtained in Appendix C. Also accounting for redshift, each photon is blue-shifted by a factor 5Γ/3(1 + z).We use z = 1 throughout this work.
We account for background and detector effects by forward folding the model spectrum through the response of Fermi, using a response file from an arbitrary burst, here taken to be GRB 150213A.We verified that the specific response does not have an impact on the results by trying several other ones from different GRBs.This means that we are also considering background in our analysis.However, to more clearly observe the features of the signal, we require the synthetic data to have a signal-to-noise ratio (SNR) above 30.The SNR for each generated spectrum varies, and the median is ∼ 45.
Once the synthetic spectrum is obtained, we perform a standard Fermi data analysis (e.g., Goldstein et al. 2012;Yu et al. 2019).In particular, we fit the synthetic data with a CPL function using the Maximum Likelihood framework in 3ML.This analysis yields values for α and E peak , where E peak = (α+2)E c is the peak energy of the CPL function in a νF ν spectrum.

RESULTS
As mentioned in Section 4, in this work we use two different jet base configurations, characterized by their different values of r 0 and Γ 0 .In the first, the free acceleration starts close to the central engine with r 0 = 10 7 cm and Γ 0 = 1.5.Similar conditions have been invoked or found through data analysis in earlier studies and are therefore interesting to investigate (Rees & Mészáros 1994;Mészáros & Rees 2000;Iyyani et al. 2013;Larsson et al. 2015).However, hydrodynamical simulations show that the jets in both long and short GRBs typically remain collimated due to the high pressure of the surrounding material (Lazzati et al. 2009;López-Cámara et al. 2013;Gottlieb et al. 2018Gottlieb et al. , 2019Gottlieb et al. , 2021)).The collimation leads to shocks that suppress the acceleration, effectively moving the base of the conical jet much further out.Motivated by these results, we use r 0 = 10 10 cm and Γ 0 = 4 for our second scenario (from a long GRB simulation in Gottlieb et al. 2019).That the acceleration starts around r 0 = 10 10 cm was also suggested theoretically by Thompson et al. (2007).We call the two scenarios the early and the delayed acceleration scenario, respectively.In both cases, we use χ = 1, however, in principle χ could be both larger or smaller than unity depending on the dynamics at the base of the jet (Rees & Mészáros 2005).

General spectral characteristics
In Figure 2, we show the photospheric spectra produced by Komrad in the early acceleration scenario in red and the delayed acceleration scenario in green for Γ s = 10 1.5 , τ s = 10 1.5 , and ψ = 10.In the early acceleration scenario, the spectrum is hard.This is partly because the photon-to-proton ratio, ñ, increases with r 0 as evident from equation (B10).Thus, for small values of r 0 , fewer photons share the same dissipated energy, leading to a steeper power-law slope (larger values of y r ).The hardness also increases due to the ratio, R, between the upper and lower cutoff shrinking with decreasing r 0 , as indicated by equation (B17).A shorter power-law segment leads to a quicker thermalization.That fewer photons share the available energy also pushes E peak to larger values, as the average photon energy in the comoving frame is higher.
In the delayed acceleration scenario, the slope of the spectrum is softer.The smooth curvature exhibits a double power-law behavior at low energies.The main power-law component (∼ 10−100 keV) is the remainder of the power-law segment generated in the RMS.The low-energy power law (≲ 3 keV) is the Rayleigh-Jeans slope from the initial thermal population, softened due to high-latitude emission and radial contributions to the photosphere, as explained in Appendix C. In the spectrum in Figure 2, the low-energy power law is below the energy range of GBM.
In Samuelsson et al. (2022), we argued that the observed peak energy could correspond either to the lowenergy cutoff or to the high-energy cutoff in the spectrum, as we were agnostic to the value of y r .However, due to the rather small values of ñ obtained when employing the internal collision framework, we find that y r > 1 is essentially always satisfied.This implies a positive slope in a νF ν spectrum, similar to those shown in Figure 2. Thus, the peak energy corresponds to the high-energy cutoff, which is the maximum energy in the spectrum after thermalization and adiabatic cooling.Therefore, in the current model, we always expect a sharp drop in the spectrum above the observed peak.To obtain a high-energy power law, additional mechanisms must be considered such as particle acceleration in subshocks in magnetic or relativistic RMSs (Lundman & Beloborodov 2019; Levinson 2020), or via shear interactions between the jet and the surrounding matter (Ito et al. 2013;Vyas & Pe'er 2022).
The additional break at lower energies in the delayed acceleration spectrum corresponds to the initial cold upstream photons, after having been subjected to thermalization and adiabatic cooling.A lower limit for the additional break is θ ph u,K , defined as the temperature a cold photon would have at the photosphere if one ignored any heating via scatterings.Accounting for adiabatic cooling to the photosphere, one gets θ ph u,K = θ u,K × (6/5)τ , where τ −2/3 i accounts for the idealised adiabatic cooling in the optically thick regime and the coefficient (6/5) appears as to not overestimate the cooling (Pe'er 2008;Beloborodov 2011, see Appendix C).This translates into an expected additional observed break in the spectrum at a temperature above where we used equations (B9), (B11), (B14), (B15).
For the spectrum in the delayed acceleration scenario shown in Figure 2, the break occurs at ∼ 5 keV.An additional break in the low-energy part of the spectrum is interesting as it has been observed in several studies (Strohmayer et al. 1998;Ravasio et al. 2019;Burgess et al. 2020;Gompertz et al. 2023, see also subsection 7.1).
In Figure 3, we show the delayed acceleration scenario spectrum from Figure 2 fitted with a CPL function.The KRA spectrum has been forward-folded through the Fermi response to account for background and detector effects.The top panel shows the count spectrum, while the bottom panel shows the νF ν spectrum.In this case, the hardening at low energies is below the sensitivity window of GBM.Thus, the CPL function provides an adequate fit to the data, as can be seen from the residuals.2, which in turn are obtained from the higher order initial parameter distributions given in Table 1.Top panels show α and bottom panels show E peak , for the early acceleration scenario (left) and the delayed acceleration scenario (right), as also indicated in the panels.The solid dark green lines show the distributions found for the CPL parameters in a time-resolved catalogue of bright Fermi GRBs during the first four years of observations (Yu et al. 2016).It is evident that unless collimation shocks delay the acceleration, the distribution obtained is harder than most observed bursts.All panels are made with χ = 1 and z = 1.

Parameter distributions
In Figure 4, we show histograms of parameter values obtained from 150 CPL function fits to synthetic KRA spectra, using the method described in subsections 4.2 and 4.3.The input distributions for Γ s , τ s , and ψ are given in Table 1 (in Appendix D, we show results for two additional input distributions).The corresponding mean values and standard deviations for the RMS and the KRA parameters are given in Table 2.The generated histograms of α (top) and E peak (bottom) are compared to the corresponding distributions found in a time-resolve analysis of bright Fermi bursts during the first four years of observations (Yu et al. 2016).
Generally, the fits obtained by the CPL function are good, with residuals similar to those shown in the top panel of Figure 3.However, the CPL function cannot quite capture the smooth curvature of the photospheric spectra around the peak.This leads to the CPL function sometimes underestimating the flux at the highest energies, but this depends on the position of the peak energy and the SNR of the input spectrum.In these cases, the Band function may be a better fit.Furthermore, the CPL function cannot account for the hardening of the spectrum that occurs at lower energies in the delayed acceleration scenario.For high values of Γ, this hardening is within the GBM energy range.In these cases, adding an additional power law to the fitting function significantly improves the fit (see subsection 7.1).

Early acceleration scenario
The obtained parameter distributions in the early acceleration scenario are shown on the left-hand-side in Figure 4.It is evident from the figure that in the absence of collimating material that delays the acceleration, the generated spectra are too hard to account for the majority of observed bursts.Furthermore, the E peak distribution is a factor of a few too high compared to the observations.Both of these findings are in agreement with the discussion in subsection 5.1 and with the single spectrum shown for the early acceleration scenario in Figure 2. It can also be deduced from Table 2, which shows that both ñ and R are significantly smaller compared to the delayed acceleration scenario.However, we note that the value of E peak is dependent on the redshift, which is here fixed to z = 1.The power-law slopes in the early acceleration scenario is consistent with roughly a quarter of the time resolved fits in Yu et al. (2016), which have α ≳ −0.6.
The early acceleration parameter distributions in Figure 4 are similar to those found for the non-dissipative photospheric model studied in Acuner et al. (2019).The reason is that the spectral slope is largely determined by the broadening effects discussed in Appendix C and any observable trace of the dissipation is mostly absent.Thus, GRB spectra that are well fitted with a non-dissipative photospheric model are not necessarily without prior dissipation.Dissipation could have been present, but its effect on the final spectrum may be obscured by the geometrical broadening effects.

Delayed acceleration scenario
The results in the delayed acceleration scenario are shown on the right-hand-side in Figure 4.When a collimation shock delays the acceleration, the range of αvalues obtained is more similar to the observed range, with typical values centered around α = −0.8.This shows that dissipative photospheric emission, when accounting for the correct shock physics below the pho-tosphere, can be much softer than what is commonly assumed for photospheric spectra in the literature.
The obtained distribution for the peak energy has a median value of ∼ 220 keV.However, it overestimates the number of bursts with low values of E peak ≲ 100.This is because we artificially increase the SNR to be above 30 in our modelling.In reality, there is an observational bias against bursts with lower values of E peak , which is unaccounted for here.Furthermore, we again note that the E peak distribution is sensitive to the value of the redshift used.

Different initial distributions
The shape of the histograms presented in Figure 4 depend on the choice of the initial parameter distributions shown in Table 1.To check the robustness of our results, we present similar figures for two different input parameter distributions in Appendix D. Furthermore, we show the effect of using the values of ñ and θ u obtained for the faster fireball instead of the slower fireball in equations (B9) and (B10).
It is evident from the figures that the detailed appearance of the histograms changes in these cases.This highlights the importance of the chosen methodology on the final results.However, the conclusions that dissipative photospheric emission can be very soft and closely resemble the spectral shapes that are observed remain robust.Indeed, when using the conditions in the faster instead of the slower fireball, values as low as α < −1.5 are obtained in a few cases.This is mainly due to the increase in the photon-to-proton ratio for the same initial internal collision parameters.

Additional emission channels
In this work, we have studied internal collisions that result in shocks in the optically thick regions.However, the optical depth of an internal collision is set by the properties of the central engine, and there is no a priori reason to believe all collisions would be subphotospheric.Including shocks above the photosphere would alter the α-histograms shown in Figure 4. Specifically, synchrotron emission from fast cooling electrons results in softer values of α ∼ −3/2 (Sari et al. 1998;Burgess et al. 2015).Additionally, photospheric emission without dissipation, or with dissipation but small values of r 0 , give hard values of α > −0.5.Moreover, under certain conditions, such as if the jet becomes transparent during the acceleration phase, values up to α ∼ 0.7 can be reached (Ryde et al. 2017;Parsotan et al. 2018;Acuner et al. 2019).Inspection of the delayed acceleration scenario in Figure 4 shows that we currently under-estimate the frequency of bursts with α-values in these regions.Therefore, it may be that all three channels combined make up the observed histogram.
A model that can properly account for collisions both above and below the photosphere would be an interesting future study.However, we note that the apparent smoothness of the α-distribution may favor a single emission mechanism.It is possible that for a suitable superposition of the histograms presented in Figure 4 together with the ones presented in Appendix D, the current framework could account for the entire α-range.Moreover, allowing r 0 to vary between bursts would also broaden the α-distribution.A more realistic model can also include different observation angles, multiple internal collisions, and/or magnetic fields, as well as allowing for data with lower significance.However, these claims require future study.

Validity of KRA
The KRA assumes negligible magnetic fields in the region where the RMS occurs.As magnetic fields likely play a crucial role in the launching of the jet, this assumption depends on how quickly the initial magnetic field is dissipated.If there exists moderate magnetization across the RMS transition layer, a collisionless subshock can form in the immediate downstream (Lundman & Beloborodov 2019).Charged particles can be accelerated in the sub-shock leading to enhanced photon production via synchrotron emission.This can shift the photon-to-proton ratio ñ to higher values, possibly increasing the viable parameter space of smaller values of r 0 in our model.
In addition, the KRA cannot account for pair production or shock breakout emission.Violent pair production can occur in the downstream of an RMS if the initial relativistic velocity of the upstream is large enough (u u ≳ 3, Lundman et al. 2018).These conditions may be difficult to obtain in an internal shock framework as it requires large differences between the terminal outflow velocities in different parts of the jet.If present, pair creation can greatly alter the spectrum generated in the downstream (Katz et al. 2010;Budnik et al. 2010;Beloborodov 2013;Lundman et al. 2018).
Shock breakout in the context of GRBs is an interesting topic that has seen a lot of recent development, specifically in the context of GRB 170817A (e.g., Bromberg et al. 2018;Gottlieb et al. 2018;Lundman & Beloborodov 2021).When the shock gets close to the photosphere, it broadens and photons start leaking out of the RMS which decreases the radiation pressure until two collisionless subshocks form (Levinson 2012;Lundman & Beloborodov 2021).Therefore, the dynamics are more complicated, which will affect the resulting spectrum.

Double smoothly broken power law
From Figure 2 and the discussion in subsection 5.1, it is evident that the spectrum in the delayed acceleration scenario has more complexity than a simple CPL function.To probe this characteristic, we employ a double smoothly broken power-law (2SBPL) function as defined in Ravasio et al. (2018).The 2SBPL is an empirical function with two low-energy power laws with indices α 1 and α 2 , smoothly connected at a break energy E break .The low-energy part of the spectrum is smoothly connected to a high-energy power law with index β at the peak energy E 2SBPL peak .The smoothness of the two breaks are determined by two smoothness parameters, which we fix to n 1 = 5.38 and n 2 = 2.69 similar to Ravasio et al. (2018) (see also Kaneko et al. 2006).
As in subsection 5.2, we fit 150 synthetic photospheric spectra with the 2SBPL function, with initial parameter distributions again given by Table 1.The distribution for the low low-energy spectral index, α 1 , and the high low-energy spectral index, α 2 , are shown in the left panel in Figure 5.In the right panel, we show the distribution for the ratio of peak energies to break energies.
By performing CPL fits to the same 150 input photospheric spectra and using the Akaike Information Criterion (AIC), we can evaluate which model is preferred.The AIC accounts for the number of free parameters in a model and a lower AIC value indicates a better model fit to the data.We find that ∼ 20% of cases have ∆AIC ≡ AIC 2SBPL − AIC CPL < 0 with ∼ 10% of total fits having ∆AIC < −4.As many of the generated synthetic spectra are well fitted with the CPL function within the GBM energy window, there is often no need for the extra degrees of freedom in the 2SBPL function and, thus, ∆AIC > 0 in a majority of cases.Still, this result shows that the 2SBPL function provides a better fit than the CPL function for a large fraction of photospheric spectra.The spectra that with ∆AIC < 0 have higher values of Γ on average, leading to a ∼ 3 times higher median value of the observed radiation luminosity, L γ (calculated using equation (11) below), as compared to the whole fitted sample.Thus, we would expect the 2SBPL to be a better fit in bright bursts.
The distributions of the parameter values obtained when fitting a 2SBPL function to KRA spectra have large dispersions.It is, however, remarkable that the obtained values are quite similar to those expected in marginally fast cooling synchrotron models (α 1 = −2/3, α 2 = −3/2, E 2SBPL peak /E break ∼ few).This sim- ilarity is especially clear in the case of larger Lorentz factors, as evident from the middle panels in Figure 12.
Obtaining values close to the values expected from synchrotron theory is not due to the KRA itself, i.e., there is no reason intrinsic to the model for why this would occur.Rather, it is an effect of the limited energy windows of the detectors and the curvature of the spectrum.If a given spectrum has more complexity than a single power law at low energies, then it is likely that the 2SBPL function will provide a better fit due to its higher flexibility.Furthermore, if the power-law index found when fitting a CPL function to the spectrum is α ∼ −1, then the two indices found in a 2SBPL fit are commonly going to satisfy α 2 < α < α 1 .Additionally, both the break energy and the peak energy are usually found within the limited detector energy window, leading to a moderate ratio between the two.The closeness of the two breaks can also be explained by the fact that the curvature of the 2SBPL function around E 2SBPL peak cannot quite capture the smooth curvature around the peak in the input spectrum.Thus, E break ≲ E 2SBPL peak to account for the hardening below the peak.In many cases, therefore, E break does not correspond to the break in X-rays, which is visible around 5 keV in Figure 3.In Figure 5, these cases correspond to the overrepresentation of fits where α 2 ∼ −2, which corresponds to a flat power law in a νF ν spectrum.
The above discussion illustrates that one has to be cautious when drawing conclusions regarding the physics based on fits of empirical functions alone, in particular over limited energy ranges.Identifying spectral indices close to −0.67 and −1.5 is not sufficient to conclude that the emission is synchrotron radiation.Such a spectrum can be produced by a dissipative photosphere, as evident from Figure 5, or by another emission mechanism that produces a curved spectrum within the energy range of the detectors.
To break this model degeneracy, complementary and independent properties of the emission needs to be assessed.For instance, Ravasio et al. (2019) fitted a 2SBPL to the time-resolved data of the large, prompt flare in GRB 160821A and found averaged power law indices of −0.74 ± 0.15 and −1.51 ± 0.17.For this burst, further independent indications corroborate that the emission is indeed synchrotron.First, the flare has a large polarisation degree (Sharma et al. 2019;Gill et al. 2020, although, see Lazzati et al. (2004)) and second, direct fits with an actual synchrotron model yield good agreement to the data and reveal an accelerated electron distribution consistent with the theory of particle acceleration in weakly magnetized flows (Ryde et al. 2022).Additional indirect arguments in favor of synchrotron origin for this burst is that the emission occurs late and has a long duration (Li et al. 2021).
Another way to break the degeneracy is to obtain early observations outside of the γ-ray band.A KRA spectrum is narrower than that predicted from synchrotron models.For instance, the asymptotic power-law slope at low energies in a broadened photospheric spectrum is much harder (α ∼ 0.4, Beloborodov 2010;Acuner et al. 2019) than that expected from synchrotron emission (α ∼ −2/3, Sari et al. 1998).Thus, optical observation have the power to differentiate between the models.Such a study was conducted in Oganesyan et al. (2019).Unfortunately, current optical observations of GRBs always occur quite late (typically ∼ 100 s after trigger).A few early optical observations may become available with the launch of the SVOM satellite in mid-2023 (Wei et al. 2016).Similarly, early high-energy data is a good way to distinguish between the models, as the KRA spectra drop rapidly in flux above the peak.

Differences between short and long GRBs
In this work, we have found that the low-energy index and the peak energy correlate with the injection radius, r 0 .Furthermore, the observed variability time is also dependent on r 0 .It is therefore interesting to speculate whether the difference in these properties between short and long GRBs could be due to differences in their initial collimation.Gottlieb et al. (2019Gottlieb et al. ( , 2021) ) found that in simulations of short GRBs, the collimation shock moves outwards and often breaks out of the surrounding neutron star ejecta.Hence, it may be difficult to prescribe a single value to r 0 .However, if we for simplicity assume that r 0 is an order of magnitude smaller on average in short GRBs compared to long ones, then one can estimate (see equations ( 9) and ( 12) below) that E peak should be a factor 1.8 times higher and the variability time should be 10 times smaller in short GRBs compared to long ones.Moreover, by generating 150 spectra using r 0 = 10 9 (still with Γ 0 = 4) and fitting them with the CPL function, we find a mean value of ⟨α⟩ = −0.55,as compared to the mean value of ⟨α⟩ = −0.76found for the distribution in Figure 4 when r 0 = 10 10 cm.All these differences are close to the values observed.Golkhou & Butler (2014) found the observed variability time to be a factor 12.5 smaller in short GRBs compared to long ones.In a time-integrated analysis of the first 10 years of Fermi data, Poolakkil et al. (2021) found that E peak was a factor of 2.6 higher in short GRBs compared to long ones on average.In the same work, Poolakkil et al. (2021) also found ⟨α⟩ = −0.59± 0.5 for short GRBs.This can be compared to the mean value ⟨α⟩ = −0.80 ± 0.311 found in the time-resolved analysis of Yu et al. (2016), where 80 of the 81 bursts analysed were long GRBs.

Additional predictions within the internal shock framework
Within the internal shock framework employed, one can make many additional predictions regarding the emission.In this subsection, we give estimates of the expected peak energy, efficiency, isotropic gamma-ray luminosity, and observed variability time.However, we stress that since these predictions are based on the idealised internal shock framework, they should be evaluated with caution.
Peak energy.In the downstream of the RMS, the photon distribution will interact with thermal electrons at temperature θ C .In thermal Comptonization, the average relative energy gain per scattering is ⟨∆ϵ/ϵ⟩ = 4θ−ϵ.Thus, high-energy photons with ϵ > 4θ lose their energy more quickly compared to the energy gained by low-energy photons.The energy loss halts once the energy approaches 4θ C , which is always greater than the mean photon energy of the distribution.Therefore, one can use εph d to get a lower-limit for the high-energy cutoff in the spectrum.Here, εph d is the average photon energy in the spectrum at the photosphere.Furthermore, as almost all spectra have a positive slope in their νF ν spectrum (see subsection 5.1), εph d should be a good lower limit for the observed peak energy.The dissipated energy per proton in the RMS is approximately (γ u − 1)m p c 2 .Hence, the average photon energy at the photosphere can be estimated as where the term (6/5)τ −2/3 i accounts for the adiabatic cooling of the distribution.With a Doppler factor of 5Γ/3, accounting for redshift, and using equations (B12) and (B14), one gets an estimated lower limit for the peak energy where the numerical factor 5.83/3 transforms the average photon energy to the peak energy in a νF ν spectrum (Ryde et al. 2017), we used (γ u − 1) ≈ ( √ ψ − 1) 2 /2 √ ψ obtained from equation (B12), we used equation (B14), and we only show the leading parameter dependence of ψ in the final line.The estimate above explains why the fitted values for E peak are higher on average when r 0 is 10 7 cm as compared to 10 10 cm (see Figure 4).
We compared E ll peak to the values of the peak energy found for the fits in the histograms in Figure 4.In the early acceleration scenario, the values of E peak are very close to E ll peak indicating almost complete thermalization of the high-energy photons.In the delayed acceleration scenario, E ll peak is commonly a factor of a few lower than the peak energy, indicating incomplete thermalization.
For instance, the fit shown in Figure 3 yielded E peak = 150 keV compared to ∼ 62 keV obtained using equation ( 9).Efficiency.As explained in Appendix B, both fireballs are assumed to have similar total masses.With the mass in each of the two fireballs denoted by M = N p m p , the total energy emitted from the central engine is E tot = (ψ + 1)Γ s N p m p c 2 .The total radiated energy in the central engine frame is E γ ≈ 2(5Γ/3)N γ εph d m e c 2 where N γ is the total number of photons in each fireball, 5Γ/3 is the Doppler boost, and an additional factor 2 comes from the fact that there are two fireballs.With N γ = N p ñ, we get an estimated radiation efficiency of where we have inserted εph d from equation ( 8).Using the more accurate value of εph d determined from the shock jump conditions, the mean efficiency is ∼ 6%, which agrees quite well with the estimate in equation ( 10).As evident from the equation above, the efficiency is almost entirely determined by the optical depth of the RMS, meaning that photospheric emission can be very efficient if the dissipation occurs close to the photosphere.
Isotropic gamma-ray luminosity.The total proton number is the same in both fireballs.In the slow fireball, it can be estimated as N p = E s /Γ s m p c 2 ≈ L s (l s /c)/Γ s m p c 2 , where E s it the total emitted energy in the slow fireball and l s is the observed width of the slow fireball.The RMS traverses the upstream in a doubling of the radius, which implies that l s ∼ δt c = χr 0 .The source frame γ-ray luminosity can be estimated as L γ ∼ E γ /(t var /(1 + z)), where t var is the observed variability time, which divided by (1 + z) is the source variability time.Thus, where we inserted L s and t var from equations (B6) and ( 12) and only show the leading dependence on ψ in the last line.
Variability time.In our delayed acceleration scenario, the values for the photospheric radius obtained in this work are of the order ∼ 10 15 cm, which is much larger than what is commonly discussed for photospheric models (∼ 10 11 −10 12 cm).The large radii are due in part to the high isotropic equivalent luminosities required to keep the shocks subphotospheric (∼ 10 54 erg s −1 ) and in part due to the low Lorentz factors (∼ 100).These two factors alone increase the photospheric radius to 6 × 10 14 cm (e.g., Pe'er 2015).Additionally, the compression in the downstream of the RMS further pushes the photospheric radius outwards by a factor of a few.
The large photospheric radius directly solves one problem with photospheric models in GRBs, which is the short variability times predicted.For an instantaneous flash in a thin shell occurring at a radius r, the emitted light reaches the observer over a timescale t var = (1 + z)r/2Γ 2 c.With a canonical photospheric radius at ∼ 10 12 cm and Lorentz factor of Γ = 100, one expects millisecond variability, compared to the commonly observed (minimum) variability time of ∼ 2.5 s in long GRBs (Golkhou & Butler 2014).However, in the delayed acceleration scenario, one gets where we used equation (B14) and set z = 1.Note that in a realistic scenario with multiple collision below the photosphere, variability on different timescales can be observed simultaneously in the light curve.
In Rees & Mészáros (2005) it was suggested that χ = θ j , where θ j is the half-opening angle of the jet.With θ j = 5 • , the luminosity, photospheric radius, and variability time would all decrease by a bit more than an order of magnitude.

SUMMARY AND CONCLUSION
We have studied what observational characteristics are expected in dissipative photospheric models of GRBs, when the dissipation is in the form of an RMS.The RMS was modeled using the KRA developed in Samuelsson et al. (2022), which allowed us to investigate the spectral properties quantitatively due to its low computational cost.To estimate relevant shock initial conditions in the context of GRBs, we employed a simple internal collision framework (Figure 1).We found that when the free fireball acceleration starts close to the central engine, the spectral slope was hard.However, in hydrodynamical jet simulations, one commonly finds that the acceleration is delayed due to collimation, resulting from high-pressure material surrounding the jet (Lazzati et al. 2009;López-Cámara et al. 2013;Gottlieb et al. 2018Gottlieb et al. , 2019Gottlieb et al. , 2021)).In this case, the obtained spectrum was softer and exhibited an additional break in X-rays (Figure 2).The spectrum was well fitted over the limited GBM energy range with a phenomenological CPL function (Figure 3).
We continued by generating 150 synthetic photospheric spectra as they would appear in the Fermi GBM detector and fitting them with the CPL function.The derived histograms for the low-energy photon index, α, and peak energy, E peak , were compared to the catalogued distributions found in a time-resolved analysis of bright Fermi burst (Yu et al. 2016).We found that when the free fireball acceleration started close to the central engine, the spectra were generally hard with α ≳ −0.75, inconsistent with a majority of observed GRBs.In the scenario where the acceleration was delayed, the average value for the power-law index was α ∼ −0.8 with a range of values consistent with the dominant fraction of catalogued bursts.This showed that photospheric emission with correctly modeled prior dissipation can be much softer than what is commonly argued.The obtained values of E peak were similarly in general agreement with the observations (Figure 4).
We also fitted 150 synthetic spectra with a 2SBPL function and found that such a fit was strongly preferred over the CPL function in ∼ 10% of cases.The bursts where the 2SBPL had lower AIC value than the CPL function had a ∼ 3 times higher median gamma-ray luminosity compared to the sample as a whole.The average values obtained for the low-energy indices α 1 and α 2 , as well as the ratio of peak energy to break energy are quite similar to those expected in marginally fast cooling synchrotron emission models (Figure 5).As this is not an intrinsic property of our dissipative photospheric model, we argued that this is an effect of the limited energy windows of the detectors and the smooth curvature of the spectrum.Thus, one must be cautious when drawing conclusions about the underlying physics based on empirical spectral fits, as different physical models can yield similar fitted parameter values.We listed several ways this degeneracy could be broken, including polarization measurements, fitting physical models directly to the data, and early optical and high-energy observations.Finally, we found that for a large range of optical depths, the radiation in the RMS does not have time to reach steady state before the dissipation ends.We called these shocks "optically shallow shocks" and found in equation ( 4) that they occur for optical depths as deep as τ ∼ 55 u −2 u .Optically shallow shocks are interesting, as signatures of the shock formation may still be imprinted in the observed spectrum.
In conclusion, we find that photospheric emission that has been energized by a subphotospheric internal collision is spectrally consistent with the majority of observed GRBs.Additionally, this type of photospheric emission can explain various observational features such as an additional break in X-rays, the observed variability time, and possibly some of the main differences between short and long GRBs.This requires that the free fireball acceleration is delayed until ∼ 10 10 cm, primarily to increase the photon-to-proton ratio (see also Table 2).The main caveats relate to the internal collision framework used and are the low efficiency (∼ 6%) and the high isotropic equivalent luminosity required (L ∼ 10 54 erg s −1 ).The luminosity requirement can be relaxed if the photon-to-proton ratio can be otherwise increased, e.g., by additional photon production at a collisionless subshock.

A. EXTENDING THE KRA TO TREAT SHOCK FORMATION
In the collision of two blobs, the central region between the blobs is initially compressed adiabatically.The compression increases the temperature and, in the opaque systems below the photosphere, the radiation pressure.When the downstream radiation pressure becomes comparable to the ram pressure of the incoming baryons, which roughly occurs when the average photon energy in the compressed region equals εd , two shocks form (Beloborodov 2017).Hence, the shocks initially contains a thermal distribution of photons with average energy εd .
The shock jump conditions dictate that the average energy of photons advected into the downstream should be εd (ignoring adiabatic cooling, which is small over the lifetime of the of the shock and does not change the reasoning).Since the escape probability from the RMS to the downstream is energy independent, the average energy in the RMS zone must also equal εd .However, this requires ⟨∆ϵ/ϵ⟩ to vary during the shock formation.Thus, in the KRA we want a varying effective electron temperature, θ v r , that keeps εr constant (which in turn keeps εd constant) and which equals θ r as given by equation (2) when the radiation in the RMS reaches steady state.
The average photon energy in the RMS zone can change due to Compton scatterings with the hot electrons and photons entering and leaving the shock region.The fraction of photons that leave the shock between an optical depth of τ and τ + dτ is given by (4θ r /y r ) dτ when the radiation is in steady state (see Samuelsson et al. 2022, for details).Assuming the number of photons that pass through the shock is the same during the shock formation as during the later steady state phases, this applies to the early stages as well.The expression for the average energy change in the RMS zone in the KRA is then given by where θ C,r is the Compton temperature of the photons in the RMS zone.Using the equation above, we wish to find the shock temperature θ v r that satisfies dεr dτ = 0. To ensure stability, we replace εr with εd in the last term in the last parenthesis, since this is the average photon energy that should be injected into the downstream.Note that the average energy of the photons injected into the downstream in the simulation is still εr ; this prescription is only used in the calculation of θ v r .Furthermore, with a fine enough radial resolution, we checked that this choice has no effect since at all times εr = εd .Solving for θ v r , one gets In Figure 6, we show how θ v r varies during the shock formation.Initially, it is lower than the asymptotic value θ r .As photons are heated, the photon distribution changes and θ v r grows.In terms of a true RMS, this can be thought of as a steepening of the velocity gradient when the shock is initially launched.The growth saturates when the influx of cold photons starts to dominate over the energy increase due to scatterings, resulting in a decrease of the Compton temperature, and, hence, of θ v r .Finally, θ v r converges to θ r .The qualitative behavior is independent on the choice of initial parameters.
In Figure 7, we show how the RMS and the downstream photon distribution evolves in a simulation by plotting ten snapshot spectra of each.From the Figure, one can see how the RMS spectrum appears before it reaches steady state.The downstream initially shows clear signs of the shock formation, but for this simulation they are washed out before the downstream reaches the photosphere.Indeed, this is the case for most of the spectra in Figure 4: less than 10% of the spectra in Figure 4 qualify as optically shallow shocks, as calculated by equation (6).

B. DETAILS REGARDING THE SUBPHOTOSPHERIC INTERNAL COLLISION FRAMEWORK
We consider a variable central engine that launches a slower and a faster fireball with a time separation δt.The two fireballs are assumed to have similar comoving densities, ρ, and similar total masses, M , but different isotropic equivalent luminosities, L s and L f , which results in different terminal Lorentz factors, Γ s and Γ f .For a collision to occur, we require Γ f /Γ s ≡ ψ > 1.That the two fireballs have similar proper densities implies that the ratio of the isotropic equivalent luminosities of the two fireballs satisfies L f /L s = ψ 2 .
Initially, the fireballs accelerate until most internal energy has been converted into bulk kinetic energy.This occurs at the saturation radius where the subscript i indicates fireball i. Above, r 0 and Γ 0 are the radius from the central engine and the Lorentz r initially grows as the Compton temperature in the shock grows.It tends to its steady state value θr, which is indicated by the horizontal dotted line.As expected, the average photon energy in the RMS, εr, remains at all times constant and equal to εd , which is indicated by the horizontal dashed line.The vertical dot-dashed line shows our estimate of τss calculated using equation ( 5).It slightly underestimates when equilibrium is reached, which is expected since equation ( 5) assumes a constant energy gain per scattering of 4θr, while the true energy gain is equal to θ v r , which is lower than θr initially.Adiabatic cooling is not included.Parameter values are τ θ = 5, θr = 0.1, R = 500, yr = 1.5.
factor at the base of the conical outflow.These could in principle be slightly different for the two fireballs but here we assume that they are the same.The radius at which the two shells collide is given by where we have parameterized δt = χr 0 /c, that is, we assume that the time separation between the two fireballs is related to the light crossing time at the radius where the free acceleration begins.In this work, we use χ = 1 in both the early and delayed acceleration scenario.However, we include it in the equations to show the parameter dependence.The collision will be subphotospheric if r coll is smaller than the photospheric radius of the slower shell.In a conical outflow where the density drops as r −2 , the optical depth toward the observer at the collision radius for a fluid element in the slower shell is given by τ s = r ph,s /r coll .As we are only interested in collisions that occur below the photosphere, we use τ s as a free parameter in our collision model.Apart from τ s , we use ψ and Γ s as free parameters.

(B6)
The comoving temperature of the plasma in fireball i at the collision radius, assuming the collision occurs above the saturation radii of both fireballs, is given by where θ 0,i is the temperature of fireball i at the base of the conical outflow.Assuming a thermodynamic equilibrium at the base, the initial temperature can be found through conservation of momentum flux as where k is the Boltzmann constant, a is the radiation constant, and β 0 = 1 − 1/Γ 2 0 .Inserting this back into equation (B7) and using equation (B6), one gets The last line in the equation above indicates that the value and parameter dependence for θ u,f is the same as that for θ u,f , apart from a factor ψ 1/6 .The photon-to-proton ratio, ñ, can be estimated assuming that photons and pairs are in thermodynamic equilibrium at the base of the jet, that all pairs have recombined before the collision radius, and that no new photons have been produced.In that case, one gets (Bromberg et al. 2011  Initially, the shock contains a Wien thermal photon distribution with mean energy εd .As the shock evolves, the initial peak moves to higher energies and broadens, while new cold photons enter the shock from the upstream.At τ ≈ 33 (snapshot seven), steady state has been reached and the final spectra are overlapping.This corresponds well with the estimate shown in Figure 6.Right: Similar to the left panel but for the downstream spectrum.This time the spectra range from τ = 50 all the way up until the photosphere at τ = 1 (from darker to brighter green).Note that the spectrum for τ = 50 is not visible since the downstream does not contain any photons at the start of the simulation.For this specific choice of parameters, the optical depth was large enough to erase all signatures of the shock formation in the final spectrum.Adiabatic cooling and the broadening effects explained in Appendix C are not included.Parameter values are the same as in Figure 6.At the collision radius, the adiabatic compression between the fireballs launches a forward shock that propagates into the slower fireball and a reverse shock that propagates into the faster fireball.From equations (B9) and (B10), we see that the upstream conditions for these two shocks are slightly different in the two fireballs.However, these differences may easily be outweighed by the uncertainties in r 0 and Γ 0 .For simplicity, in this work we assume the upstream conditions in both shocks are similar.For the results presented in Section 5, we used the values obtained for the forward shock, i.e., the upstream conditions in the slow fireball.However, in the main text, the subscript s has been dropped.Using the upstream conditions in the faster fireball instead, the photon-to-proton ratio increases, which in turn softens the power-law and leads to to lower values of both α and E peak .The results obtained when using θ u,f and ñf are shown in Appendix D.
The Lorentz factor of the causally connected downstream in between the two shocks, given similar proper densities of the two fireballs, is given by (Kobayashi et al. 1997;Samuelsson et al. 2022) The assumption of similar proper densities in the two shells assures that the upstream relativistic velocity, u u , is the same for the reverse and the forward shock.To determine the upstream relativistic velocity one must solve the shock jump conditions, which is done for all the results presented in main text.However, it is useful with analytical estimates, which we give below.The Lorentz factor of the incoming upstream can be estimated by where Γ fs is the Lorentz factor of the forward shock as measured in the lab frame, and we used Γ fs ≈ Γ for the last approximation.This gives where f (ψ) = 1 2 (1 − 1/ψ) satisfies 1/4 < f (ψ) < 1/2, as long as ψ > 2. This estimate shows that the dimensionless momentum is expected to be mildly relativistic with u u ≲ 3 for a wide range of ψ.
Although the Lorentz factor of the downstream is larger than the Lorentz factor of the slower fireball, the optical depth of a downstream fluid element toward the observer is actually slightly larger due to compression of the plasma (Levinson & Nakar 2020).(This is only true if one ignores diffusion out of the downstream.Indeed, diffusion in a thin shell might make photons decouple from the downstream before reaching the photosphere (Ruffini et al. 2013;Bégué et al. 2013).Therefore, an implicit assumption of our model is that the specific collision considered is embedded in a jet where many such collisions occur.We note, however, that the diffusion time becomes comparable to the dynamical time close to the photosphere and that all our shocks have finished dissipation before this happens.Indeed, the shape of the spectrum is determined at early times before diffusion becomes relevant, as the number of scatterings decreases with density.Including escape via diffusion would only slightly decrease the amount of adiabatic cooling of the spectrum.)The optical depth for the downstream at the collision radius is given by where the compression is given by the ratio u u /u d , coming from the conservation of baryon number, and the factor 1/ √ ψ accounts for the increased Doppler boost.The function f is monotonically decreasing with ψ and depends only weakly on the other parameters.Its numerical value is f ≈ 4 for ψ = 2 and f ≈ 2.6 for ψ = 30.The compression also implies that the photospheric radius moves outwards by a factor f , as compared to the photospheric radius of the slower fireball.By looking at equation (B5), we see that the photosphere can be as far out at r ph ≳ 10 15 cm.Some implications of this is discussed in subsection 7.3.
We continue with analytical estimates of the KRA parameters.The upstream temperature is given by where θ u is given in (B9).The steady state value of the effective electron temperature can be estimated as where we used εd /ε u = 100 in the last equality, neglecting their parameter dependence for simplicity as it only enters logarithmically and therefore has a small effect.
The ratio R = θ r /θ u,K is then estimated as and τ θ = τ i θ r is given by (B18) Although we cannot estimate the parameter dependence on y r analytically, its numerical value is y r ≈ 1.6 for the values used in this subsection.

C. HIGH-LATITUDE AND RADIAL CONTRIBUTIONS TO THE PHOTOSPHERIC EMISSION
Three ingredients are needed in order to post-process the transition to transparency: the local comoving spectral shape at the photosphere, the distribution of the last scattering radius and angle, and the radial cooling of radiation as the jet transitions to transparency.The spectral shape is obtained from Komrad.The probability density for the last scattering as a function of radius, r, and comoving propagation direction, µ = cos θ, where θ is the angle between the comoving radial direction and the line-of-sight, is given by equation (B16) in Beloborodov (2011).The radial cooling profile is obtained from the full scale transfer simulation code radshock (Lundman et al. 2018).The broadened spectrum is then obtained by integrating the local spectrum over radius (or optical depth) and angle.This prescription assumes that there is no structure to the jet, i.e., that the comoving spectrum does not change with angle within 1/Γ to the line of sight of the observer.
Assume that adiabatic cooling is described by a cooling function ϕ(τ ), such that the energy of a photon trapped in its fluid element as measured by a comoving observed at optical depth τ can be written as where ϵ ph is the corresponding photon energy as measured by a comoving observer at the photosphere on the line-of-sight.This prescription neglects the evolution of the spectrum close to the photosphere due to thermalization, i.e., it only accounts for adiabatic cooling.However, we expect the thermalization to be inefficient at small optical depths where the photons are released.Indeed, Beloborodov (2011) found that ∼ 80% of photons are released at optical depths τ < 3. Calculating a Compton y-parameter defined as y = released and the low-energy part of the observed spectrum should not be affected by this prescription.Since all photons are affected similarly by adiabatic cooling, i.e., ϕ is not a function of ϵ, the (logarithmic) spectral shape is constant as the photons cool.Specifically, we have that where ϵ i = ϕ(τ )ϵ ph,i as given by equation (C19), and N ϵ = dN/dϵ and N ϵ ph = dN ph /dϵ ph are the photon number spectra at optical depth τ and the photosphere, respectively.Since this is true for any ϵ 1 and ϵ 2 , one gets that As photons make their last scattering over a range of radii where the average photon energies are different, the observer receives a spectrum that is broader compared to the comoving photospheric spectrum.Furthermore, the observer receives photons from different scattering angles, which have different Doppler boosts.The observed photon energy is related to the emitted one by ϵ obs = D(µ)ϵ(τ ) = D(µ)ϕ(τ )ϵ ph , where D(µ) = Γ(1 + βµ) ≈ Γ(1 + µ).As we are in the coasting phase, Γ is constant.The observed spectrum will be a superposition of spectra emitted at different optical depths and angles.Let f (τ, µ)dτ dµ be the probability that an observed photon was emitted at an optical depth between τ and τ + dτ and at an angle between µ and µ + dµ.The observed spectrum consisting of these specific photons has the same (logarithmic) spectral shape as the emitted comoving spectrum, just Doppler boosted by a factor D(µ).Using a similar argument to that of equation (C20), one gets N ϵ obs dτ dµ = f (τ, µ)dτ dµ N ϵ ph /D(µ)ϕ(τ ), where 1/D(µ)ϕ(τ ) = ∆ϵ ph /∆ϵ obs .The observed spectrum is obtained by integrating over optical depth and angle as where N ϵ ph = dN ph /ϵ ph is everywhere evaluated at ϵ ph = ϵ obs /D(µ)ϕ(τ ) (note that ϵ obs in equation ( C21) is a single, fixed observed photon energy).Thus, it cannot be taken outside of the integrals since it is a function of τ and µ.By integrating over energy, one can confirm that the prescription conserves particle number: where this time, N ϵ ph dϵ ph = N ph can be evaluated first and taken outside of the integral, since the total photon number at the photosphere is not a function of τ or µ.Equation (C22) also uses the fact that f (τ, µ)dτ dµ = 1.
Informed by the output of the full radiative transfer simulation code radshock, we find that ϕ(τ ) = C(τ 2/3 + 0.2) with C being a constant, reproduces the spectrum very well.Using the fact that ϕ(τ ) = 1 for τ = 1 as required by equation (C19), we find that C = 1/1.2.Thus, This expression is similar to the one found for the first moment of radiation intensity, I 0 , in Beloborodov (2011).However, there is a small discrepancy, which may stem from the assumption of isotropy in Beloborodov (2011).
The expression for how the photon energy evolves with optical depth due to adiabatic cooling given in equation (C23) is not compatible with the prescription used in Samuelsson et al. (2022).In Samuelsson et al. (2022), we numerically stopped the adiabatic cooling in Komrad at an optical depth of τ = 3.This was based on Figure 3 in Beloborodov (2011), where it is shown that the total adiabatic cooling of the photon spectrum as calculated by radiative transfer is equal to an ideal cooling of the photon spectrum (as τ −2/3 ) if the ideal cooling is stopped at τ = 3.However, in Figure 3 in Beloborodov (2011), the adiabatic cooling shown is for the luminosity in the lab frame.Within the broadening prescription described in this Appendix, the output from Komrad represents the comoving photon distribution at τ = 1.The cooling history (and future cooling for the photons that decouple at τ < 1) of these photons are recreated via equation (C19), and integrating over the probability distribution f (τ, µ)dτ dµ generates the superposition that is the observed spectrum.Thus, it is more accurate to numerically stop the ideal cooling at τ = 1.2 3/2 = 1.31, as evident from equation (C23), which has been incorporated throughout this paper. .Two spectra before (solid lines) and after (dashed lines) the broadening effects have been added.The solid lines correspond to comoving spectra at the photosphere and the dashed lines to the shapes of the observed spectra, but normalized to have the same average photon energy as the photospheric spectra (see the text for details).The red spectrum has a smaller ratio between the peaks (R = 50) than the green (R = 1000).Both initial spectra have been normalized such that θu,K = 1.It is evident that the effect of the broadening scheme is most important for the narrower spectrum.Other parameter values are τ θ = 5 and yr = 1.5.
The probability density function f (τ, µ) is found by transforming equation (B16) of Beloborodov (2011)  With ϕ(τ ) and f (τ, µ) given by equations (C23) and (C24), we find that the double integral in the last line of equation (C25) equals 5/3.Thus, the average photon energy in the observed spectrum (not accounting for redshift) is 5Γ/3 times larger compared to εph .Since the broadening is linear with Γ, we can model it as a post-processing parameter.We normalize the output broadened spectrum such that the average photon energy in the broadened spectrum equals εph .Thus, to obtain the observed spectrum, the energy grid of the broadened spectrum should be multiplied by 5Γ/3, which motivates the choice of the Doppler factor used in Section 5.In Figure 8, we show the effect of the broadening scheme on two spectra with different values of R.
The effect of the broadening is more prominent for the narrower spectrum.

D. ADDITIONAL CONFIGURATIONS
In this appendix, we present results similar to those shown in the main text for three additional configurations of the internal shock parameters.The first two (Figures 9 and 10) have different distributions for the higher-order initial parameters.For Figure 9, uniform distributions have been used, as given in Table 3.In Figure 10, the distributions are similar to the main text except that the average value of Γ s has been increased, as given in Table 4.The methodology is the same as when generating the histograms presented in Figure 4. Specifically, the cuts outlined in subsection 4.2 still apply.This leads to a removal of ∼ 30% of the generated initial parameter triplets in the case of Figure 9, mostly due to many parameter combinations resulting in τ θ > 50.For the third additional configuration (Figure 11), the distributions for the higher-order initial parameters is given in Table 1, i.e., the same as for the results presented in the main text.However, it uses the upstream temperature and photon-to-proton ratio obtained for the reverse shock instead of the forward shock in equations (B9) and (B10).This leads to softer spectra on average, with values for α reaching as low as −1.5 in the delayed acceleration scenario.In Figure 12, we show the obtained 2SBPL parameter distributions for the three additional configurations, similar to Figure 5.  3.This choice of initial parameters results in higher values of τi on average, leading to harder spectra.Many of these burst would likely be missed in observations due to the increase in adiabatic degrading and decrease in luminosity.4, but generated from the higher order initial parameter distribution given on the in Table 4.A higher initial value of Γs has a small effect on α but increases the observed peak energy.
Figure1.A schematic of the subphotospheric internal collision framework.Two fireballs are accelerated from the base in a conical outflow, starting at radius r0 and with an initial Lorentz factor Γ0. The slower of the two fireballs has a bulk Lorentz factor Γs and the faster fireball has a Lorentz factor Γ f = ψΓs.The collision occurs at an optical depth (as measured for the slower fireball) τs, with the requirement τs > 5 to avoid the dynamical changes that occur in the RMS at small optical depths.

Figure 3 .
Figure 3. Count spectrum (top) and νFν spectrum (bottom) obtained by fitting a CPL function to the (forwardfolded) photospheric spectrum in the delayed acceleration scenario from Figure 2. As the hardening of the input spectrum at low energies is outside of the GBM sensitivity range (indicated by the shaded area), the fit is adequate.In the bottom panel, the un-folded photospheric spectrum from Figure 2 has been overlaid for comparison.Best fit parameters are α = −0.92± 0.05 and E peak = 150 +13 −12 keV.

Figure 4 .
Figure4.Histograms obtained by fitting a CPL function to 150 synthetic photospheric spectra, using the method described in subsections 4.2 and 4.3.The mean and standard deviations for the input RMS and KRA parameters are given in Table2, which in turn are obtained from the higher order initial parameter distributions given in Table1.Top panels show α and bottom panels show E peak , for the early acceleration scenario (left) and the delayed acceleration scenario (right), as also indicated in the panels.The solid dark green lines show the distributions found for the CPL parameters in a time-resolved catalogue of bright Fermi GRBs during the first four years of observations(Yu et al. 2016).It is evident that unless collimation shocks delay the acceleration, the distribution obtained is harder than most observed bursts.All panels are made with χ = 1 and z = 1.

Figure 5 .
Figure5.Histograms obtained by fitting a 2SBPL function to 150 synthetic photospheric spectra in the delayed acceleration scenario, obtained using the method described in subsections 4.2 and 4.3.The insets show the subset with ∆AIC < 0. Left panel shows the low-energy indices α1 and α2 and the right panel shows the ratio of the peak energy to the break energy.The lines are smoothed kernels of each distribution.The excess of fits with α2 ∼ −2 is due to bursts where the hardening of the spectrum is below the detector energy window (see main text).These spectra are not better described by 2SBPL function compared to a CPL function as evident from the inset.When the Lorentz factor increases, more of the curvature is visible in the detector and the excess disappears, which can be seen in the top middle panel of Figure12.Both panels are made with χ = 1 and z = 1.

Figure 6 .
Figure 6. Figure showing θ v r (left axis, blue line) and εr (right axis, red line) as a function of optical depth for a simulation.The effective electron temperature θ vr initially grows as the Compton temperature in the shock grows.It tends to its steady state value θr, which is indicated by the horizontal dotted line.As expected, the average photon energy in the RMS, εr, remains at all times constant and equal to εd , which is indicated by the horizontal dashed line.The vertical dot-dashed line shows our estimate of τss calculated using equation (5).It slightly underestimates when equilibrium is reached, which is expected since equation (5) assumes a constant energy gain per scattering of 4θr, while the true energy gain is equal to θ v r , which is lower than θr initially.Adiabatic cooling is not included.Parameter values are τ θ = 5, θr = 0.1, R = 500, yr = 1.5.

FFigure 7 .
Figure7.Left: Evolution of the RMS spectrum in a simulation.Ten snapshot spectra are shown, taken at even intervals of τ ranging from the initial value τ = 50 (black) until the shock has crossed the upstream at τ = 25 (red).Initially, the shock contains a Wien thermal photon distribution with mean energy εd .As the shock evolves, the initial peak moves to higher energies and broadens, while new cold photons enter the shock from the upstream.At τ ≈ 33 (snapshot seven), steady state has been reached and the final spectra are overlapping.This corresponds well with the estimate shown in Figure6.Right: Similar to the left panel but for the downstream spectrum.This time the spectra range from τ = 50 all the way up until the photosphere at τ = 1 (from darker to brighter green).Note that the spectrum for τ = 50 is not visible since the downstream does not contain any photons at the start of the simulation.For this specific choice of parameters, the optical depth was large enough to erase all signatures of the shock formation in the final spectrum.Adiabatic cooling and the broadening effects explained in Appendix C are not included.Parameter values are the same as in Figure6.
τ, µ)D(µ)ϕ(τ ) dϵ ph dτ dµ = N ϵ ph dϵ ph f (τ, µ)dτ dµ = N ph f (τ, µ)dτ dµ = N ph ≡ N,(C22) Figure8.Two spectra before (solid lines) and after (dashed lines) the broadening effects have been added.The solid lines correspond to comoving spectra at the photosphere and the dashed lines to the shapes of the observed spectra, but normalized to have the same average photon energy as the photospheric spectra (see the text for details).The red spectrum has a smaller ratio between the peaks (R = 50) than the green (R = 1000).Both initial spectra have been normalized such that θu,K = 1.It is evident that the effect of the broadening scheme is most important for the narrower spectrum.Other parameter values are τ θ = 5 and yr = 1.5.

Figure 9 .
Figure 9. Similar to Figure 4, but generated from the higher order initial parameter distribution given in Table3.This choice of initial parameters results in higher values of τi on average, leading to harder spectra.Many of these burst would likely be missed in observations due to the increase in adiabatic degrading and decrease in luminosity.

Figure 10 .
Figure10.Similar to Figure4, but generated from the higher order initial parameter distribution given on the in Table4.A higher initial value of Γs has a small effect on α but increases the observed peak energy.

Table 1 .
Mean and standard deviations for the normal distributions used to generate values for the internal collision parameters.The logarithms used for Γs and τs are base 10.

Table 3 .
Lower and upper limits for the uniform distributions used to generate the histograms shown in Figure9

Table 4 .
Mean and standard deviations for the normal distributions used to generate the histograms shown in Figure10.The logarithms used for Γs and τs are base 10