{'fraction_non_alphanumeric': 0.03089606141395806, 'fraction_numerical': 0.012744152627299714, 'mean_word_length': 4.418852459016393, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
An efficient method for fitting radiation-mediated shocks to gamma-ray burst data: The Kompaneets RMS approximation

Shocks that occur below a gamma-ray burst (GRB) jet photosphere are mediated by radiation. Such radiation-mediated shocks (RMSs) could be responsible for shaping the prompt GRB emission. Although well studied theoretically, RMS models have not yet been fitted to data due to the computational cost of simulating RMSs from first principles. Here, we bridge the gap between theory and observations by developing an approximate method capable of accurately reproducing radiation spectra from mildly relativistic (in the shock frame) or slower RMSs, called the Kompaneets RMS approximation (KRA). The approximation is based on the similarities between thermal Comptonization of radiation and the bulk Comptonization that occurs inside an RMS. We validate the method by comparing simulated KRA radiation spectra to first-principle radiation-hydrodynamics simulations, finding excellent agreement both inside the RMS and in the RMS downstream. The KRA is then applied to a shock scenario inside a GRB jet, allowing for fast and efficient fitting to GRB data. We illustrate the capabilities of the developed method by performing a fit to a non-thermal spectrum in GRB 150314A. The fit allows us to uncover the physical properties of the RMS responsible for the prompt emission, such as the shock speed and the upstream plasma temperature.


INTRODUCTION
The launching, propagation and collimation of a highly supersonic jet unavoidably leads to immense shock formation inside the jet and its surroundings (see e.g., López-Cámara et al. 2013, 2014Gottlieb et al. 2019). Shocks that occur deep inside gamma-ray burst (GRB) jets are mediated by radiation. Such radiationmediated shocks (RMSs) fill the jet with hot, nonthermal radiation, which is advected toward the jet photosphere where it is released. The released spectra will range from strongly non-thermal to thermal, depending on whether the radiation has had time to thermalize via scatterings before reaching the photosphere or not.
GRBs are observed to have strong spectral evolution, both in terms of peak energy (Golenetskii et al. 1983) and shape (e.g., the width of the spectrum Wheaton et al. 1973). In around one quarter of GRB pulses, the narrowest, time resolved spectrum is consistent with a thermal spectrum, which strongly suggests that the Corresponding author: Filip Samuelsson filipsam@kth.se whole pulse is of a photospheric origin (Yu et al. 2019;Acuner et al. 2020;Dereli-Bégué et al. 2020;Li et al. 2020). It is plausible, therefore, that the wider, nonthermal spectra in such pulses have undergone subphotospheric dissipation (Rees & Mészáros 2005;Ryde et al. 2011). Even though RMSs are a natural cause of this dissipation, so far, RMS models have not been fitted to the data. 1 The main reason for this is that RMSs are expensive to simulate from first principles. RMSs have previously been considered in one spatial dimension (Levinson & Bromberg 2008;Nakar & Sari 2012;Beloborodov 2017;Lundman et al. 2018;Ito et al. 2018;Lundman & Beloborodov 2019;Levinson 2020;Ito et al. 2020;Lundman & Beloborodov 2021). The 1D simulations illustrate the main features of the nonthermal RMS radiation expected inside GRB jets: a broad power-law spectrum for up to mildly relativistic shocks (that have relative relativistic speed between the up-and downstreams of βγ few, where β is the speed in units speed of light c and γ is the Lorentz factor), while faster shocks have more complex spectral shapes due to Klein-Nishina effects, anisotropic radiation in the shock, and photon-photon pair production. Once advected into the downstream, the RMS spectrum gradually thermalizes through scatterings.
Currently, these 1D simulations are not fast enough to build a table model of simulated RMS spectra over the relevant parameter space, which is an efficient way to test models against data. However, model testing is of crucial importance to further develop our understanding of the prompt emission in GRBs. With this motivation, in this work we explore an alternative path to connecting RMS theory and GRB observations. In Section 2, we construct an approximate, but very fast, method called the Kompaneets RMS approximation (KRA). The approximation is based on the strong similarities between bulk Comptonization of radiation inside an RMS, and thermal Comptonization of radiation on hot electrons, the latter being described by the Kompaneets equation. The KRA is appropriate to use for mildly relativistic (and slower), optically thick RMSs. We validate the KRA by comparing simulated radiation spectra to those produced by the full radiation hydrodynamics simulations, finding excellent agreement. The KRA is then applied to a minimal model of a shock inside a GRB jet in Section 3, which generate synthetic photospheric spectra, accounting for both adiabatic cooling and thermalization of the photon distribution. The KRA simulations are about four orders of magnitude faster to run than the corresponding 1D simulations, allowing for table model construction. As an illustration of the model capabilities, we use the table model to perform a fit to the prompt emission of GRB 150314A in Section 4. We conclude by summarizing and discussing our results in Section 5.

THE KOMPANEETS RMS APPROXIMATION
In this section, we develop the KRA and compare the resulting spectra to full scale, special relativistic RMS simulations in planar geometry. The approximation is valid for RMSs where the photons inside the shock do not obtain energies exceeding the electron rest mass energy, as the transfer problem then becomes more complicated, including Klein-Nishina scattering effects, anisotropy, and γγ-pair production. This typically corresponds to shocks that are mildly relativistic 2 or slower, 2 We illustrate this point later with a shock simulation that has an upstream four-velocity of βγ = 3. inside plasma where the downstream radiation pressure dominates over the magnetic pressure. As is appropriate for RMSs inside GRB jets, the RMS is assumed to be photon rich (Bromberg et al. 2011), i.e., the photons inside the RMS are mainly supplied by advection of upstream photons, as opposed to photon production inside the RMS and in the immediate downstream. The RMS is also assumed to be in an optically-thick region, which is appropriate deep below the photosphere. The approximation will therefore not hold for shocks that dissipate most of their energy close to the photosphere; such shocks require full radiation hydrodynamics simulations. 3 We note that Blandford & Payne (1981) showed that the shape of a photon spectrum traversing a photon rich, nonrelativistic RMS can be obtained analytically. Although their analytical calculation is impressive, it ignores photon energy losses due to electron recoil. The photon spectrum, therefore, lacks a high-energy cutoff, which makes it accurate only for soft spectra where the bulk energy is not carried by the high-energy photons. Due to this limitation, their solution is not applicable here.

Bulk Comptonization inside the RMS
The following treatment assumes a nonrelativistic shock, but comparison to full RMS simulations show that the approximation is valid also for mildly relativistic shocks with γβ few (see also Section 2.4).
In the shock rest frame, the incoming speed of the upstream is greater than the outgoing speed of the downstream, leading to a speed gradient inside the shock. The photons that diffuse inside the RMS speed gradient directly tap the incoming kinetic energy by scattering on fast electrons. If the photon mean free path is λ, the velocity difference of the plasma over a scattering length is λ(dβ/dx) = dβ/dτ , where x is the spatial coordinate and dτ = dx/λ measures the optical depth along the x-coordinate. Doppler boosting the photon to the frame of the scatterer, performing a scattering and averaging over the scattering angles, one finds a relative energy gain of per scattering, e.g., a first order Fermi process. Here, the photon energy is given in units of electron rest mass, = hν/m e c 2 , where h is Planck's constant and ν is the photon frequency. Equation (1) is valid for a relative energy gain ∆ / < 1, where ∆ is the energy gain in a scattering for a photon with initial energy . Note that dβ/dτ is a local quantity that changes continuously across the RMS transition region, and vanishes in the far up-and downstreams, where the plasma velocity is constant (see Figure 1 for a schematic of the velocity profile across the shock). The dβ/dτ -profile is self-consistently determined by the radiation feedback onto the plasma: the photons gain precisely the available kinetic energy such that the Rankine-Hugoniot shock jump conditions are satisfied. Since plasma is advected through the RMS, so are the photons that scatter inside the plasma. However, photons also diffuse within the flow, and a fraction of the photons will stay inside the RMS much longer than the advection time across the RMS, accumulating more scatterings and therefore more energy. As is always the case when both the probability of escaping the shock and the relative energy gain per scattering, ∆ / , are energy independent, a power-law spectrum develops. The power-law extends up to energies where the energy gain per scattering is balanced by energy losses due to electron recoil, which occurs when ∆ / ≈ . This gives a maximum photon energy inside the shock of where the brackets on the right-hand side indicates a weighted average across the shock. The exact expression for ∆ / is difficult to determine from first principles. With¯ u and¯ d as the average photon energies in the up-and downstreams, respectively, and β u γ u the four-velocity of the upstream evaluated in the shock rest frame, we empirically find in Appendix A that with ξ = 55 is a good approximation across the relevant shock parameter space. Although Equation (3) contains the relativistic four-velocity, it is only valid while ∆ / ≈ max 1.

Modelling an RMS as thermal Comptonization
The energy gain process described in Section 2.1 looks strikingly similar to thermal Comptonization on hot electrons (see e.g., Rybicki & Lightman 1979). Consider a hot cloud of nonrelativistic electrons at a constant temperature θ = kT /m e c 2 1 where k is the Boltzmann constant and T is the temperature, with injection of low energy photons ( θ) into the cloud, and an escape probability that is energy independent. The low-energy photons will gain a relative energy per scattering of ∆ / ≈ 4θ, and the energy gain continues until balanced by recoil losses at max ≈ ∆ / ≈ 4θ. Such Comptonization is described by the Kompaneets equation, with a source term s for the photon injection and escape from the cloud where t sc = λ/c is the Thompson scattering time and n is the photon occupation number. Stimulated scattering (∝ n 2 ) has been omitted in Equation (4), as this effect is insignificant as long as the occupation number is small, n 1, which is true for the non-thermal emission considered here.
Motivated by the similarities between the two systems, our aim is to construct an approximate RMS model based on the Kompaneets equation. The plasma is split into three discrete zones: the upstream zone, the RMS zone and the downstream zone. The time evolution of the radiation spectrum inside each zone is computed using the Kompaneets equation. Each zone has an effective electron temperature, and the zones are connected via source terms. A schematic of the KRA is shown in Figure 1.
The upstream zone feeds thermal radiation into the RMS zone, which passes radiation on to the downstream zone. Dissipation occurs only inside the RMS zone. This is achieved by prescribing an effective electron temperature θ r ≈ 1 4 ∆ / , found from Equation (3) as The subscript r here and henceforth denote quantities in the RMS zone, and the subscripts u and d will be used to denote quantities in the upstream and downstream zones, respectively. Equation (5) assures that the maximum photon energy and the energy gain per scattering in the RMS zone mimic those of a real RMS. By matching how long photons stay in the shock such that the average downstream energies in the two systems become equal, the evolution of the photon distribution in the KRA will closely match that of a real shock. This is achieved by using appropriate source terms (see Section 2.3).
The up-and downstream zones do not dissipate energy. Therefore, the temperatures inside these zones equal the radiation Compton temperature θ C , defined as the electron temperature with which there is no net energy transfer between the photon and electron populations. It is given by where the integrals are taken over all photon energies. The upstream temperature of the radiation 4 , θ u , is a free parameter of the model, but the downstream Compton temperature is not a free parameter, as it is determined by the upstream temperature and the amount of dissipation in the shock. The non-thermal radiation that streams from the RMS zone is accumulated in the downstream zone, where it gradually thermalizes 5 via scatterings. The downstream zone contains all photons that passed through the RMS zone, and the degree of thermalization of the radiation inside the downstream zone increases with time.

The KRA source terms
The three zones in the KRA are coupled by source terms. Denoting the source of photons that stream into the RMS by s in and the source that streams out of the RMS by s out , one gets The probability for a photon to escape the RMS into the downstream is independent of the photon energy. Thus, s out = kn r , where k is a constant and n r is the occupation number inside the RMS zone. In this scenario, one can show (e.g., Rybicki & Lightman 1979) that the steady state solution to the Kompaneets equation inside the RMS zone is a power-law distribution, n ∝ −α , with α = 3/2 ± (9/4 + k/θ r ) 1/2 . In analogy with Rybicki & Lightman (1979), we identify the RMS zone y-parameter as y r = 4θ r /k. Therefore, The RMS zone y-parameter determines how much time photons spend inside the shock and it is, therefore, a measure of the average photon energy gain inside the RMS. As such, y r sets the hardness of the non-thermal spectrum that is injected into the downstream. A value of y r = 1 corresponds to a flat νF ν -spectrum, with larger values of y r giving harder slopes. The value of y r in the KRA is chosen such that the average downstream photon energy obtained equals that of a real RMS. The full conversion between the parameters that specify the RMS and the corresponding KRA parameters is shown in Appendix A.
Requiring that the photon number inside the RMS zone is conserved one finds from Equation (8) where the integrals are again taken over all photon energies.

Estimating the KRA upper speed limit
The energy gain inside the RMS qualitatively changes when the relative energy gain per scattering, ∆ / , becomes comparable to unity. This is because the upper photon energy inside the shock is max ≈ ∆ / , and the radiative transfer is different for such high-energy photons. In particular, Klein-Nishina effects modify the scattering cross-section and γγ-pair production is triggered. Furthermore, shocks with such a high energy gain per scattering will have a narrow width (comparable to a few photon mean free paths), which makes the radiation anisotropic. The KRA is therefore limited to modelling shocks with ∆ / 1. This corresponds to (see Equation (3)) with ξ ≈ 55.
The ratio of average downstream to upstream photon energies can vary significantly but enters Equation (12) only as a logarithmic factor. For a typical energy ratio of¯ d /¯ u = 10 2 , one gets an upper velocity limit of β u γ u ≈ 3.5. Thus, the KRA is expected to be applicable to shocks with β u γ u 3, with the exact value only marginally dependent on the shock parameters.

Quasi-thermal RMS spectra
The radiation in the downstream of an RMS becomes quasi-thermal if the energy dissipation per photon is either very low or very high. In the former case, the upstream photon distribution is largely unaltered, while in the latter case, the photons gain so much energy that they pile up in a thermal Wien distribution around max (i.e., saturated Comptonization). In both cases, the radiation relaxes to a near thermal distribution after a few scatterings in the downstream, at which point the information from the shock is all but lost. When fitting to data, such shocks are almost indistinguishable from each other and from outflows where no dissipation occurred. As such, they are less interesting from an observational perspective.
Shocks with small photon-to-proton ratios, n γ /n p , where n γ and n p are the photon and proton number densities, respectively, tend to have more thermal-like spectra inside the RMS. This is because the average downstream photon energy¯ d is inversely proportional to the photon-to-proton ratio (i.e., more photons sharing the same shock kinetic energy), while the maximum photon energy max is proportional to the logarithm of¯ d (see Equation (3)). Thus, as the photon number shrinks,¯ d increases faster than max , until the spectrum appears quasi-thermal with¯ d ∼ max .
For¯ d ¯ u , the average downstream photon energy is¯ where m p is the proton mass. Equating¯ d to max ≈ ∆ / and solving for n γ /n p , one gets Consider the limit β u 1, which gives (γ u − 1)/(β u γ u ) 2 ≈ 1/2. With ξ = 55, and a typical energy ratio of¯ d /¯ u = 10 2 , one finds a critical photon-to-proton ratio of n γ /n p ≈ 1.1 × 10 4 . Shocks that have photonto-proton ratios close to this value will result in quasithermal radiation spectra. Below, we illustrate this fact with a simulation that has n γ /n p ≈ 4 × 10 4 .

Comparing the Kompaneets RMS approximation to full RMS simulations
In this subsection, we will compare the spectrum inside the RMS and downstream regions as computed by the two codes radshock and Komrad, the latter implementing the KRA. The radshock code is a special relativistic, Lagrangian radiation hydrodynamics code (Lundman et al. 2018). The radiation field is computed using the Monte Carlo method, which self-consistently connects to the hydrodynamics via energy and momentum source terms. The RMS is set up by smashing plasma into a wall boundary condition, and allowing the code to relax into an RMS that propagates steadily away from the wall. For the case of a thermal upstream radiation spectrum, the RMS solution is fully specified by three parameters. These can be taken to be the temperature of the upstream radiation, θ u , the speed of the upstream plasma relative to the shock, β u , and the photon-to-proton ratio, n γ /n p , inside the upstream (Lundman et al. 2018).
Komrad implements the KRA described in the previous subsections, evolving the radiation in the RMS zone and the downstream zone 6 using Kompaneets solvers (e.g., Chang & Cooper 1970). We choose the following three parameters to describe the RMS in Komrad: the temperature of the upstream photons, θ u,K (where the subscript K indicates Komrad), the effective electron temperature inside the RMS zone, θ r , and the Compton y-parameter of the RMS zone, y r . As mentioned above, the conversion from the KRA parameters to the corresponding RMS parameters is given in Appendix A. A non-trivial point is that the two codes will have somewhat different upstream temperatures. This is because plasma compression inside the RMS will increase the internal energy density and shift the upstream spectral peak, and no analogue to this compression exists for the KRA.
A simulation is fully specified by the three shock parameters and the total simulation time t/t sc . The simulation time affects the degree to which the downstream has been thermalized. To highlight the similarities between the downstream spectra, we do not include the radiation produced during the initial RMS formation.  This is because the formation of the shock is different between the simulations. Therefore, the simulation time starts when the RMS is already in steady state. The RMS transition region in radshock is continuous and it is not obvious a priori what part should be compared to the RMS zone in Komrad. For the comparison, we chose the radiation in radshock that is located at the point where the shock has just finished dissipating all incoming energy, as this represents the spectrum that is injected into the downstream. This location is determined as the point where the average photon energy has reached its downstream value. The plasma that has passed through this location after the start of the simulation time belongs to the downstream. At the end of the simulation, the radiation inside the downstream is collected and its spectrum computed. Figures 2 and 3 show the results of six different shock simulations, labeled Run A-Run F. The Komrad parameters for the six runs are shown in Table 1 and the radshock parameters are found in Table 2. The simulation parameters were chosen to test the KRA in different regions of the shock parameter space, resulting in differently shaped radiation spectra inside the shock. The Komrad parameters for the six runs are calculated from the corresponding radshock parameters using the method described in Appendix A. The only free parameter in the conversion is ξ from Equation (3). All Komrad runs are made with ξ = 55, as we empirically found this value gave good agreement across the parameter space.
In Figure 2, Run A-Run E are shown in five different panels. The spectra produced by the two codes are remarkably similar, highlighting the close analogy between bulk and thermal Comptonization. We conclude that the KRA can accurately capture the RMS radiation physics.
In Figure 3, we show the spectra for Run E, which is a mildly relativistic shock with upstream speed β u γ u = 3 in the shock rest frame. The KRA neglects relativistic effects such as Klein-Nishina suppression and pair production. Furthermore, as shown in e.g., Ito et al. (2018), anisotropy starts to become important when the shock becomes relativistic. However, Komrad can still capture the behavior of mildly relativistic shocks, as long as the photon energies inside the RMS do not exceed the electron rest mass, i.e., as long as max 1 as discussed in Section 2.4.
In this run, the relative energy gain per scattering is close to unity and the shock is only a few Thomson optical depths wide in the radshock simulation. Hence, photons have a high probability of diffusing in and out of the different regions and there are no sharp "zone boundaries" in radshock. Therefore, the comparison to a discrete RMS zone in Komrad is less accurate. Furthermore, Klein-Nishina suppression starts to become important. This can be seen in the high-energy tail of the photon distribution in the downstream. The highenergy photons in radshock have cooled less than those in Komrad, due to their lower scattering cross section. However, this effect will likely be unimportant when the KRA is fitted to actual data, as radiation with 1 have time to downscatter to lower energies before reaching the photosphere, even with Klein-Nishina suppression. Furthermore, the high-energy part of the spectrum is often given little weight in the fitting process due to the lower photon counts at high energies in the GBM energy channels (Yu et al. 2019). Overall the approximation is surprisingly accurate even in this case when β u γ u = 3, especially in the downstream zone, which contains the radiation that will later be observed. This indicates that any anisotropy of the radiation field within the shock does not have a major impact on the shape of the spectrum in this case. We conclude that the limit of the KRA is when max ≈ ∆ / ≈ 4θ r starts to approach unity.

APPLYING THE KOMPANEETS RMS APPROXIMATION TO A GRB JET
RMSs come with a variety of dynamical behaviors. Explosions, such as supernovae, generate an outward going shockwave. The shockwave propagates through the star until it either breaks out of the stellar surface (i.e., the photosphere), or dissolves as the downstream pressure becomes too small, due to the limited explosion energy budget. A different dynamical behavior is seen in recollimation shocks, which arise as the jet propagates in a confining medium. Such shocks can be approximately stationary with respect to the star and might therefore never break out. However, the radiation that is advected through the recollimation shock is energized, and the emission released at the photosphere can be non-thermal. Yet another behavior is seen in shocks that arise due to internal collisions of plasma inside the jet. When two plasma blobs collide, the plasma in between the blobs is compressed, increasing the pressure adiabatically until the pressure profile is steep enough to launch two shocks. The shocks propagate in opposite directions into the two colliding blobs while sharing a causally connected downstream region. Such shocks cease when they have dissipated most of the available kinetic energy in the two blobs. The time it takes the shocks to cross the blobs is roughly a dynamical time (as they cross the causally connected jet ejecta). Dynamical effects on the shock structure are important if the shock reaches the jet photosphere. For instance, part of the RMS can transform into a pair of collisionless shocks at the point of breakout when the photons mediating the shock start to leak out toward infinity (Lundman & Beloborodov 2021). The KRA is not able to handle such dynamical effects. On the other hand, the KRA is well suited for simulating plasma that is shocked while still being optically thick, and then ad- Figure 4. Schematic showing four stages in the evolution of the minimal shock model from left to right. Green, red, and purple indicate the upstream, RMS, and downstream zones, respectively. At some optical depth τi (small r), the shock is initiated. Photons start to diffuse through the shock region and gain energy. The photons that pass through the RMS are collected in the downstream region, which gain more and more photons (stages 1 and 2). When the RMS has crossed the upstream it dissolves, leaving only the downstream (stage 3). The photons in the downstream continues to scatter until they are released at the photosphere (stage 4).
vected toward the jet photosphere where the emission is later released.

A minimal sub-photospheric shock model
One shock scenario that can be modeled by the KRA is the collision inside the jet of two blobs of similar mass and density, but different speeds. This is a minimal scenario with few free parameters. More complex models with additional parameters can be considered in the future, if the current model fails to fit prompt GRB data. In Appendix B we compute the speed of the two shocks, the energy dissipated, and the radius at which the shocks have finished dissipating most of the available energy. For blobs of similar properties (i.e., similar mass and density), we show that also the properties of the shocks are similar. In that case, only one of the two shocks have to be explicitly simulated, and the number of model parameters are kept at a minimum.
The KRA is valid for shocks where ∆ / 1, as discussed in Section 2.4. In Appendix B we translate this limit to jet quantities in the context of internal shocks. It corresponds to a Lorentz factor ratio of Γ 2 /Γ 1 30, where Γ 2 and Γ 1 are the lab frame Lorentz factors of the fast and slow blobs, respectively. As an example, Γ 1 = 50 and Γ 2 = 1000 produces shocks that the KRA can accurately model.
A schematic illustration of the minimal shock model at four different stages of its evolution is shown in Figure 4. In the first stage, the two blobs have recently collided. The RMS has started to propagate into the upstream and a few photons have had time to diffuse into the downstream. In the second stage, the shock has almost crossed the upstream. The shock has finished crossing the upstream and dissolved in the third stage, with all photons accumulated in the downstream, and the fourth stage shows the radiation being released at the photosphere. Each of the three zones account for thermalization of the photon spectrum via scatterings and adiabatic cooling.
The time over which shocks dissipate their energy is not a free parameter; it should be found self-consistently from, e.g., hydrodynamical simulations. However, the shock crossing time is always comparable to the dynamical time scale of the jet, which corresponds to a doubling of the jet radius. Therefore, we let the KRA dissipate energy over a doubling of the jet radius in our shock model. This also corresponds to a halving of the optical depth as we consider the case of a conical outflow. After the shock stops dissipating, the downstream plasma containing the shocked radiation is advected toward the photosphere, while it gradually thermalizes through scatterings and cools adiabatically. The simulation ends when the shocked radiation reaches the photosphere.
As mentioned above, we omit photon production by the plasma (i.e., the shock is photon rich). This is a valid assumption, as the advected flux of upstream photons already existing inside the GRB jet is much larger than the number of photons produced by bremsstrahlung or double Compton scattering (e.g., Bromberg et al. 2011;Lundman & Beloborodov 2019). Photon production will occur in the downstream but the time scale for such production is long; the photon spectrum thermalizes into a Wien spectrum via scatterings long before photon production acts to modify the Wien spectrum into a Planck spectrum. As shown by Levinson (2012), photon production has time to modify the spectrum if the shock occurred at optical depths of ∼ 10 5 . In that case, the radiation will have lost essentially all its energy to adiabatic expansion before reaching the photosphere, and is, therefore, of little interest.

KRA implementation in spherical geometry
A conical jet appears locally as spherically symmetric. The Kompaneets equation inside a steady state, spherical relativistic outflow (with outflow bulk velocity β → 1) is given by Equation (3) wherer = r/R ph is a normalized radius and R ph is the radius of the photosphere. The normalized radius equalsr = 1/τ , where τ = σ T n e r/Γ is the optical depth of the jet (not to be confused with the optical depth of the RMS) with σ T being the Thomson cross section and n e the electron number density. The comoving time coordinate in the Kompaneets equation has here been re-written into a lab frame radial coordinate (using t = r/Γc). The last term in the curly brackets of Equation (15) accounts for adiabatic cooling of the spectrum. 7 In the optically thick regime, this causes the average energy of the photon distribution to decrease asr −2/3 , while the shape of the spectrum is preserved. When the photons start to decouple close to the photosphere (τ 1), the evolution changes and the idealized cooling ofr −2/3 is no longer valid (Pe'er 2008;Beloborodov 2011). To account for this, we numerically stop the cooling at an optical depth of τ = 3. The total adiabatic cooling of the photon distribution is then similar to that of a real spectrum where the proper radiation transfer is taken into account (see, Beloborodov 2011). Scattering is incorporated until τ = 1.

Lab frame transformation of the simulated radiation spectrum
The simulation outputs the comoving radiation spectrum at the jet photosphere. The simplest approximate transformation to the lab frame involves multiplying all photon energies by a factor Γ (the Doppler boost for a typical photon), where Γ is the Lorentz factor of the downstream zone. Since Γ does not explicitly enter the simulation, it is effectively a post-processing parameter.
In reality, the radiation spectrum broadens somewhat as it decouples from the plasma at the jet photosphere (Pe'er 2008;Beloborodov 2010;Lundman et al. 2013). This is because individual photons decouple at different angles to the line-of-sight, which affects their Doppler boosts, and also at different radii, which affects their energy losses due to adiabatic expansion. These effects are important to take into account when performing spectral fits to data, specifically to narrow bursts (Ryde et al. 2017). This spectral broadening can be approximately computed in a post-processing step, under the assumption that the jet Lorentz factor Γ is constant at the photosphere. The post-processing calculation is fairly long, and will be described in full detail elsewhere.
The Kompaneets equation without the induced scattering term is linear in the photon occupation number n. Therefore, the total photon number of the simulation is also a free parameter, which effectively makes the normalization of the GRB luminosity a post-processing parameter.

New parameters based on parameter degeneracy
In the case of planar geometry described in Section 2, the KRA has three parameters: θ u , θ r , and y r . In the case of a jet, the optical depth of the jet where the shock is initiated, τ i , is an additional parameter. Together, they determine the shape of the photospheric spectrum in the rest frame of the outflow. Additionally, as described in the subsection above, the bulk Lorentz factor Γ and the luminosity L γ are two post-processing parameters that shifts the observed spectrum in energy and flux.
However, there exists a degeneracy in the current parameter set. As long as the product τ i θ r , the ratio θ r /θ u , and y r remain unchanged, the spectral shape in the rest frame of the outflow is identical. With the translational freedoms given by Γ and L γ , this becomes degenerate. Therefore, a more suitable set of parameters are τ θ ≡ τ i θ r , R ≡ θ r /θ u , and y r . This brings the number of simulation parameters down from four to three, which simplifies the process of table model building.
The degeneracy can be understood as follows. Imagine τ i is increased by some factor f but θ r and θ u are decreased by the same amount. The evolution of the photon distribution with optical depth is slower, as the energy transfer per scattering is proportional to the electron temperature (when the photon gains energy in the scattering) or the photon energy (when the photon loses energy), both of which have decreased by a factor f . However, the number of scatterings is f times larger, so the relative energy transfer is the same, i.e., the spectrum evolves similarly but is a factor f lower in energy. The net effect is that the evolution of the whole system is equivalent. The shape of the photospheric spectrum is identical, but shifted down in energy by a factor f 5/3 , where an additional factor f 2/3 comes from the increased adiabatic cooling. The degeneracy in number of scatter-ings and energy gain per scattering is not unique to our model. Indeed, it is inherent to all jetted RMS models.
To see how each parameter influences the shape of the released photospheric spectrum, in Figure 5 we vary τ θ (top panel), R (middle panel), and y r (bottom panel), while keeping the other parameters constant. The value of the parameter being varied increases from black to red. As can be seen in the figure, the combined parameter τ θ determines the amount of thermalization after the shock has finished dissipating its energy. A higher τ θ implies a higher number of scatterings and/or higher energy transfer per scattering, leading to a faster thermalization. For large τ θ the downstream spectrum relaxes to a Wien spectrum, in which case the original shock parameters cannot be retrieved. The ratio R ≡ θ r /θ u determines the separation between the lower and upper cutoff in the spectrum. A large R leads to a long power-law segment in the downstream. The slope of the power-law depends on the Compton y-parameter y r , which is a measure of how much energy is dissipated in the shock. Higher values of y r lead to harder spectra, with y r = 1 corresponding to a flat νF ν -spectrum. The spectral broadening discussed in Section 3.3 has been omitted in Figure 5, so that the effect of each parameter on the final spectrum is more clearly seen.

FITTING GRB DATA WITH THE KOMPANEETS RMS APPROXIMATION
As a proof of concept of fitting an RMS model to prompt GRB emission data, we present an analysis of a time-resolved spectrum in GRB 150314A. This luminous burst was observed by the Fermi Gamma-ray Space Telescope and its Gamma-ray burst monitor (GBM), which covers the energy range of 8 keV-40 MeV. GRB 150314A is an example of a GRB pulse in which the spectrum becomes very narrow during a portion of its duration (Yu et al. 2019). The low-energy photonindex, α, of the Band-function (Band et al. 1993) reaches a very large value, α max = −0.27. Such a large α-value strongly suggest a photospheric origin of the emission during the analysed time-bin . It is further natural to assume that the same emission mechanism operates throughout a coherent pulse structure such as the one in GRB 150314A (Yu et al. 2019). Therefore, the whole pulse can be argued to be photospheric, even though most of the other time bins have non-thermal spectra (α ∼ −1). As described in the introduction, these non-thermal spectra must then have been formed by subphotospheric dissipation (Rees & Mészáros 2005;Ryde et al. 2011) with RMSs as the most probable source of dissipation (e.g., Levinson & Bromberg 2008;Lundman & Beloborodov 2019). The value of the varying parameter increases from black to red, being evenly log-spaced from 1.5 to 50 for τ θ, 10 to 10 3 for R, and 0.5 to 3 for yr. Increasing τ θ increases the thermalization. The ratio R determines the separation between the lower and upper cutoff in the spectrum and yr determines the slope of the power-law segment. All spectra have been normalized to unity in photon number (Nν = 1) at = 1.
In order to perform fast and efficient fits with the RMS model, we generate synthetic photospheric spectra over a large parameter space using Komrad within the minimal shock model as explained above. With the photospheric spectra, we construct a table model in the Multi-Mission Maximum Likelihood Framework (3ML; Vianello et al. 2015). Here, we include the broadening effect described in Section 3.3 by post processing of the spectra. Our initial table model consists of 125 spectra. Figure 6 shows a fit to one of the non-thermal spectra in GRB 150314A with our RMS model. The data is from a narrow time-bin at around 4.6 seconds after the trigger. The fit indicates that there are two spectral breaks present, one at around 30 keV and the other at around 400 keV. The high-energy break corresponds to the energy the high-energy photons, those that reached the maximum energy max ∼ 4θ r in the RMS, have downscattered to before decoupling. Conversely, the low-energy break depends on how much the low-energy photons, those that entered the downstream with energy ∼¯ u , are heated by scatterings before they reach the photosphere (see also Section 5.1). A corresponding fit with the Band function yields α = −0.73 ± 0.06 and a high-energy index β = −2.47 ± 0.25, which are typical values for non-thermal spectra in GRBs (e.g., Yu et al. 2016). Both the Band function and the RMS model have AIC = 1610, and can therefore equally well describe the data. The model parameters of the best fit of the RMS model are τ θ = 11.3 ± 2.9, R = 290 ± 50, and y r = 1.72 ± 0.14. The initial separation between θ r and θ u was thus relatively large and the thermalization moderate, which allow for the broad, non-thermal shape of the spectrum. The slope of the power-law segment at around 200 keV reveals that quite a lot of energy has been dissipated in the shock, causing a large value of y r . Assuming a Lorentz factor of Γ = 300, the parameters can be decoupled. 8 This gives θ r = 0.055 and τ i = 206. Translating this into the physical RMS parameters yields the upstream 4-velocity as β u γ u = 1.89, an upstream temperature of θ u = 8.81 × 10 −5 , and a photon to proton ratio of n γ /n p = 2.01 × 10 5 . This fit thus illustrates that our model can be used to study the flow properties and the shock physics in observed GRBs.

DISCUSSION AND CONCLUSION
Dissipation in the optically thick regions of a GRB jet has the potential to generate a wide variety of released photospheric spectral shapes, and it is, therefore, . Time-resolved spectrum in GRB 150314A from a narrow time-bin at around 4.6 seconds after the GBM trigger. Upper panel: νFν -spectrum of the best fit RMS model. Two breaks are present at around 30 and 400 keV. The best fit model is depicted by the black line and the grey region is its statistical uncertainty. The data points are derived from the counts fit, and correspond to three of the triggered GBM detectors. Note that the non-linearity of the GBM response matrix means that the data points will not be accurate in a νFν -spectrum, they are here shown only for visual purpose. Lower panel: The best fit model to the observed count data, including the residuals, which show random variation, indicating a good fit. a promising candidate for the prompt emission. Although RMSs are a natural dissipation mechanism, so far, no such model has been fitted to data. In this paper, we have for the first time performed a fit to a timeresolved spectrum of the prompt emission in a GRB using an RMS model. This allowed us to determine the physical properties of the initial shock, such as its speed and the upstream photon temperature.
The main reason for the previous lack of fitted prompt spectra within an RMS framework is that RMSs are computationally expensive to simulate from first principles. To overcome this obstacle, we developed an approximate model (KRA; Kompaneets RMS approximation, see Figure 1 for a schematic) based on the similarities between the bulk Comptonization of photons crossing an RMS and thermal Comptonization of photons on hot electrons. By comparing the simulated spectra from Komrad, a code employing the KRA, to those generated by a special relativistic radiation hydrodynamics code, we verified that the KRA can indeed accurately reproduce the RMS and downstream spectra from the full simulations in a wide parameter range (see Figures  2 and 3).
We connected the KRA to GRB prompt observations by creating a minimal shock model considering a single RMS occurring well below the photosphere. The downstream of the shock is allowed to thermalize and cool adiabatically as it advects to the photosphere, where its radiation is released (see Figure 4 for a schematic). The model has only three free parameters determining the shape of the released spectrum: the combined parameter τ θ that determines the amount of thermalization, R that determines the extent of the power-law segment, and y r that determines the hardness of the power-law (see Figure 5). Additionally, there are two, post-processing parameters for the normalization and the frequency-shift. We generated 125 spectra using the model and, after accounting for broadening of the observed spectrum due to high-latitude effects and a radially varying photosphere, performed a fit to a broad spectrum in a narrow time bin of GRB 150314A as a proof of concept (see Figure  6).

Qualitative spectral features
Within the minimal shock model developed in this paper, there are some clear observational predictions. The spectra will consist of smooth low-and high-energy cutoffs, with a power-law segment in between. The lowenergy cutoff is very smooth due to the broadening effects discussed in Section 3.3, while there is an exponential cutoff at the highest energies. Typically, a singlebreak function is used to fit the spectral GRB data, e.g., the Band function, with its peak at E p . Depending on the hardness of the power-law segment, E p can either correspond to the low-energy cutoff (when y r < 1) or the high-energy cutoff (when y r > 1). Breaks both above and below the peak energy have been detected. Additional breaks at low energies (< 10 keV) were reported in, e.g., Strohmayer et al. (1998), while additional highenergy breaks are discussed in, e.g., Barat et al. (1998). Within our model, bursts that have an additional lowenergy break should be well fitted with an exponential cutoff above E p or very soft values of the high-energy power-law index β in the Band-function. Conversely, bursts that have reported high-energy breaks above E p should produce hard low-energy slopes, as long as the smooth curvature is within the detector energy range. However, we note that many of our generated spectra will appear as a single smooth curvature over a large range of energies (see for instance the best fit in Figure  6), due to the effects of thermalization on the downstream spectrum, as well as the broadening effects of high-latitude emission and radially varying emission.
As it propagates toward the jet photosphere, the downstream spectrum will tend toward a Wien spectrum at the Compton temperature. The more thermalized the downstream spectrum becomes, the more difficult it will be to retrieve the original shock parameters. Once the spectrum has relaxed into a Wien spectrum, the shock information is lost. This is an inherit degeneracy in photospheric models that is important to keep in mind when drawing conclusions about the physics from the parameter estimation.

Optical emission
Although the curvature of the spectrum is very smooth and although it may be outside the observable energy range of the prompt detectors, the Rayleigh-Jeans limit always exists in our spectra at low energies. Therefore, our minimal model with a single shock cannot account for low-energy observations such as optical during the prompt phase. Thus, we must conclude that any early optical emission is part of the afterglow. Early optical observations are rare, and very few are within ∼ 100 s after trigger (see Oganesyan et al. (2021) for a recent example). Optical detections are commonly reported as prompt as long as they are observed within the T 90 of the GRB 9 (Yost et al. 2007;Klotz et al. 2009). This definition disregards whether the GRB had a quiescent period within the active phase or not. Given that most optical observations occur quite late, they could be the onset of the afterglow. This highlights the need for early optical observations in GRBs, which have the power to discriminate between the current models of the prompt emission (see also Oganesyan et al. 2019). Early optical observations also have the potential to discern if GRBs are significant contributors to the observed ultrahigh-energy cosmic ray flux, as discussed in Samuelsson et al. (2019Samuelsson et al. ( , 2020. Recollimation shocks below the photosphere and their connection to the prompt emission in GRBs have been investigated by several authors (e.g., Gottlieb et al. 2019). Although not discussed in this paper, we expect the KRA to be able to model recollimation shocks as well. Such shocks have different dynamics, but bulk Comptonization is still responsible for the energy dissipation, leading to the same spectral features. That is, a power-law segment with a cutoff at high energies and a Rayleigh-Jeans slope a low energies. Indeed, oblique shocks, such as recollimation shocks, can be transformed into parallel shocks with a suitable Lorentz transformation (Henriksen & Westbury 1988). Therefore, a recollimation shock could plausibly be responsible for the subphotospheric dissipation in a GRB whose spectra can be well fitted with the KRA. The minimal shock model considered here consists of a single RMS, dissipating energy over a dynamical time. It is easy to imagine a more complex jet structure with multiple shocks and turbulence. However, although the dynamics below the photosphere are complicated, it is not inconceivable that the shape of a time-resolved spectrum is dominated by a single, strong dissipation event. The good fit to the emission in GRB 150314A shows that the current minimal shock model can plausibly explain the data. Additional model complexity should be considered only if the current model is found inadequate to explain the observations. Further investigation will tell if the current minimal model is sufficient when applied to a larger sample of GRBs.
From a Komrad simulation, we know θ u,K , θ r , and¯ d . Given the equation above, the system can be solved. Numerically, one can start by guessing u d . Then u u and θ u can be found from equations (A6) and (A7). With equations (A4) and (A5), the only unknown left is n γ /n p , which can be solved from Equation (A1). If the original guess of u d was correct, Equation (A2) should be satisfied.
If one wishes to instead go from the RMS parameters to the Komrad parameters, one can find w d and u d numerically through equations (A1) and (A2), using equations (A4) and (A5). Then θ u,K and θ r are found from equations (A6) and (A7). The parameter y r can be found iteratively by requiring that the downstream energy¯ d should be equal in both models. In practice, a qualitative first guess of y r can be made from a plot of the RMS spectrum from radshock: by comparing the power-law slope to the spectra in the lower panel of Figure 5, the value of y r can be estimated.

B. WHEN CAN THE KRA MODEL INTERNAL SHOCKS?
Consider a part of the jet that consists mainly of two masses: a slower and a faster mass with lab frame Lorentz factors Γ 1 1 and Γ 2 Γ 1 , respectively. The masses are assumed to be initially separated by a lab frame distance δr δl 1 , δl 2 , where δl 1 and δl 2 are the corresponding initial lab frame widths of the slower and faster mass, respectively. The faster mass catches up to the slower mass at radius R i ≈ 2Γ 2 1 δr. By then, the plasma between the masses has been highly compressed, increasing its pressure adiabatically until a forward and a reverse shock forms. The forward and reverse shocks propagate into the slower and faster masses, respectively. The speed of the shocked region (i.e., the shared downstream, which is bounded by the forward and reverse shocks) is found by balancing the momentum flux in the rest frame of the shocked region, β 2 Γ 2 h + p, from both sides. Here, h ≡ ρc 2 + e + p is the specific enthalpy, where ρ is the mass density and p is the pressure, all of which are measured in the respective rest frames of the unshocked masses. If we suppose that the initial pressure inside the two masses is small (such that h ≈ ρc 2 ), then we can solve for the lab frame Lorentz factor Γ of the shocked material as where ρ 1 and ρ 2 are the proper densities of the respective mass, before being shocked. In a wide range of density ratios, (Γ 2 /Γ 1 ) 2 ρ 2 /ρ 1 (Γ 1 /Γ 2 ) 2 holds and the above expression can be simplified to The condition (Γ 2 /Γ 1 ) 2 ρ 2 /ρ 1 (Γ 1 /Γ 2 ) 2 also ensures that Γ 2 Γ Γ 1 . The radii where the two shocks have crossed their respective masses are then R 2 ≈ 2Γ 2 δl 2 and R 1 ≈ 2Γ 2 1 δl 1 respectively. The masses are related to their widths and densities by δm ≈ 4πr 2 Γρδl, and so the ratio of the radii where the reverse and forward shocks have crossed the respective masses can be written as The relative Lorentz factors between the upstream (moving with Lorentz factor Γ 1 or Γ 2 for the forward and reverse shocks, respectively) and the downstream (Lorentz factor Γ) gives a measure of how relativistic the two shocks are. They can be computed using Equation (B9), and their ratio isΓ The energy dissipated for each mass (in the rest frame of the shocked plasma) is E ≈ (Γ − 1)δmc 2 , whereΓ is the relative Lorentz factor between the up-and downstream, so that Based on the analysis above, we see that the collision of two masses with similar properties, δm 1 ∼ δm 2 and ρ 1 ∼ ρ 2 results in shocks of similar strengths,Γ 1 ∼Γ 2 . The shocks also dissipate roughly the same amount of energy, E 1 ∼ E 2 , and finish dissipating roughly at the same time, R 1 ∼ R 2 . The heated radiation is located in the shared downstream between the two shocks. Since the shocks have similar strengths, and the heated radiation from both shocks sits inside plasma that propagates with the same Lorentz factor (here called simply Γ), modeling of only one shock is necessary.
The KRA can accurately model shocks as long as the relative energy gain per scattering is small less than unity, ∆ / 1. In Appendix A we found that ∆ / ≈ 0.018(γ u β u ) 2 log(¯ d /¯ u ), which means the approximation is valid up to β u γ u ∼ 3-4. Such a scenario is shown in Figure 3. For two blobs with ρ 1 ∼ ρ 2 , γ u β u = 3 translates toΓ ≈ 3. Using Equations (B11) and (B12), we find Γ 2 Γ 1 36.
As an example, two masses of similar properties that propagate with initial Lorenz factors of Γ 1 ≈ 50 and Γ 2 ≈ 1000 would give rise to two RMSs with ∆ / < 1, which can be modelled by the KRA.