© The Authors 2026

1. Introduction

The emission mechanism responsible for the prompt phase in gamma-ray bursts (GRBs) has been a matter of discussion for many decades. The competing models can be divided into two families depending on the value of the optical depth at the emission radius towards a distant observer, τ. In the optically thin family (τ < 1), the models predict that the observed γ-rays are due to synchrotron radiation from high-energy electrons. The electrons can be accelerated either by internal shocks (Narayan et al. 1992; Rees & Mészáros 1994; Kobayashi et al. 1997; Daigne & Mochkovitch 1998) or by magnetic reconnection (Spruit et al. 2001; Drenkhahn & Spruit 2002; Uhm & Zhang 2014). More recently, synchrotron emission from protons has also been discussed (Ghisellini et al. 2020; Florou et al. 2021; Bégué et al. 2022). The optically thick family of models suggests that the observed radiation is emitted when the initially trapped photons decouple at the photosphere (τ = 1, Paczyński 1986; Goodman 1986).

The typical observed GRB spectrum is much broader compared to thermal spectra (e.g. Yu et al. 2016). The released emission must therefore be out of thermal equilibrium with the plasma for photospheric models to be viable. Since the radiation should be in thermodynamic equilibrium with the plasma close to the central engine, this requires that something disturbs the equilibrium close to the transparency radius (Rees & Mészáros 2005). Numerical simulations of GRB jets suggest that the jet is highly variable, which should lead to internal dissipation (Lazzati et al. 2009; López-Cámara et al. 2013; Ito et al. 2015; Gottlieb et al. 2019).

Shocks that occur below the photosphere in a GRB jet are mediated by radiation (Levinson & Bromberg 2008). The photon distribution downstream of a radiation-mediated shock (RMS) is very different compared to a synchrotron spectrum radiated downstream of an optically thin collisionless shock. Therefore, it is important to capture the correct shock physics if one wants to compare model predictions to data (see Levinson & Nakar 2020, for a review on RMSs). RMSs in GRBs have been studied extensively in the past (Eichler 1994; Levinson & Bromberg 2008; Budnik et al. 2010; Katz et al. 2010; Bromberg et al. 2011; Levinson 2012; Keren & Levinson 2014; Beloborodov 2017; Ito et al. 2018; Lundman et al. 2018; Samuelsson et al. 2022; Samuelsson & Ryde 2023; Wistemar et al. 2025a,b). However, the main focus in these studies has been the spectral shape rather than the time-dependent signal that one might expect from these events.

In this paper, we study for the first time the time-resolved signal expected from a GRB photosphere including subphotospheric dissipation due to RMSs. To achieve this, we used a combination of different numerical codes, as outlined in the flowchart in Figure 1. First, we used a 1D special relativistic hydrodynamical code (GAMMA, Ayache et al. 2022) to simulate an internal collision. The details regarding the hydrodynamical implementation are given in Section 2. The internal collision produces a reverse and a forward shock, and the two RMSs are modelled using the Kompaneets RMS approximation (KRA) developed in Samuelsson et al. (2022). However, as explained in Appendix A, the radiation field in the shocks are not necessarily in steady state, and we have made several modifications and improvements to the KRA to tackle the problem at hand. The modifications, outlined in Section 3, include the modelling of forward and reverse shocks, a dynamical (time-dependent) treatment of the two shocks, and a more realistic approach at low optical depths. The new version of the KRA is modelled using the code Komrad dynamical. The ejecta dynamics from the hydro simulation, together with the comoving photon distribution as a function of optical depth from the RMS modelling, are then used as input for a third code called Raylease (Alamaa 2024). As explained in Section 4, it calculates the observed time-resolved signal, accounting for a probabilistic photon decoupling as a function of both optical depth and angle.

Fig. 1. Flowchart showing the methodology of the paper. Each box shows one step of the chain, with the treated physics given in light blue and the simulation code used given in red. The final result is the time-resolved signal in the observer frame. The hydrodynamics are explained in Section 2, the RMS modelling in Section 3, and the photon decoupling in Section 4. The observed time-resolved signal is shown in Section 5.

Our results are presented in Section 5. We find that the peak energy remains constant during the pulse rise time and starts to decrease at later times. At very early times, the spectrum has a complex spectral shape, with clear signatures from both the reverse and the forward shock. Later, the spectrum above the peak is well described by a power-law with an exponential cut-off, with both the spectral index and the cut-off energy decreasing as a function of time. The low-energy part consists of a smooth curvature, which may be approximated as a double power law. The results are discussed in Section 6, with emphasis placed on some key assumptions. Finally, we conclude in Section 7.

Temperatures in this paper denoted by θ are comoving and dimensionless and are given in units of electron rest mass energy as θ = k_BT/m_ec². Similarly, photon energies denoted by ϵ are comoving and dimensionless, defined as ϵ = hν/m_ec². Energies denoted by E are in the observer frame and given in units kiloelectronvolt. The outermost material in the ejecta is referred to as the front edge, and the innermost material, closest to the central engine, is referred to as the back edge.

2. Hydrodynamical implementation

2.1. GAMMA

The dynamics of the relativistic ejecta was computed in this paper using the publicly available code GAMMA (Ayache et al. 2022). GAMMA is a moving mesh special relativistic hydrodynamical code that is well suited for the treatment of GRB outflows. The code uses an arbitrary Lagrangian-Eulerian approach, which means that the mesh edges can be moved at any velocity during the dynamical evolution. Matching the velocity of the mesh edges to that of the flow makes the code pseudo-Lagrangian, which increases the numerical resolution (Duffell & MacFadyen 2011, 2013; Duffell 2016). In 1D, the code is perfectly Lagrangian.

The code uses a finite-volume Gudonov scheme, and the fluxes at the cell interfaces are calculated using an approximate Riemann solver (HLLC; Mignone & Bodo 2006). GAMMA has been developed to follow the dynamics of internal shocks and of the jet deceleration by the ambient medium and to compute the associated radiation in the optically thin regime. It is well adapted to highly variable outflows (see e.g. Ayache et al. 2020). Therefore, we used it here to get a handle on the GRB dynamics at lower radii, starting well below the photosphere.

A very important aspect when modelling the optically thick part of a GRB jet is to allow for the feedback of the radiation on the hydrodynamical properties. Unfortunately, GAMMA does not provide such feedback. This is a delicate point that is further discussed in Section 6.4. Therein, we argue that the feedback is not important in the current framework since we are not interested in exactly how the pressure, density, and velocity varies across the ejecta at any given time, but rather we care about the average values of these quantities in the two upstreams and the downstream. These average quantities should generally be accurate since they are determined by the continuity equations (shock-jump conditions).

2.2. Initial conditions

We considered a central engine that ejects relativistic material during a lab frame duration, t_e. The variability of the jet launching process at the central engine and the interaction of the jet and the surrounding material during the early propagation leads to variable observed injected power, Ė, and observed injected mass flux, Ṁ. As a consequence, the terminal Lorentz factor, η ≡ Ė/Ṁc², also varies. The fluctuations of η leads to shocks developing inside the jet, which reprocesses the kinetic energy of the outflow into internal energy, thereby shaping the observed emission. These internal shocks can occur both above and below the photosphere depending on the profiles of Ė(t) and Ṁ(t).

How Ė(t) and Ṁ(t) vary across the jet is not known, and their profiles are likely complex; see further discussion in Section 6.1. In this paper, we focused on one single internal collision leading to a single short pulse in the light curve. The mass flux was assumed constant over the ejecta duration. The injected energy flux was assumed to be a constant low value, Ė_s, during an initial time, t_s, to smoothly increase during a transition time, t_t, and to remain constant at a high value, Ė_f, during the remaining time, t_f = t_e − (t_s + t_t). Specifically, the profile was given by

$$\begin{array}{cc}\hfill \dot{E}(t)& ={\dot{E}}_s,\mathrm{if}0\le t\le t_s\hfill \\ \hfill \dot{E}(t)& =\frac{{\dot{E}}_s+{\dot{E}}_f}{2}-\frac{{\dot{E}}_f-{\dot{E}}_s}{2}\hfill \\ \hfill & \times cos\left(\pi \frac{t-t_s}{t_t}\right),\mathrm{if}{t}_s<t<t_s+t_t\hfill \\ \hfill \dot{E}(t)& ={\dot{E}}_f,\mathrm{if}{t}_s+t_t\le t\le t_e.\hfill \end{array}$$(1)

In the equation above, t is the lab frame time that has elapsed since the front edge was emitted from the central engine.

With the mass flux constant, the terminal Lorentz factor is directly proportional to Ė. Specifying a minimum terminal Lorentz factor, η_s, fixes the constant mass flux as Ṁ = Ė_s/η_s c² and the maximum terminal Lorentz factor as η_f = η_sĖ_f/Ė_s. Since η ∝ Ė, Equation (1) implies a velocity stratification across the ejecta, which eventually leads to the formation of a reverse and a forward shock. We limited ourselves in the present study to the case where both shocks had crossed their respective masses before the ejecta reached the photosphere.

We assumed that the ejecta accelerates freely from an initial radius, R₀, where the Lorentz factor is Γ₀. The value of R₀ could potentially be comparable to a few times the Schwarzschild radius of the black hole (e.g. Iyyani et al. 2013). However, numerical simulations suggest that the free acceleration is significantly delayed to R₀ ∼ 10⁹ − 10¹⁰ cm due to high-pressure material surrounding the jet (Gottlieb et al. 2019). We assumed a pure fireball evolution, i.e. negligible magnetisation above R₀. At the start of the hydro simulation, the back edge had just reached the saturation radius, R_s, for the fast fireball as R_min = R_s(η_f)≡η_fR₀/Γ₀. Thus, the simulation start time was t_start ≈ R_min/c + t_e. The front edge was set as R_max = R_min + ct_e¹. The simulation used 3000 grid cells equally spaced in radius.

In each cell, the comoving density was set as (e.g. Piran 1999)

$$\begin{array}{c}\hfill \rho (R)=\frac{\dot{E}}{4\pi R^2{c}^3\mathrm{\Gamma }\eta },\end{array}$$(2)

the comoving pressure as

$$\begin{array}{c}\hfill p(R)=\frac{\rho c^2}{4}{\left(\frac{R}{R_s}\right)}^{-2/3}{\left(\frac{\mathrm{\Gamma }}{\eta }\right)}^{-1/3},\end{array}$$(3)

and the Lorentz factor as

$$\begin{array}{c}\hfill \mathrm{\Gamma }(R)=\frac{\eta }{1+4\frac{p}{\rho c^2}},\end{array}$$(4)

where R_s ≡ ηR₀/Γ₀ is the saturation radius for a fireball with terminal Lorentz factor η². These equations were obtained by assuming η ≫ 1, a relativistic equation of state, and an adiabatic expansion of the fireball. Equations (2)–(4) result in a cubic equation for (Γ/η)^1/3, which close to the saturation radius has the real valued solution

$$\begin{array}{c}\hfill \frac{\mathrm{\Gamma }}{\eta }=\frac{1}{3}(3-C_1^2+C_2-\frac{6C_1^2-C_1^4}{C_2}),\end{array}$$(5)

where C₁ ≡ R_s/R and

$$\begin{array}{c}\hfill C_2\equiv {(-C_1^6+9C_1^4-\frac{27}{2}{C}_1^2+\frac{3\sqrt{3}}{2}\sqrt{27C_1^4-4C_1^6})}^{1/3}.\end{array}$$(6)

The above solution is valid when $C_1<\sqrt{27/4}$ (corresponding to the regime where 4p/ρc² < 3). Equations (2), (3), (5), and (6), together with $v=c\sqrt{1-1/{\mathrm{\Gamma }}^2}$, determined the initial state in each cell. Thus, the free parameters of the problem are t_e, t_s, t_t, Ė_s, Ė_f, η_s, R₀, and Γ₀.

2.3. Boundary conditions

In this paper, we performed two simulations that differed only by their boundary conditions. The first simulation assumed that outside of the ejecta region on both sides, there existed a low density and low pressure material that propagated with a lower Lorentz factor. In practice, this was achieved by setting p_ghost/p_inside = ρ_ghost/ρ_inside = 10⁻³ and Γ_ghost = Γ_inside/2 + 0.55. Here, the subscript “ghost” refers to properties in the ghost cells and the subscript “inside” refers to quantities in the grid cell closest to the ghost cell that is inside of the simulated region. The expression for Γ_ghost assured that it was always larger than unity. We used this expression instead of a vacuum at rest since it significantly increased the computational speed. These boundary conditions were supposed to mimic a vacuum outside of the ejecta region. This choice is motivated once the jet has escaped the surrounding high density material, originating from the parent star in the case of a collapsar progenitor or the tidally disrupted neutron star(s) in the case of a merger progenitor.

The second simulation was motivated by numerical studies of collapsars, where it is commonly found that the jet at break out is highly variable (e.g. López-Cámara et al. 2014; Ito et al. 2015, 2019; Gottlieb et al. 2019, 2020). This variability would lead to multiple internal collisions occurring in the jet, with the single collision studied in this paper being one of these. This second simulation used a more appropriate periodic boundary condition, where the properties of the ghost cells were determined by the properties in the cells on the other side of the simulated region. In practice, this was achieved by assuming that ρR², pR^8/3, and Γ were periodic.

As discussed in Section 6.1, we found that the two simulations gave very similar results. Thus, we present only the results from the first simulation in Section 5.

2.4. Geometry of the ejecta

The hydrodynamical setup described in Section 2.2 generated one reverse shock propagating to the left and one forward shock propagating to the right. Here, left and right are in terms of the Lagrangian mass coordinate, m(R, t), defined as the total mass in the ejecta between the central engine and R at time t as

$$\begin{array}{c}\hfill m(R,t)={\int }_0^{R}4\pi R^2\mathrm{\Gamma }\rho dR.\end{array}$$(7)

One of the main purposes of the hydrodynamical simulation in this paper was to obtain the evolution of the shock properties as a function of radius, which served as input in the RMS modelling. In this regard, there were five different regions of interest in each snapshot of the hydrodynamical simulation. From left to right in mass coordinate, these were: i) the upstream of the reverse shock, ii) the immediate downstream of the reverse shock, iii) the far downstream region (common for both shocks), iv) the immediate downstream of the forward shock, and v) the upstream of the forward shock. Throughout the paper, hydrodynamical quantities are denoted by a subscript u in the upstream regions i) and v), by a subscript d in the immediate downstream regions ii) and iv), and by a subscript * in the common far downstream region iii). A sketch of the Lorentz factor profile across the ejecta after both shocks have formed can be seen in the top panel in Figure 2, where the five regions have been marked.

Fig. 2. Top panel: Sketch of the Lorentz factor profile across the ejecta with the positions of the five different regions mentioned in subsection 2.4 marked. The dashed vertical lines indicate the position of each region relative to the KRA zones given in the bottom panel. Bottom panel: Geometry of the KRA with both reverse and forward shocks included. Green, red, and purple colours indicate the upstream zones, the RMS zones, and the common downstream zone, respectively. The zones are coupled via source terms as indicated in the figure. The RMS and downstream zones are evolved using the Kompaneets equation (Equation (C.2)) while the upstream zones are assumed to be in a thermal Wien distribution. The relevant temperature in each zone is given. This panel is adapted from Samuelsson et al. (2022), which shows the geometry of the KRA with three zones. In between the two panels, the subscripts used to refer to the different regions and zones are shown for clarity.

2.5. Output of the hydrodynamical simulation

This subsection lists the hydrodynamical properties that are of interest for the rest of the paper. To model the evolution of the comoving photon distribution in Komrad dynamical, we needed information regarding the structure of the two shocks, the scattering timescale, and the overall cooling of the ejecta due to adiabatic expansion. The latter was quantified using the parameter φ, defined as

$$\begin{array}{c}\hfill \phi \equiv -\frac{1}{3}\frac{dlog\rho }{dlogR}.\end{array}$$(8)

In an idealised outflow where the comoving density decreases as ρ ∝ R⁻², one gets φ = 2/3.

The scattering timescale, t_sc, is important since it determines the time it takes for the radiation in each shock to reach steady state, as well as the level of thermalisation experienced by the downstream radiation field. In the observer frame, it is given by

$$\begin{array}{c}\hfill t_{\mathrm{sc}}=\frac{\mathrm{\Gamma }}{\kappa \rho c},\end{array}$$(9)

where κ is the opacity, and we assumed Γ ≫ 1. In this paper, we use κ = 0.4 cm²/g.

When the scattering time is much shorter than the dynamical time, the radiation in the RMS is in steady state as explained in Appendix A. Each steady-state RMS is characterised by three parameters, which can be taken to be the dimensionless upstream velocity as measured in the shock rest frame, β_u, the comoving upstream temperature, θ_u, and the photon-to-proton ratio, ζ, (e.g. Lundman et al. 2018). The exact calculation of these three quantities from the hydrodynamical output is detailed in Appendix B. In short, the upstream velocity was obtained from the shock jump conditions since the increase in pressure and density across the shock was known. The upstream temperature was given by the ratio of the upstream pressure to the upstream particle number density. The photon-to-proton ratio at mass coordinate m was obtained by assuming a thermodynamic equilibrium at R₀, after which ζ was calculated as a function of Ė and η. The three RMS parameters were at any given time different for the two shocks, and they varied as a function of radius for both shocks.

To obtain the time-resolved signal in the observer frame in Section 4, we needed to compute the evolution of the optical depth within the jet. Commonly, the radial optical depth from R towards a distant observer is estimated as τ = κρR/Γ (e.g. Beloborodov 2011). However, for the vacuum boundary case considered in Section 5, the ejecta was very thin, making photon escape at the front edge important. Thus, we calculated τ as (Abramowicz et al. 1991; Daigne & Mochkovitch 2002)³

$$\begin{array}{c}\hfill \tau (R,t)={\int }_R^{R_{\mathrm{esc}}}\frac{\kappa \rho }{\mathrm{\Gamma }}dR,\end{array}$$(10)

where R_esc is the radius where a radially moving photon injected at R and t reaches the front edge of the ejecta.

With the definition in Equation (10), each mass coordinate had one associated value for the optical depth at each time. However, to obtain the observed signal, we wished to associate the whole ejecta at each lab frame time with a single value of the optical depth. This value was taken to be the mass averaged value of the optical depth. This was an approximation since the optical depth actually varied significantly across the ejecta at any given time. Specifically, the front edge had a lower optical depth compared to the rest of the ejecta, which meant that in reality one might expect a thermal precursor to the main signal. Considering such effects are unfortunately outside the scope of the current paper.

Using the mass averaged value, the optical depth became a function of the lab frame time only, τ(t). Thus, we could also obtain the lab frame time when the ejecta had reached a specific optical depth, t(τ). This quantity is used in the Section 4. Each snapshot of the hydro also needed to be associated with a single radius, which was taken to be the mass averaged value over the ejecta. This was valid since the width of the ejecta, ΔR, was at all times much smaller than the current radius, ΔR/R ≪ 1.

3. The Kompaneets RMS approximation in a dynamical flow

The hydrodynamical simulation described in the previous section gave us a handle on the dynamics of the GRB ejecta. As we considered initial conditions leading to the formation of shocks at high optical depth, we needed to model the evolution of the photon distribution below, at, and above the photosphere to get accurate predictions for the observed time-resolved signal. Crucially, we needed to capture how the photon distribution changes due to the subphotospheric reverse and forward shocks, where photons can be repeatedly scattered in a steep velocity gradient, significantly altering their energies.

The Kompaneets RMS approximation (KRA, Samuelsson et al. 2022) is an efficient way to model radiation in this context, without the need for computationally expensive simulations. The KRA was developed in Samuelsson et al. (2022) assuming a steady jet (constant Γ and ρ ∝ R⁻²). Furthermore, the shock was assumed to dissipate all of its energy during one doubling of the radius. Therefore, the shock life-time was comparable to the expansion time of the jet, and potential dynamical effects on the shock itself were ignored. However, as derived in Appendix A, the radiation in the shock goes from being in a steady state to being dynamically evolving already at an optical depth of τ = 1/β_u² for a steady jet. Thus, this transition can happen far below the photosphere in non-relativistic RMSs. To account for this, as well as the fact that the ejecta properties evolved significantly, i.e. the jet was not steady, we developed Komrad dynamical, which is a numerical code that utilises the KRA to simulate RMSs in a dynamical environment. All corresponding modification are detailed in Appendix C. In the following, we start by giving a brief description of the KRA (the interested reader is referred to Samuelsson et al. (2022) for more details), after which the most significant modifications are summarised.

3.1. A brief summary of the KRA

When photons traverse an RMS, they Compton scatter in the velocity gradient across the shock. This is a Fermi-like process in which the photons gain energy on average with each scattering. For non-relativistic and mildly relativistic RMSs, the relative photon energy gain in a scattering is independent of the pre-scattered photon energy. This leads to the development of a power-law photon distribution in the shock region, extending between two characteristic energies, ${ϵ}_{\min }^{\mathrm{RMS}}$ and ${ϵ}_{\max }^{\mathrm{RMS}}$ (Ito et al. 2018; Lundman et al. 2018; Samuelsson et al. 2022).

The above process is very similar to thermal Comptonisation of photons with hot electrons, which is also characterised by Compton scattering and an energy-independent average increase of the photon energy per scattering. This is the idea behind the KRA. The KRA models thermal Comptonisation using the Kompaneets equation, which is numerically easy to solve. The RMS transition region is split into three zones, the upstream zone, the RMS zone, and the downstream zone. Thermal photons are injected from the upstream zone into the RMS zone. In the RMS zone, the photons interact with hot electrons, which increases the average photon energy. The energy gain in the RMS zone mimics the energy gain a photon would experience due to bulk Comptonisation in a real RMS. Photons have a non-zero probability per scattering time to be advected from the RMS zone into the downstream zone. Since the escape probability is energy independent, a power-law spectrum develops in the RMS zone between two characteristic energies, ${ϵ}_{\min }^{\mathrm{KRA}}$ and ${ϵ}_{\max }^{\mathrm{KRA}}$.

The steady state RMS zone spectrum in the KRA is characterised by three parameters: the comoving upstream zone temperature θ_u, K⁴, the effective temperature in the RMS zone, θ_r, and the Compton y-parameter of the RMS zone, y_r. Here, the subscript u when used for KRA quantities refers to the upstream zone, corresponding to regions i) and v) in the hydro notation, and the subscript r refers to the RMS zone, which is situated in between the upstream and the immediate downstream regions, i.e. between regions i) and ii), and iv) and v) in the hydro notation. The subscript * is used when referring to quantities in the downstream zone, which corresponds to region iii) in the hydro notation.

The photons in the upstream zone are assumed to be in a thermal Wien distribution at temperature θ_u, K. Since the average photon energy in a Wien distribution equals three times the temperature, the lower characteristic energy is given by ${ϵ}_{\min }^{\mathrm{KRA}}=3{\theta }_{u,\mathrm{K}}$. The maximum energy in the RMS zone is reached when the energy gain in a scattering event is balanced by the energy loss due to electron recoil. This gives ${ϵ}_{\max }^{\mathrm{KRA}}=4{\theta }_r$ (Rybicki & Lightman 1979; Samuelsson et al. 2022). The parameter y_r is a measure of the average photon energy increase within the RMS zone. With ${ϵ}_{\min }^{\mathrm{KRA}}$ and ${ϵ}_{\max }^{\mathrm{KRA}}$ set by the other two parameters, the effect of y_r is to regulate the slope of the power-law segment in between the two characteristic energies.

To get the KRA to produce similar spectra to an RMS, one requires ${ϵ}_{\min }^{\mathrm{KRA}}={ϵ}_{\min }^{\mathrm{RMS}}$, ${ϵ}_{\max }^{\mathrm{KRA}}={ϵ}_{\max }^{\mathrm{RMS}}$ and that the average energy of the radiation advected into the downstream is similar in both frameworks. Via analytical arguments and an empirical study, an accurate conversion between the KRA parameters and the physical RMS parameters was found in Samuelsson et al. (2022) (for completeness, it is also given in Appendix D.1). With this conversion, the steady-state output of the KRA agreed very well with the output from a special relativistic, Lagrangian radiation hydrodynamical code (radshock, Lundman et al. 2018) (see Figures 2 and 3 in Samuelsson et al. 2022). However, the simulation time was drastically decreased by as much as 4–5 orders of magnitude, highlighting the strength of the KRA. The parameter correspondence is one-to-one, which means that each steady state RMS spectrum corresponds to one unique set of KRA parameters.

The KRA is applicable to RMSs in weakly magnetised environments, where the incoming upstream plasma velocity satisfies γ_uβ_u ≲ 3.5 as measured in the stationary shock frame⁵. As expected for RMSs in GRBs, the shock is assumed photon rich, which means that the main source of photons is advection of existing upstream photons rather than photon generation inside the shock and the immediate downstream (Bromberg et al. 2011).

3.2. Major modifications in KRA dynamical

In this subsection, we summarise some of the most significant modifications to the KRA in a dynamical flow. All of the details are given in Appendix C.

The KRA has been updated such that it can model two shocks if both a reverse and a forward shock are present. In case of two shocks, the KRA consists of five different zones that are coupled via source terms. Radiation flows from the upstream zones, via the two shock zones, to the common downstream zone. The common downstream is treated as a single fluid with one collective electron temperature. This temperature is not a free parameter. Rather, the temperature is assumed to equal the downstream Compton temperature, θ_*, C, defined as the temperature the electrons get in Compton equilibrium with the current downstream radiation field. This is a good approximation since the electrons are tightly coupled to the radiation and their thermalisation timescale is incredibly short compared to the photon scattering time (Lundman et al. 2018). A schematic of the geometry with the five zones, and their positions relative to the hydrodynamical ejecta, is given in Figure 2.

The previous version of the KRA assumed that the jet was steady and that the upstream properties did not vary significantly over the shock life-time. In practice, this meant that β_u and ζ were both constant, and φ = 2/3 such that all photon energies and the upstream temperature decreased as ϵ, θ_u, K ∝ R^−2/3. In the modified version, β_u, ζ, θ_u, and φ can all vary with radius. In the current paper, the evolution of each parameter was given by the hydrodynamical simulation, as explain in Section 2.5.

When an RMS reaches lower optical depths, it broadens (Levinson 2012). An analytical description of this is given in Appendix A but one can get an intuitive idea for the origin of this broadening from the following argument. For an RMS in steady state, there is a necessary number of scatterings that the photons have to make in order to mediate the shock and stop the incoming flow (as seen in the shock stationary frame). As the scattering length increases at lower optical depths when the plasma density decreases, the shock transition region must increase. The dynamical version of the KRA accounts for this broadening by varying the source terms as a function of radius such that the RMS zone always contains the appropriate number of photons.

The upstream zone contains a finite amount of photons, N_u, tot. When this zone is depleted, the RMS dissolves. In Komrad (the name of the original simulation code), the photons in the RMS zone were instantaneously removed following the upstream depletion. However, if this happens at low optical depths, the shocks can contain a large number of photons due to the aforementioned RMS broadening (Levinson 2012). This progressive dissolution of the RMS is taken into account in Komrad dynamical: the photon distribution in the RMS zone continues to evolve even when the upstream zone is depleted, while the remaining photons are slowly advected from the RMS into the downstream zone.

3.3. Input and output

The input parameters required to run Komrad dynamical is θ_u and N_u, tot for the upstream zones, β_u and ζ for the RMS zones, and φ and t_sc for both the RMS zones and the downstream zone. All of these parameters were obtained from the hydro simulation as explained in Appendix D.

The output of Komrad dynamical is the photon spectral number density as a function of radius, N_ϵ(R), in the five different zones (see Appendix C for details). The photon spectral number density, together with the information regarding the dynamics from GAMMA, such as τ(R) and Γ(R), was used to obtain the time-resolved signal as explained in the next section.

4. Observed time-resolved signal

Photons make their last scattering at different optical depths and different angles compared to the line-of-sight due to the probabilistic nature of the decoupling process. Thus, the signal at any given observer time is a superposition of comoving spectra at different optical depths and different Doppler boosts (Pe’er 2008; Beloborodov 2011). The first photons to reach the observer are those that decouple at high optical depths and at the line-of-sight. These are followed by photons emitted at progressively lower optical depths (larger radii) and larger angles (Pe’er & Ryde 2011). If the comoving photon field evolves significantly across the region where photons make their last scattering, the observed signal will evolve as a function of time (Alamaa 2024).

The photon spectral flux, F_ϵ(t_obs), at observer time t_obs was obtained by integrating the comoving photon distribution over all optical depths and angles, accounting for the probability of a photon making its last scattering at a given optical depth and angle, as well as the predicted arrival time as (see e.g. eq. (A5) in Alamaa 2024)

$$\begin{array}{cc}\hfill F_{ϵ}(t_{\mathrm{obs}})& =\frac{1+z}{4\pi d_l^2}{\int }_0^{\infty }d\tau {\int }_{-1}^{1}d{\mu }^{\prime }{\int }_{R_{\mathrm{back}}}^{R_{\mathrm{front}}}dR\hfill \\ \hfill & \frac{f(\tau ,{\mu }^{\prime })}{\mathrm{\Delta }R}{ϵ}^{\prime }{N}_{\mathrm{tot},ϵ}^{\prime }\delta (\frac{t_{\mathrm{obs}}}{1+z}-t_{\delta }).\hfill \end{array}$$(11)

In the equation above, z is the redshift, d_l is the luminosity distance, μ′=cos(θ′), with θ′ being the comoving angle between the radial direction and the line-of-sight, and f(τ, μ′)dτdμ′ is the probability for a photon to make its last scattering between optical depth τ and τ + dτ and angle μ′ and μ′+dμ′. The probability distribution f(τ, μ′) was derived in Beloborodov (2011) and rewritten as a function of optical depth in Samuelsson & Ryde (2023)⁶. R_front and R_back are the radial position of the front and back edge of the ejecta, respectively, both evaluated at time t(τ), and ΔR = R_front − R_back is the instantaneous ejecta width. The photon spectral number density N′_tot, ϵ is that of the whole ejecta (the sum of the photon number densities in all five KRA zones) and it is evaluated at an energy $ϵ\prime =(1+z)ϵ/\mathcal{D}$, where $\mathcal{D}=\mathrm{\Gamma }(1+\beta \mu \prime )$ is the Doppler factor. The time t_δ in the Dirac delta function is given by

$$\begin{array}{c}\hfill t_{\delta }=t(\tau )-\frac{R}{c}\left[\frac{{\mu }^{\prime }+\beta }{1+\beta {\mu }^{\prime }}\right].\end{array}$$(12)

To solve Equation (11), we used the code Raylease developed in Alamaa (2024) (see Appendix A therein). The parameters ΔR, 𝒟, β, and t(τ) were readily taken from the hydrodynamical simulation. The Lorentz factor Γ and velocity β used were those for the downstream region, which contained essentially all of the radiation at the relevant optical depths where the photons escaped. The comoving photon distribution N_tot, ϵ was obtained from the KRA.

5. Results

The results presented in this section are for the simulation using the vacuum boundary conditions, as explained in Section 2.3. The results from the simulation using the periodic boundary conditions were very similar.

5.1. Ejecta evolution

In Figures 3 and 4, we show the evolution of the physical conditions in the ejecta. The parameters for the run are given in Table 1. Figure 3 shows the Lorentz factor (top) and density (bottom) at different optical depths. The distributions at the start of the simulation (τ ≈ 430) are given by the black dashed lines. At the start of the simulation the ejecta was still accelerating. Notably, the initial Lorentz factor in the fast part of the ejecta was less then half of its terminal value, η_f = 400. The acceleration was gradual and the terminal value was never quite reached due to presence of the shocks. A maximum value of Γ ≃ 370 was reached at the reverse shock crossing. The forward shock was weaker and it reached the edge of the ejecta first. Once the forward shock had crossed, the downstream region accelerated freely and the average Lorentz factor increased at the expense of the ejecta internal energy.

Fig. 3. Lorentz factor (top) and the comoving density (bottom) at different optical depths as obtained from the hydrodynamical simulation using the vacuum boundary conditions. The dashed black line gives the corresponding profile at the start of the simulation. The comoving photon distributions in the five different zones of Komrad dynamical at the corresponding times are given in Figure 4. The simulation parameters are given in Table 1.

Fig. 4. Comoving photon distribution in the five different zones as obtained by Komrad dynamical. The line coding is shown by the legend in the top-left panel, where RS and FS stand for reverse shock and forward shock, respectively. The four different panels show the distributions at different optical depths, corresponding to the optical depths in Figure 3. The observed photon distribution consists of a superposition of comoving spectra, leading to a smoothening and broadening effect on the observed spectrum. The simulation parameters are given in Table 1.

The density (ρR²) initially increased drastically in the intermediate region between the slow and the fast part due to adiabatic compression. The compression led to an increase in temperature and pressure, and the shocks formed once the immediate downstream pressure had become comparable to the incoming ram pressure. As the left-hand side started off hotter than the right-hand side, the reverse shock formed slightly earlier than the forward shock. The density continued to increase after the shocks had been launched albeit more slowly. When the forward shock had crossed, the acceleration of the downstream led to a decrease in the density, which in turn led to a quick drop in the optical depth. This means that there can be a tendency for RMSs to finish their dissipation close to the photosphere, see further discussion in Section 6.2.

Figure 4 shows the comoving photon distribution in the five zones at the corresponding optical depths, as obtained by Komrad dynamical. In the first snapshot, most of the photons were situated in the two upstreams zones. Bulk Comptonisation in the two shocks had just started and the photon distributions in the RMSs were still quasi-thermal. At optical depth τ = 30, the reverse and forward shocks had developed further, leading to a non-thermal shape. There were still a significant fraction of photons in the two upstreams zones, leading to a thermal bump around ϵ ≳ 10⁻³. The downstream spectrum had begun to broaden but was still quite narrow.

At τ = 10, the reverse shock had developed a power-law from bulk Comptonisation that stretched more than two orders of magnitude in energy. Its upper cut-off energy continued to increase, which was due to two reasons: i) as the density decreased, the shock front sped up (Sakurai 1960), which increased the maximum photon energy in the shock (e.g. Samuelsson et al. 2022). ii) The photons had a progressively harder time to mediate the shock when the ejecta approached the photosphere. As the scattering frequency decreases, the relative energy transfer per scattering must increase for the radiation to decelerate the incoming flow. Since the maximum energy is proportional to the relative energy gain per scattering (Rybicki & Lightman 1979; Samuelsson et al. 2022), this further increased the upper cut-off energy in the RMS. This latter behaviour is apparent from Equation (C.11) and became important when the RMS entered the dynamical stage.

That the maximum energy in the RMS increases at lower optical depths is in direct contrast to the behaviour of a photon poor RMS at shock breakout in a stellar wind (Ioka et al. 2019; Ito et al. 2020). In this case, the escaping radiation at low optical depths leads to a higher compression behind the shock, which in turn enhances photon production via bremsstrahlung emission. Thus, more photons share the dissipated energy, which decreases the temperature. However, in the very photon rich regime considered here (ζ ∼ 10⁵), photon production via bremsstrahlung would remain negligible even at high radiative losses (Ioka et al. 2019). Thus, we obtained the opposite trend for the maximum energy.

From the figure it is clear that the forward shock upstream had already been depleted at this stage, and the forward RMS had begun to dissolve. Note that the upstream was depleted before the hydro shock reached the front edge, as can be seen from the τ = 10 profile in Figure 3. This was because the RMS transition region became broad at low optical depth, and the front of the large transition region reached the front edge of the ejecta earlier than the thin hydro shock. The reverse shock upstream was not yet depleted. Its temperature had dropped by a factor of ∼3 since τ = 50 due to adiabatic cooling. The downstream region contained about half of the total photon population and consisted of a broad spectrum. The broad spectrum was partly due to the increasing scattering time, leading to a less efficient thermalisation process, and partly due to the broad reverse shock injecting both lower and higher energy photons than before.

At the photosphere (τ = 1), the downstream contained ≈95% of all photons. It had a smooth, soft curvature below the peak energy at ϵ_p ≈ 4 × 10⁻³ and a clear high-energy component. The high-energy part exhibited a slight spectral hardening before it cuts off sharply, which is a characteristic signature of a photon distribution in the “fast Compton regime” (Alamaa 2024).

5.2. Light curve

The observed light curve for a redshift of z = 1 is shown in the upper-right panel in Figure 5. The grey dashed line shows the bolometric flux. The blue and red solid lines shows the light-curves in the Fermi NaI (8 − 1000 keV) and Fermi BGO (0.4 − 40 MeV) band, respectively, the Fermi NaI and BGO detectors being part of Fermi Gamma-ray Burst Monitor (GBM, Meegan et al. 2009). The energy output in the NaI and BGO bands were comparable and the bolometric flux was completely dominated by the emission within these bands (note that F_bol < F_NaI + F_BGO due to the overlap between the detectors). At peak flux, the burst reached ∼5 × 10⁻⁶ erg s⁻¹ cm⁻². This is about an order of magnitude higher than the median peak flux of Fermi GRBs (compared to the BEST sample in Poolakkil et al. 2021). Thus, these parameters produced a fairly bright GRB at the current redshift of z = 1.

Fig. 5. Observed spectral evolution (top left), light curve (top right), count spectrum of a Band fit to the time-integrated spectrum (bottom left), and the evolution of the best-fit Band parameters (bottom right). The dashed grey line in the top-left panel shows the time-integrated spectrum taken as the average over the first second of observations, and purple shading shows the energy range of the GBM, with the darker colour indicating the region with the highest effective area. The dot-dashed red line in the light curve shows the light curve in the BGO band scaled to the peak of the NaI light curve to highlight the difference between the two. The red, blue, cyan, and green data points in the count spectrum show the three NaI detectors and the one BGO detector used for the fit, respectively. The best-fit peak energy is given in units Ep, 2 = Ep/100 keV, and the mock dataset for the fit had a S/N = 50. The best-fit parameters with errors for the time-integrated spectrum were α = −0.96 ± 0.06, β = −2.23 ± 0.09, and Ep = 156 ± 16 keV. All errors given represent 1σ. The simulation parameters are given in Table 1.

Integrating the bolometric flux, we obtained the radiated gamma-ray efficiency as

$$\begin{array}{c}\hfill {\eta }_{\gamma }=\frac{E_{\mathrm{rad}}}{E_{\mathrm{tot}}}=\frac{\int 4\pi d_l^2{F}_{\mathrm{bol}}\frac{d{t}_{\mathrm{obs}}}{(1+z)}}{{\int }_0^{t_e}\dot{E}dt}\approx 23\%,\end{array}$$(13)

where the numerator has been integrated over the first second of observation. The value found is quite high. RMSs are very efficient at converting kinetic energy into radiation energy since almost all dissipated energy is taken by the photons. This is different compared to collisionless shocks where the radiation ends up with no more than the fraction of dissipated energy given to electrons, which is usually of the order ∼0.1 (Sironi & Spitkovsky 2011). However, in contrast to collisionless shocks, the radiation downstream of an RMS suffers from adiabatic degrading, which can severely decrease the radiated gamma-ray efficiency. That the reverse shock terminates at low optical depth was the reason for the high η_γ in this case. As a consistency check, we found that kinetic energy at the photosphere E_kin(R_ph) = Γ(R_ph)M_totc² was ≈79% of the total energy, where M_tot = Ṁt_e is the total ejected mass, implying that the kinetic plus radiated energy was equal to E_tot within 2%. This means that we did not artificially inject too much energy into the photons during the RMS modelling.

The duration of the pulse was short, with a duration comparable to the ejecta duration, t_e. It was of the same order as the observed geometrical timescale at the photosphere t_geo = (1 + z)R_ph/2cΓ², which, for the obtained R_ph ≈ 3 × 10¹³ cm and Γ(R_ph)≈130, gives t_geo ≈ 0.06 s. In this paper, we focused on a single internal collision. A more complex environment would lead to a more complex light curve consisting of many, potentially overlapping pulses, (see further discussion in Section 6.1). Furthermore, additional surrounding material can keep the downstream region compressed, which increases the photospheric radius and the observed geometrical timescale. While the aforementioned complexities can moderately increase the pulse duration, many GRBs have pulse durations that are much longer than what we found here (Golkhou & Butler 2014). These pulses require emission at larger radii, potentially originating from shocks above the photosphere.

The light purple dotted line shows the hardness as a function of time, calculated as the ratio between the flux in the 100–300 keV band to the flux in the 50–100 keV band. It rose by a factor of four from very early times to the light curve peak, after which it decayed smoothly. Although the rise in the hardness ratio is drastic, it would likely be difficult to observe as the GRB was very faint at that time. The decay happened once photons from radii above the reverse shock crossing radius started to reach the observer. When the injection of high-energy photons from the reverse shock ceased and adiabatic cooling started to dominate in the comoving frame, the hardness ratio started to decrease (see further discussion in Section 5.4).

The dot-dashed red line shows the BGO flux scaled to peak at the NaI maximum flux to highlight the differences between the two light curves. The light curve was very similar in the two bands up until the peak but decreased more rapidly in the BGO band. The explanation for the faster decay in the BGO band is the same as that for the decreasing hardness ratio.

5.3. Time-integrated spectrum

The observed time-integrated spectrum, taken as the average over the first second of observation, is shown by the dashed line in the top left panel of Figure 5. Qualitatively, it was quite similar in shape to the total comoving spectrum at the photosphere as shown in the bottom right panel of Figure 4. However, the superposition of radiation emitted at a range of optical depths (different spectral shapes) and angles (different Doppler boosts) led to a time-integrated spectrum that was broader and softer at low energies and had a more extended high-energy power law. The spectrum peaked at E_p ≈ 200 keV, which is a typical value for GRBs (Poolakkil et al. 2021). At low energies, the spectrum had a smooth curvature from E_p down to ∼3 keV, below which the spectrum became very steep (tending asymptotically towards α ∼ 0.4, Beloborodov 2010; Lundman et al. 2013). Above the peak, the spectrum exhibited a high-energy power law extending for ∼1.5 orders of magnitude with a slope β ≈ −2.3 cutting off at ∼5 MeV.

To get a more qualitative idea of how such a spectrum compared to actual observations, we used a similar analysis to Samuelsson & Ryde (2023), Alamaa (2024). The time-integrated spectrum was forward-folded through a detector response matrix. This generated a mock dataset that accounted for both detector and background effects. The background was obtained by using the response file from a real GRB, in this case GRB 150213A, and the background strength was determined by setting a user specified signal-to-noise ratio (S/N). The mock dataset could then be fitted with any standard GRB function such as the Band function (Band et al. 1993). In this paper, we used the Fermi Gamma-ray Burst Monitor (GBM, Meegan et al. 2009) as the detector. The procedure automatically accounted for the detector effective area, which meant that the fit could only consider data within the GBM energy sensitivity window (8 keV to 40 MeV, Meegan et al. 2009).

A Band function fit to the time-integrated spectrum is shown in the bottom-left panel in Figure 5. The top shows the count spectrum and the bottom shows the residuals. The mock dataset had S/N = 50 and the best-fit parameters for the fit were α = −0.96 ± 0.06, β = −2.23 ± 0.09, and E_p = 156 ± 16 keV where α, β, and E_p are the low-energy spectral index, the high-energy spectral index and the peak energy, respectively.

It is clear from the residuals in the bottom panel that the Band function was a good fit when S/N = 50. However, the time-integrated spectrum had more complexity below the peak than a single power law. To test this, we generated a new mock dataset with S/N = 300. The Band function fit to this dataset is shown in the top panel of Figure 6. The evident oscillations in the residuals in this case indicated that the Band function was not a good fit at this high S/N. In contrast, a phenomenological function that allows for two breaks below the peak (double smoothly broken power law (2SBPL), Ravasio et al. 2018) fitted the spectrum very well. Using the Akaike Information Criterion (AIC, Akaike 1974) to evaluate the performance of each model, we found ΔAIC = AIC_Band − AIC_2SBPL = 114, indicating that the 2SBPL was greatly favoured by the data. In contrast, we found no statistically significant improvement using the more complex 2SBPL function when S/N = 50.

Fig. 6. Comparison of a Band function fit (top) and a 2SBPL fit (bottom) to the time-integrated spectrum with S/N = 300. The red, blue, cyan, and green data points in the count spectrum show the three NaI detectors and the one BGO detector used for the fit, respectively. It is clear from the uneven residuals in the top panel that the Band function did not produce a satisfactory fit. The 2SBPL function performed much better in this region of high S/N, being statistically preferred with ΔAIC = 114. The best-fit Band parameters were α = −0.92 ± 0.02, β = −2.33 ± 0.02, and Ep = 152 ± 4 keV, and the best-fit 2SBPL were α1 = −0.69 ± 0.06, α2 = −1.45 ± 0.04, β = −2.49 ± 0.04, Ep = 178 ± 7 keV, Eb = 27 ± 3 keV.

It is interesting to note that the values for the two low-energy slopes were very close to the theoretical values expected in synchrotron emission models, α₁ = −0.69 ± 0.06, α₂ = −1.45 ± 0.04. As pointed out already in Samuelsson & Ryde (2023), there is no a priori reason for why this would be the case. Rather it is a consequence of fitting a function that allows for two power laws below the peak to a curved spectrum that, when fitted with the Band function, produces α ∼ −0.9. Then, for most fits, one gets α₂ < α < α₁, which may by chance give α₂ ≈ −3/2 and α₁ ≈ −2/3. This shows that one must be cautious in the interpretation of the underlying physics based on empirical fits.

5.4. Time-resolved spectrum and spectral evolution

The instantaneous spectral flux at four different observer times are shown in the top-left panel in Figure 5. It is evident that the average energy of the observed radiation decreased with time. Photons that make their last scattering at small radii outrun the ejecta, travelling with a slightly larger velocity compared to the bulk plasma. This means that photons from deeper layers and on the line-of-sight are the first to arrive at the observer. These photons have a high Doppler boost and have suffered less from adiabatic cooling, and therefore they have a higher average energy compared to the radiation that arrives at later times.

The very early observation (t_obs = 0.15 s, yellow line) had a more complex spectral shape compared to the others, consisting of a main component peaking at ∼300 keV with several additional bumps. The radiation observed at early times during the light curve pulse rise come from the same local region in the jet (τ ∼ 10, μ′∼1, Alamaa 2024). Thus, the observed spectrum at this time should be very similar to the Doppler boosted comoving spectrum at τ = 10. This can be confirmed by comparing it to the total comoving photon distribution at τ = 10 shown in the bottom-left panel in Figure 4. The comparison shows that the main component was due to shocked photons sitting in the downstream region, while the highest energy component originated from the reverse shock.

The peak flux spectrum (t_obs = 0.3 s, orange line) was quite similar to the time-integrated spectrum, with a slightly higher peak-energy, a harder low-energy slope, and a curvier high-energy part with a bump-like signature at 5 MeV. At later times, the peak energy decreased and the low-energy part became softer. The high-energy cut-off decreased due to adiabatic cooling, but also due to the depletion of high-energy photons in the comoving frame. Once the injection of fresh photons from the reverse shock ceased, the high-energy part quickly downscattered on the colder electrons, leading to a decrease of the cut-off energy in the comoving frame.

To get a more clear view of the parameter evolution across the GRB, we forward-folded twenty instantaneous spectra through the Fermi GBM response matrix and fitted the mock datasets with the Band function. The resulting best-fit parameters and their uncertainties are shown in the bottom-right panel of Figure 5. The S/N was set to be 300 for the peak flux bin, while the other bins had their S/N determined by their relative brightness. The high maximum S/N was chosen such that the fitted parameter evolution could be seen more clearly. However, we note that a time-resolved analysis with twenty bins of such high quality in a short pulse similar to the one studied here would be extremely rare.

The fitted values are compared to theoretical estimates showed by dashed lines. The theoretical peak energy, E_p, th, was taken to be the maximum of the ϵF_ϵ spectrum. The theoretical indices were estimated as α, β = dlog F_ϵ/dlog ϵ − 1 evaluated at 0.1E_p, th and 10E_p, th for α and β, respectively.

The fitted low-energy index decreased from α ≈ −0.5 to α ≈ −1.0 between 0.2 s and 0.4 s. After 0.4 s, it remained ∼ constant although all fitted values showed a large scatter at t_obs > 0.5 s due to the low S/N. Comparing with the theoretical estimate, we can see that the fit generally underestimated the slope at 0.1E_p. This was because the slope just below the peak was softer and since it was the brightest part of the spectrum, it dominated the fit. At late times, the hardest low-energy curvature was below the GBM detector window, which also led to the fit underestimating the value of the slope. This tendency becomes more and more pronounced as the peak energy decreases further (Acuner et al. 2019).

The fitted values of the peak energy increased slightly during the pulse rise time. However, it did not represent the real peak energy very well. This was due to the complex spectral shape that did not resemble a Band function at early times. The real peak energy had a high (∼300 keV), constant value during the pulse rise, after which it started decreasing due to adiabatic cooling. That the peak energy was not decreasing during the pulse rise time was because of the existence of the reverse shock at low optical depths. Unless there is some ongoing dissipation in the decoupling region (τ < 10), the peak energy decreases throughout the pulse (Alamaa 2024).

The fitted values of the high-energy index decreased during the pulse rise time from β = −2 to β = −2.4 and remained ∼constant thereafter, albeit with a large scatter. Although subtle, such a signature is quite revealing and arises from photons with high comoving energies being downscattered in the comoving frame. Since the Compton cooling is proportional to the photon energy, these high-energy photons are the only ones that have time to be downscattered in the decoupling region. Energy transfer via Compton scattering is only effective during the pulse rise time before adiabatic cooling becomes the dominant energy loss channel (Alamaa 2024). One of the ways that this can manifest itself in the data is by a rapid decrease of the high-energy index during the pulse rise followed by a constant value during the pulse decay.

6. Discussion

6.1. Single pulse vs highly variable light curve

In this paper, we focused on a single internal collision with an idealised initial Lorentz factor profile. This simple setup is unlikely to be realistic for a real GRB. However, investigating a single collision is a good starting point before moving on to more realistic, highly variable ejecta.

The single collision considered here produced a light curve consisting of a ∼0.1 s duration pulse. This pulse could be one building block within a complicated, highly variable longer GRB light curve. Numerical simulations of collapsar GRBs show that multiple shocks usually form within the relativistic outflow (“shock train”, e.g. López-Cámara et al. 2014; Ito et al. 2015, 2019; Gottlieb et al. 2019, 2020). Specifically, Gottlieb et al. (2019, 2020) found that mixing around the collimation shock can result in strong short term variability for the terminal Lorentz factor, with typical values close to what we considered here. This variability led to internal collisions, some of which occured below the photosphere (Gottlieb et al. 2019).

If the single collision considered here occurred within a broader variable ejecta with multiple other shocks, the periodic boundary condition used by the second simulation performed are more appropriate (see Section 2.3). Remarkably, we found that the results from that simulation were very similar to those presented for the vacuum boundary case in Section 5⁷. The boundary conditions only became important at low optical depth when the shocks had crossed the considered region of the ejecta, and therefore, they had little effect on the radiation. Furthermore, the optical depth itself was not much affected either, since the previous high-velocity material had left a low-density cavity in its wake, meaning that photons made their last scattering at roughly the same radii as in the vacuum case. We conclude that the results presented in Section 5 can be applied to short pulses in complex light curves as well, while long duration pulses may require shocks occurring above the photosphere. Note that the spectral evolution obtained here would then be for a single pulse within the light curve; the complex GRB as a whole can have a vastly different evolution (see e.g. Ito et al. 2019).

The light curve shown in Figure 5 could also constitute a whole short GRB, as expected in a merger event. For instance, Pais et al. (2024) simulated short GRBs by injecting high power jets into a realistic neutron star merger environment (obtained from a 1 s long 3D GR neutrino-radiation MHD simulation, Kiuchi et al. 2023). In one of their simulations, they found that the jet material clustered in velocity space around Γ ≈ 40 and Γ ≈ 400 at the time of jet break out. The slower material consisted in part of jet material shocked by the jet head reverse shock as it propagated through the neutron star debris, while the faster material, situated further inwards, consisted of unshocked material. This would lead to the formation of shocks at a larger distance and gives some support to the Lorentz factor profile used in this paper.

6.2. Non-thermal photospheric emission from RMSs

Dissipative photospheric models require that the dissipation occurs not too deep below the photosphere for the emitted spectrum to be non-thermal (outside of the Wien zone, τ ≲ 100, Beloborodov 2013). In the current framework, there are two physical reasons for why RMSs that start at high optical depth will have a tendency to finish their dissipation close to the photosphere. The first is that while both forward and reverse shocks exist, the downstream stays compressed, keeping the optical depth high. Once either shock has crossed, the high internal energy in the downstream region leads to an acceleration of the ejecta (a second fireball), which rapidly decreases the comoving density and, subsequently, the optical depth. Thus, it is plausible that the transition to transparency happens quickly after the RMS has crossed.

The second is that number of photons in the RMS transition region increases when the shock reaches low optical depths (Levinson 2012, see also Appendix A). Indeed, it is clear from the approximate expression in Equation (A.5) that the photon number increases exponentially when $\tau /{\mathcal{N}}_{\mathrm{sc}}\ll 2$, where ${\mathcal{N}}_{\mathrm{sc}}$ is the average number of scatterings a photon performs in the RMS. This means that the shock front effectively “swallows” the upstream when the shock comes close to the photosphere. To what extent this latter phenomenon would occur in reality needs to be checked with more detailed numerical simulations.

6.3. High-energy power law

The time-integrated spectrum had a high-energy power law that extended for ∼1.5 orders of magnitude before it cut off at ∼5 MeV. It was generated by the injection of high-energy photons from the shock into a region of relatively cold plasma, where they suffered Compton losses (Alamaa 2024). However, it is not straight forward in photospheric models to get a high-energy power law extending beyond ∼Γm_ec² due to the high γγ-opacity (although see Beloborodov 2010). Thus, it is difficult to account for the sample of GRBs that have a power law extending above several hundred MeV as those presented in Ravasio et al. (2024), Macera et al. (2025). Potentially, such high-energy photons could be generated at the transition from an RMS to a collisionless shock around the photosphere.

6.4. Neglecting radiation feedback

In this paper, we have used a pure hydrodynamical simulation to get a handle on the GRB ejecta dynamics. However, below the photosphere, the photons dominate the pressure, meaning that the hydrodynamical profile will be determined by the radiation. For instance, the detailed shock structure is dependent on the photon scattering length and it can only be accurately determined if radiation feedback on the hydrodynamics is included. However, in this paper, the hydrodynamical simulation was used to obtain general properties regarding the shocks, such as the compression, the downstream enthalpy, and the downstream Lorentz factor. Such global quantities are determined by the shock jump conditions, which must be satisfied regardless of the shock microphysics. The photon distributions in the two RMSs were then modelled by the KRA, which has been compared to detailed numerical simulations where radiation feedback was included with good results in Samuelsson et al. (2022).

The above described picture changes at low optical depths. Close to the photosphere, photons decouple, which means a loss of internal energy that is currently unaccounted for. Indeed, photon escape changes the dynamics of RMSs as shown in Ito et al. (2020, 2026). Here, both shocks had crossed at moderately thick optical depths (τ > 7), where photon escape is not yet highly significant (Beloborodov 2011). Thus, we argue that the main conclusions are unaffected by this approximation.

6.5. Klein-Nishina and pair production

Klein-Nishina (KN) effects and pair production becomes important once the photon energies in the comoving frame become comparable to the electron rest mass. Both of these effects were neglected here. Including KN effects, the spectrum should soften above ϵ ≳ 0.1 (Ito et al. 2018; Lundman et al. 2018), meaning that the simulation overestimated the number of photons at the highest energies. This would have in turn affected the detailed shape of the high-energy part of the spectrum. For instance, the high-energy bump visible at ≳10 MeV at t_obs = 0.15 s in Figure 5 would have become less pronounced or disappeared. Furthermore, there would have been a precursor of high-energy photons in the light curve at energies above $\sim 10\frac{{\mathrm{\Gamma }}_2}{(1+z)}$ MeV, where Γ₂ = Γ/100. However, KN effects would likely not have had a big impact on the estimated best-fit parameters, which were determined by the spectral shape around the comoving peak energy at ϵ_p ≈ 4 × 10⁻³ where KN effects are negligible.

RMSs can lead to violent pair production if the shock velocity is high (Budnik et al. 2010; Katz et al. 2010). In the current simulation, very few photons had energies high enough for pair production. In addition, this number was overestimated as KN effects were ignored. Furthermore, the spectrum of the reverse shock at τ = 10 was quite similar to a detailed RMS simulation performed in Ito et al. (2018) (Figure 4, bottom-left panel, similar maximum energy and similar extension above ϵ > 1), which was found to produce no pairs in the RMS transition region. We conclude that including KN effects and pair production would not have altered the main conclusions of the paper.

7. Conclusion

In this paper, we have used a series of three numerical codes to go from initial parameter distributions deep below the photospheric radius in a GRB jet all the way to the best-fit parameter evolution of the time-resolved observed signal (see Figure 1 for a flowchart of the methodology). A 1D special relativistic hydrodynamical simulation code (GAMMA, Ayache et al. 2022, Section 2) was used to model an internal collision that generated a reverse and a forward shock below the photosphere. The simulation was set up such that both shocks finished their dissipation in the optically thick regime close to the photosphere. The evolution of the Lorentz factor and the comoving density as a function of optical depth was shown in Figure 3.

Since the shocks were subphotospheric, they were mediated by radiation, and thus a designated RMS model was needed. For this purpose, we used the Kompaneets RMS Approximation (KRA) developed in Samuelsson et al. (2022). However, to be able to model long-lasting, dynamically evolving shocks, and to be able to account for the contribution from both the reverse and the forward shock, several additions and modifications had to be made. This led to the development of the dynamical version of the KRA, summarised in Section 3 and described in detail in Appendix C (see Figure 2 for a schematic of the geometry). The dynamical version of the KRA was modelled by the simulation code Komrad dynamical and the obtained comoving photon distribution as a function of optical depth was shown in Figure 4.

Finally, the time-resolved signal was obtained using Raylease (Alamaa 2024). It used the ejecta dynamics from GAMMA and the subphotospheric photon distribution from Komrad dynamical as input to generate the spectrum as a function of time in the observer frame. Raylease accounts for the probabilistic nature of the photosphere, including emission from different optical depths and angles to the line-of-sight.

The results were shown in Figures 5 and 6 and can be summarised as follows:

1. The light curve consisted of a single pulse with a sharp peak and had a duration that was comparable to the geometrical timescale at the photosphere, which in this case was ∼0.1 s. It could either constitute a whole GRB in itself or be a part of a more complex longer light curve. The radiative efficiency was high at ∼23% due to the reverse shock existing at low optical depths. The hardness peaked at the light curve peak, and then slowly decayed due to adiabatic cooling in the comoving frame and high-latitude emission.

2. The time-integrated spectrum consisted of a smooth curvature below the peak energy a high-energy power law that extended to ∼5 MeV. For a redshift of z = 1 and S/N  = 50, it had best-fit Band parameters of α = −0.96 ± 0.06, β = −2.23 ± 0.09, and E_p = 156 ± 16 keV. It was well fitted with a Band function in the low S/N region, but was much better fitted with a 2SBPL function in the high S/N region.

3. A high-energy power law extending for 1.5 orders of magnitude in energy was generated due to the injection of high-energy photons from the reverse shock. The maximum energy in the reverse shock increased at lower optical depths partly because the shock sped up in the decreasing density profile and partly because the decreasing scattering time led to a steeper velocity gradient across the shock (in units of optical depth).

4. The peak energy remained ∼ constant during the pulse rise due to the existence of the reverse shock. After the reverse shock crossing, adiabatic cooling led to the observed peak energy decreasing. The high-energy index showed the opposite behaviour with a rapid decrease during the pulse rise, after which it stayed ∼ constant. The decrease was due to a rapid downscattering of high-energy photons, a process that becomes inefficient at low optical depths (τ ≲ 5).

5. The very early spectrum probed the deeper regions of the jet and showed a more complex behaviour, where both reverse and forward shocks were visible.

In conclusion, we find that RMSs show several promising features with regards to the prompt GRB emission. However, we have only studied a single example in this paper. It remains to be tested if the model can reproduce the diversity of observed light curves and pulse durations, and under which conditions some shocks can still propagate above the photosphere with a different contribution to the emission.

Acknowledgments

We thank the anonymous referee for their insightful comments. This research has made use of the High Energy Astrophysics Science Archive Research Center (HEASARC) Online Service at the NASA/Goddard Space Flight Center (GSFC). We thank the GBM team for providing the necessary tools and data. F.A. is supported by the Swedish Research Council (Vetenskapsrådet; 2022-00347). F.D. acknowledges the Centre National d’Études Spatiales (CNES) for financial support in this research project.

References

Appendix A: Shock broadening and the dynamical stage of an RMS

The width of the RMS is set by the number of scatterings the photons have to do to decelerate the incoming flow. As the photon mean free path increases with time due to the decreasing lepton density in the jet, the physical width of the RMS grows.

The instantaneous optical depth across the shock transition region, Δτ_r, (not to be confused with the optical depth towards the observer, τ) is given by an integral over the shock transition region as

$$\begin{array}{c}\hfill \mathrm{\Delta }{\tau }_r=\int n_e{\sigma }_{\mathrm{T}}d{R}^{\prime }\approx n_e{\sigma }_{\mathrm{T}}\mathrm{\Delta }{R}_r^{\prime },\end{array}$$(A.1)

where n_e is the comoving electron number density,⁸σ_T is the Thomson cross section, ΔR′_r is the width of the shock as measured in the downstream comoving frame, and we assumed ΔR′_r ≪ R and a non-relativistic RMS. The number of photons in the shock transition region is given by N_r = 4πR²n_γΔR′_r, where n_γ is the comoving photon number density. Thus, the photon number in the RMS is given by

$$\begin{array}{c}\hfill N_r(R)=\frac{4\pi R^2\zeta \mathrm{\Delta }{\tau }_r}{{\sigma }_{\mathrm{T}}},\end{array}$$(A.2)

where the photon-to-proton ratio, ζ, is equal to the photon-to-electron ratio in the absence of pairs.⁹ The photon-to-proton ratio remains ∼ constant across the photon transition in a photon rich shock. Equation (A.2) implies that the photon number in the RMS increases as N_r ∝ R², with instantaneous change in the lab frame as

$$\begin{array}{c}\hfill \frac{d{N}_r}{\mathit{dR}}=\frac{8\pi R\zeta \mathrm{\Delta }{\tau }_r}{{\sigma }_{\mathrm{T}}}=\frac{2}{R}{N}_r,\end{array}$$(A.3)

assuming that ζ and Δτ_r vary slowly with R.

Lets assume that the probability for a photon to leave the shock per scattering time is inversely proportional to the average number of scatterings in the shock, ${\mathcal{N}}_{\mathrm{sc}}$. Then the outgoing photon number is given by

$$\begin{array}{c}\hfill \frac{d{N}_{\mathrm{out}}}{\mathit{dR}}=\frac{N_r}{{\mathcal{N}}_{\mathrm{sc}}}{({\beta }_{r}c{t}_{\mathrm{sc}})}^{-1}\approx \frac{N_r}{{\mathcal{N}}_{\mathrm{sc}}}\frac{\tau }{R},\end{array}$$(A.4)

where β_r ≈ 1 is the dimensionless shock speed in the lab frame, and we used t_sc = Γ/ρκc and τ ≈ ρκR/Γ for the last approximation. Combining Equations (A.3) and (A.4), we get the ingoing photon number as

$$\begin{array}{c}\hfill \frac{d{N}_{\mathrm{in}}}{\mathit{dR}}=\frac{d{N}_r}{\mathit{dR}}+\frac{d{N}_{\mathrm{out}}}{\mathit{dR}}\approx \frac{N_r}{R}(2+\frac{\tau }{{\mathcal{N}}_{\mathrm{sc}}}).\end{array}$$(A.5)

From the equation above, it is clear that in the optically thick regime τ > 1, there exists two different limiting cases in the shock evolution. When ${\mathcal{N}}_{\mathrm{sc}}\ll \tau /2$, the outgoing photon number roughly equals the incoming photon number. In a steady outflow with ρ ∝ R⁻², the optical depth is halved during one dynamical time. Thus, when ${\mathcal{N}}_{\mathrm{sc}}\ll \tau /2$, photons cross the RMS well within a dynamical time and the shock has time to reach a steady state in response to the slowly varying upstream properties. However, once ${\mathcal{N}}_{\mathrm{sc}}\sim \tau /2$, the dynamical stage of the shock evolution starts. At this stage, the shock crossing timescale equals the dynamical time and the shock will struggle to respond to the changing upstream properties. Due to adiabatic cooling in the upstream, the average number of scatterings in the shock is not enough to dissipate the necessary energy and decelerate the incoming flow. The velocity gradient (in units of photon mean free paths) across the shock must steepen, which in turn will increase the average energy gain per scattering for the photons until the shock jump conditions are satisfied. If the shock jump conditions are not satisfied even with the steeper velocity gradient, a collisionless subshock must form (Ito et al. 2020; Lundman & Beloborodov 2021). Note that at this stage, Δτ_r is no longer slowly varying with R and Equation (A.5) breaks down.

For non-relativistic RMSs, the number of scatterings follows ${\mathcal{N}}_{\mathrm{sc}}\propto {\beta }_u^{-2}$ (Blandford & Payne 1981). This implies that the transition to a dynamically evolving shock can happen deep below the photospheric radius, much earlier than the shock breakout time.

The same result can be obtained using Equation (A.1), the approximate scaling Δτ = β_u⁻¹ (e.g. Levinson & Nakar 2020), and n_e = Γτ/σ_TR valid for a steady wind in the absence of pairs. Then, one finds that

$$\begin{array}{c}\hfill t_{\mathrm{cross}}=\frac{t_{\mathrm{dyn}}}{{\beta }_u^2\tau },\end{array}$$(A.6)

where t_cross = ΔR′_r/β_uc is the shock crossing time and t_dyn = R/Γc is the dynamical time. Thus, the condition t_cross ≪ t_dyn requires τ ≫ 1/β_u². Shock breakout occurs once the optical depth towards the observer is of the order Δτ_r, which means that at breakout, the crossing time is a factor β_u⁻¹ larger than the dynamical time according to the equation above.

Appendix B: Obtaining the RMS parameters from the hydro output

B.1. Upstream, immediate downstream, and far downstream properties

The hydrodynamical simulation produced two shocks, which separated the ejecta into five regions (see Section 2.4 and Figure 2). GAMMA comes with its own shock finding algorithm. However, here we instead used Equation (72) in Martí & Müller (1996) (see also Colella & Woodward 1984), which considers cell j to be within a shock if the relative pressure increase across the cell satisfies

$$\begin{array}{c}\hfill \frac{|p_{j+1}-p_{j-1}|}{min(p_{j+1},p_{j-1})}>\delta ,\end{array}$$(B.1)

where δ is a constant that we set to 1 (in accordance with Martí & Müller 1996). Equation (72) in Martí & Müller (1996) also includes a velocity condition that has been neglected here, since we wanted to identify both left-moving and right-moving shocks.

Since the boundary condition could lead to spurious shock detections close to the boundary, Equation (B.1) was only evaluated between Lagrangian mass coordinates m_ll = 0.01M_tot and m_ul = 0.99M_tot.

Numerically, the discontinuity corresponding to a shock is always smoothed over a few cells. Therefore, Equation (B.1) identified several consecutive cells as being within the shock, for both the forward and the reverse shock (of the order of three-four cells for the used grid resolution). The upstream properties in each saved snapshot were taken to be the properties in the cell that was five cells upstream of the shocked cells. Similarly, the immediate downstream properties were taken at five cells downstream of the shocked cells. As an example, say that Equation (B.1) identified cells 293, 294, and 295 as within the reverse shock and cells 741, 742, 743 and 744 as within the forward shock. The upstream (immediate downstream) properties for the reverse shock and the forward shock were then taken to be the properties of cell 288 (300) and cell 749 (736), respectively. The choice of five was found to give a good representation of the pressure, density, and velocity in the upstream and the immediate downstream for the current grid resolution.

For the properties in the downstream region, we set the value of a quantity Q at time t to equal the mass-averaged value between the two shocks as

$$\begin{array}{c}\hfill Q(t)=\frac{1}{M_{\mathrm{tot}}}\int Q(t,m)dm,\end{array}$$(B.2)

where the lower (upper) limit to the integral is the mass coordinate of the reverse (forward) shock at time t. The bulk Lorentz factor, the scattering time, and the adiabatic cooling term for the downstream region were all calculated this way.

B.2. Obtaining the RMS parameters

Using the shock jump conditions from Beloborodov (2017) (see also Samuelsson et al. 2022), we solved for the upstream Lorentz factor of the incoming material in the shock rest frame as

$$\begin{array}{c}\hfill {\gamma }_u^2=\frac{1-{({\rho }_u/{\rho }_d)}^2}{{(1+w_u)}^2/{(1+w_d)}^2-{({\rho }_u/{\rho }_d)}^2},\end{array}$$(B.3)

where w = 4p/ρc² is the dimensionless enthalpy obtained directly from the hydro.

The photon-to-proton ratio was estimated by assuming a thermodynamic equilibrium at R₀, a negligible energy in magnetic fields compared to the internal energy of radiation and pairs at R₀, and that all pairs had recombined at the collision radius. This gives (Bromberg et al. 2011; Levinson & Nakar 2020; Alamaa 2024)

$$\begin{array}{c}\hfill \zeta =\frac{1}{4{\theta }_0}\frac{\eta }{{\mathrm{\Gamma }}_0}\frac{m_p}{m_e},\end{array}$$(B.4)

where θ₀ is the temperature at R₀ in units electron rest mass given by

$$\begin{array}{c}\hfill {\theta }_0=\frac{k_{\mathrm{B}}}{m_e{c}^2}{\left[\frac{3}{11}\frac{\dot{E}}{4\pi R_0^{2}ac{\mathrm{\Gamma }}_0^2{\beta }_0}\right]}^{1/4}.\end{array}$$(B.5)

In the equation above, a is the radiation constant and ${\beta }_0=\sqrt{1-{\mathrm{\Gamma }}_0^{-2}}$. Since ζ depends on η and Ė, it varied across the ejecta. However, as given in Equation (B.4), it remains constant in time for a given mass coordinate. Therefore, in our Lagrangian framework, each cell had one constant value of ζ associated with it. This is an approximation, as photons would flow between different mass elements in reality. However, in the high optical depth regions treated here, photons are tied to the plasma and this approximation was likely quite good.

Since the photons completely dominated the particle number, the upstream temperature in units electron rest mass energy was given by θ_u = p_u/n_u, γm_ec². Thus,

$$\begin{array}{c}\hfill {\theta }_u=\frac{w_u}{4}\frac{m_p}{m_e}\frac{1}{{\zeta }_u}.\end{array}$$(B.6)

Finally, the average photon energy in the immediate downstream in units electron rest mass energy was obtained from the relation ${\overline{ϵ}}_d=3p_d/n_{d,\gamma }{m}_e{c}^2$ where we used a relativistic equation of state and assumed that the photons dominate the downstream pressure. This gives

$$\begin{array}{c}\hfill {\overline{ϵ}}_d=\frac{3w_d}{4}\frac{m_p}{m_e}\frac{1}{{\zeta }_d}.\end{array}$$(B.7)

The downstream average energy, ${\overline{ϵ}}_d$, can be obtained from γ_u, ζ_u, and θ_u using the shock jump conditions, i.e. it is not a free parameters if the other three are known.

Appendix C: Details regarding the KRA in a dynamical flow

In this appendix, the interested reader can find all necessary details regarding the dynamical version of the KRA.

C.1. Geometry of the KRA with a reverse and a forward shock

The original version of the KRA only treated one RMS. However, the dynamical version of the KRA can model both the forward and the reverse shock resulting from a collision. It is possible to run the dynamical version of the KRA with only one RMS as well. In other words, all the equations written in this appendix are applicable in the case of either one or two RMSs, unless otherwise stated. Furthermore, the reverse and forward shocks are modelled identically. Thus, we do not generally state which shock a specific quantity, such as θ_u, K or θ_r, belongs to.

The material flowing through the forward and reverse shocks is situated in a common downstream with ∼ constant pressure and velocity. The observed width of the common downstream marginally satisfies ΔR ≲ R/Γ² (Levinson 2012), which means that the downstream light crossing time is shorter than the dynamical time. Therefore, we choose to treat the photons in the downstream as a single fluid, regardless of which shock they originally crossed. For the simulation presented in Section 5, the condition ΔR ≲ R/Γ² broke down when τ < 2. However, at these small optical depths, the photon distribution changes mostly due to adiabatic cooling and the single zone approximation makes no difference.

To evolve both the forward and the reverse shock, the KRA has five different zones: 1) the reverse shock upstream zone, 2) the reverse shock zone, 3) the common downstream zone, 4) the forward shock zone, and 5) the forward shock upstream zone. This geometry is shown in Figure 2.

In each zone, photons interact with the plasma electrons via thermal Comptonisation. The photons in the upstream zones are assumed to be in a thermal Wien distribution at the corresponding upstream temperature, N_u, ϵ ∝ ϵ²exp(−ϵ/θ_u, K). The effect of the dissipation of kinetic energy in a true RMS is modelled by enforcing the electrons in the shock zones to have a high effective temperature. The downstream temperature is assumed to equal the Compton temperature, which is analytically obtained from the downstream photon distribution, N_*, ϵ, as (Beloborodov & Illarionov 1995)

$$\begin{array}{c}\hfill {\theta }_{\ast ,\mathrm{C}}=\frac{1}{4}\frac{\int {ϵ}^2{N}_{\ast ,ϵ}dϵ}{\int ϵN_{\ast ,ϵ}dϵ}.\end{array}$$(C.1)

Since θ_*, C is a function of the downstream photon distribution, it needs to be continuously updated as the distribution evolves.

The different zones are coupled via source terms as indicated in Figure 2. These are given in subsection C.3.

C.2. New version of the Kompaneets equation

In Samuelsson et al. (2022), the Kompaneets equation was written as a function of the normalised radius $\overline{R}=R/R_{\mathrm{ph}}$. This is not a suitable choice when the downstream properties, such as the bulk Lorentz factor, can vary, since the estimated value of the photospheric radius will change over time. Thus, we rewrite the Kompaneets equation as a function of R. Using the generic subscript z to refer to an arbitrary KRA zone, and neglecting induced scatterings and photon creation and absorption, it can be written as (Vurm & Beloborodov 2016)

$$\begin{array}{c}\hfill \begin{array}{cc}& \frac{\partial {\mathcal{F}}_z}{\partial R}=\frac{s_z(ϵ,R)}{{ϵ}^2}+\hfill \\ \hfill & \frac{1}{{ϵ}^2}\frac{\partial }{\partial ϵ}[\frac{{ϵ}^4}{t_{\mathrm{sc},z}c}({\theta }_z\frac{\partial {\mathcal{F}}_z}{\partial ϵ}+{\mathcal{F}}_z)+{\phi }_z\frac{{ϵ}^3{\mathcal{F}}_z}{R}].\hfill \end{array}\end{array}$$(C.2)

In the equation above, s_z is a source term (further discussed in the next subsection) and ${\mathcal{F}}_z$ is related to the photon spectral number density in zone z as ${\mathcal{F}}_z\equiv N_{z,ϵ}/{ϵ}^2$. The total number of photons in the zone is given by

$$\begin{array}{c}\hfill N_z={\int }_0^{\infty }{N}_{z,ϵ}dϵ.\end{array}$$(C.3)

The Kompaneets equation is solved using the numerical scheme outlined in Chang & Cooper (1970).

The second term in the brackets in Equation (C.2) accounts for photons cooling adiabatically as the ejecta expands (see Equation (8)). It can be shown using Equation (C.2) that when sources are absent, the average photon energy in zone z, ${\overline{ϵ}}_z$, decreases due to adiabatic cooling as $dlog{\overline{ϵ}}_z/dlogR=-{\phi }_z$. However, the cooling of the photon field should cease when the radiation decouples from the plasma (Pe’er 2008; Beloborodov 2011). To account for this, the KRA shuts off adiabatic cooling above a specific radius, R_cool, evaluated such that the total adiabatic cooling in the KRA equals the total cooling the photon field would experience in reality. The value of R_cool depends weakly on φ, but is roughly equal to R_cool ≳ R_ph/1.2.

C.3. Source terms in Komrad dynamical

In Komrad, the number of photons in the RMS zone was assumed to be constant, i.e. dN_out/dR = dN_in/dR. From the discussion in Appendix A, this implies an implicit assumption of ${\mathcal{N}}_{\mathrm{sc}}\ll \tau /2$. For Komrad dynamical, this assumption is relaxed.

The outgoing number of photons is assumed to be given by Equation (A.4) (the expression before the last approximation). Since the probability to be advected from the shock zone into the downstream zone is assumed to be independent of the photon energy, the energy distribution of the escaping photons is the same as that of the photon distribution in the shock zone. The outgoing source term is thus

$$\begin{array}{c}\hfill s_{\mathrm{out}}=\frac{d{N}_{\mathrm{out}}}{\mathit{dR}}\frac{N_{r,ϵ}}{N_r}=\frac{N_{r,ϵ}}{c{t}_{\mathrm{sc},r}{\mathcal{N}}_{\mathrm{sc}}}.\end{array}$$(C.4)

To get a more accurate description of the incoming photon number, we do not assume that ζ_r and Δτ_r vary slowly with R as in Equation (A.3). The more general expression for the instantaneous change of photon number in the RMS is

$$\begin{array}{c}\hfill \frac{d{N}_r}{\mathit{dR}}=\frac{2}{R}{N}_r+\frac{4\pi R^2}{{\sigma }_{\mathrm{T}}}\frac{d(\zeta \mathrm{\Delta }{\tau }_r)}{\mathit{dR}}.\end{array}$$(C.5)

In this paper, the evolution of ζΔτ_r with R is deduced from the hydro simulation, assuming ζ = ζ_u.

To assure the correct photon number in the RMS zone, the incoming photon number is given by

$$\begin{array}{c}\hfill \frac{d{N}_{\mathrm{in}}}{\mathit{dR}}=\frac{d{N}_r}{\mathit{dR}}+\frac{d{N}_{\mathrm{out}}}{\mathit{dR}},\end{array}$$(C.6)

where dN_r/dR is given by Equation (C.5) and $d{N}_{\mathrm{out}}/dR=N_r/{\mathcal{N}}_{\mathrm{sc}}/c{t}_{\mathrm{sc},r}$.

The incoming photon number as given in Equation (C.6) is not necessarily positive. If the term is negative, this corresponds to the RMS width shrinking, and photons being “advected” into the upstream from the shock. Komrad dynamical allows for this by using the source term

$$\begin{array}{c}\hfill s_{\mathrm{in}}=\frac{d{N}_{\mathrm{in}}}{\mathit{dR}}\frac{N_{(u/r),ϵ}}{N_{(u/r)}},\end{array}$$(C.7)

where the subscript (u/r) implies upstream zone quantities if dN_in/dR > 0 and RMS zone quantities if dN_in/dR < 0. The source terms used in Equation (C.2) are then

$$\begin{array}{c}\hfill s_r=s_{\mathrm{in}}-s_{\mathrm{out}}\end{array}$$(C.8)

for the RMS zone(s) and

$$\begin{array}{c}\hfill s_{\ast }=s_{\mathrm{out}}^{\mathrm{rs}}+s_{\mathrm{out}}^{\mathrm{fs}}\end{array}$$(C.9)

for the common downstream zone. In the case of a single shock, the source term for the downstream will only contain one term.

That photons are allowed to flow from the RMS to the upstream means that the assumption of a thermal upstream photon distribution is not completely valid. However, in practice, this happens very rarely in the simulation and in this paper, we ignore this effect.

To avoid injecting too many photons, the source term linking the upstream zone and the RMS zone is set to zero once the upstream has been depleted. The upstream is said to have been depleted once the condition

$$\begin{array}{c}\hfill \int \frac{d{N}_{\mathrm{in}}}{\mathit{dR}}dR\ge N_{u,\mathrm{tot}}\end{array}$$(C.10)

is satisfied. After this time, the RMS zone still advects photons into the downstream zone according to Equation (C.4). It also continues to evolve according to the Kompaneets equation.

C.4. Energy conservation across the RMS zone

The shock jump conditions dictate what the average energy of the photon distribution in the immediate downstream, ${\overline{ϵ}}_d$, must be depending on the upstream conditions (see Equation (B.7)). To ensure energy conservation in the KRA model, the photon distribution that is advected into the downstream zone should have an average energy that equals ${\overline{ϵ}}_d$. According to Equation (C.4), this implies that ${\overline{ϵ}}_r={\overline{ϵ}}_d$ at all times, where ${\overline{ϵ}}_r$ is the average energy in the RMS zone. In principle, this is achieved by choosing the appropriate value for y_r as explained in Section 3.1.

The appropriate value of y_r when the RMS is in steady state can be obtained from the empirical relation given in Bagi et al. (2025) (see discussion in Section D.1). However, at early times during the shock formation and at late times during the dynamical stage of the shock, the radiation in the RMS will not be in steady state. To ensure energy conservation even at these times, the value of y_r must be continuously adapted.

The Compton y-parameter for the RMS zone is defined as $y_r=4{\theta }_r{\mathcal{N}}_{\mathrm{sc}}$. The value of y_r can thus be changed by changing either ${\mathcal{N}}_{\mathrm{sc}}$ or θ_r. In this paper, we assume that the change in y_r necessary to satisfy the shock jump conditions is entirely due to a corresponding change in θ_r. The real picture is likely more complicated with a correlated change to both quantities. Indeed, an increase in the energy gain per scattering across an RMS implies a steeper velocity gradient (in units of photon mean free paths), which in turn must effect the average number of scatterings. Such nonlinear effects require detailed numerical simulations, which are outside the scope of the current paper.

The RMS zone is initiated with a thermal Wien distribution at temperature ${\theta }_r={\overline{ϵ}}_d/3$ (see Appendix A in Samuelsson & Ryde 2023, for motivation). Thus, the average photon energy is initially what it should be according to the shock jump conditions. Therefore, it is enough to ensure that $\partial {\overline{ϵ}}_r/\partial R=\partial {\overline{ϵ}}_d/\partial R$ to ensure ${\overline{ϵ}}_r={\overline{ϵ}}_d$ for all R. Using Equation (C.2), together with the source terms in Equations (C.4) and (C.7), it can be shown that the average photon energy in the RMS zone varies with radius as

$$\begin{array}{c}\hfill \frac{\partial {\overline{ϵ}}_r}{\partial R}=\frac{4{\overline{ϵ}}_r}{c{t}_{r,\mathrm{sc}}}({\theta }_r-{\theta }_{r,\mathrm{C}})-\frac{{\phi }_r{\overline{ϵ}}_r}{R}+\frac{1}{N_r}\frac{d{N}_{\mathrm{in}}}{\mathit{dR}}({\overline{ϵ}}_{(u/r)}-{\overline{ϵ}}_r),\end{array}$$(C.11)

where ${\overline{ϵ}}_{(u/r)}$ equals ${\overline{ϵ}}_u$ if dN_in/dR > 0 and ${\overline{ϵ}}_r$ if dN_in/dR < 0, and θ_r, C is the Compton temperature in the RMS zone. Solving the above equation for θ_r and inserting $\partial {\overline{ϵ}}_r/\partial R=\partial {\overline{ϵ}}_d/\partial R$ gives the desired behaviour. We have ensured that the prescription works by comparing the input ${\overline{ϵ}}_d(R)$ with ${\overline{ϵ}}_r(R)$ as obtained in Komrad dynamical and the agreement is excellent.

For most of the simulation, θ_r is obtained from Equation (C.11) as explained above. However, at the final stages when the upstream has been depleted and the RMS zone is being drained of its remaining photons, θ_r is set equal to the Compton temperature in the RMS zone.

C.5. Parameters for the dynamical KRA

Running a KRA simulation as described in this section requires input of the following parameters:

For the upstream zone(s),

• θ_u,K(R), the temperature.

• N_u,tot, the total photon number.

For the RMS zone(s),

• 𝒩_sc(R), the steady state number of scatterings.

• φ_r(R), the adiabatic cooling term.

• t_rsc(R), the observed scattering time.

• ${\overline{ϵ}}_d(R)$, the average photon energy in the immediate downstream.

For the downstream zone,

• φ_r(R), the adiabatic cooling term.

• t_∗,sc(R), the observed scattering time.

In addition to the parameters above, one also needs to know $\frac{d(\zeta \mathrm{\Delta }{\tau }_r)}{\mathit{dR}}(R)$. However, this value can be obtained from the other parameters and is thus not a free parameter.

The above list results in eight parameters in the single shock case and fourteen parameters in the reverse plus forward shock case. Furthermore, most of these are functions of R, meaning that we need to now how they evolve with radius. However, it is important to note that these parameters may be coupled in non-trivial ways. Indeed, the hydrodynamical simulation in this paper only had eight free parameters but fixed all sixteen input parameters as a function of radius for the KRA.

Compared to its previous version, it is clear that the dynamical version of the KRA is much more complex. However, the lower complexity of the previous version can be recovered with the following simplifying assumption. Under the assumption of ρ ∝ R⁻², one gets φ_r = φ_* = 2/3. Assuming that the shock(s) dissipates all of its energy during a doubling of the radius, N_u, tot is not a free parameter. Lastly, it was previously assumed that t_r, sc(R) = t_*, sc(R) and that Γ was constant. Under these assumptions, the KRA with all its new modifications can be run, and the subphotospheric evolution of the photon distribution obtained, knowing only the initial values of θ_u, K, ${\mathcal{N}}_{\mathrm{sc}}$, ${\overline{ϵ}}_d$, and τ, which is the same amount of free parameters as before.

Appendix D: Obtaining the KRA input from the hydro output

D.1. Converting between the RMS and the KRA parameters

As explained in Section 3.1, a steady state RMS zone spectrum in the KRA is characterised by three parameters: the comoving upstream temperature θ_u, K, the effective temperature in the RMS zone, ${\stackrel{\sim }{\theta }}_r$, and the Compton y-parameter of the RMS zone, ${\stackrel{\sim }{y}}_r$. Here, the tilde indicates steady state. The number of scatterings in steady state is then obtained as ${\mathcal{N}}_{\mathrm{sc}}={\stackrel{\sim }{y}}_r/4{\stackrel{\sim }{\theta }}_r$.

The upstream temperature θ_u, K is larger than θ_u by a factor of (ρ_d/ρ_u)^1/3 to account for the adiabatic compression that occurs across a real RMS (Blandford & Payne 1981). Thus, it is given by

$$\begin{array}{c}\hfill {\theta }_{u,\mathrm{K}}={\left(\frac{{\rho }_d}{{\rho }_u}\right)}^{1/3}{\theta }_u.\end{array}$$(D.1)

The compression factor above is ∼ two for the relevant parameter space.

In Samuelsson et al. (2022), the effective temperature ${\stackrel{\sim }{\theta }}_r$ was found by equating ${ϵ}_{\max }^{\mathrm{RMS}}={ϵ}_{\max }^{\mathrm{KRA}}$. Via analytical arguments and an empirical study, it was determined that

$$\begin{array}{c}\hfill 4{\stackrel{\sim }{\theta }}_r=\frac{{({\gamma }_u{\beta }_u)}^{2}ln({\overline{ϵ}}_d/{\overline{ϵ}}_u)}{\xi },\end{array}$$(D.2)

and Samuelsson et al. (2022) found that ξ ≈ 55 gives good agreement across a wide range of initial parameters.

The parameter ${\stackrel{\sim }{y}}_r$ is a measure of the energy increase of a typical photon within the RMS zone. It is obtained by requiring that the average energy of the photons that are advected into the downstream zone in the KRA equals the average energy immediately downstream of the shock in the hydrodynamical simulation, i.e. ${\overline{ϵ}}_r={\overline{ϵ}}_d$. In the KRA model, the value of ${\stackrel{\sim }{y}}_r$ corresponding to a specific average photon energy ${\overline{ϵ}}_r$ is not known a priori. In principle, a simulation of the Kompaneets equation can be run in steady state iteratively until the unique value of ${\stackrel{\sim }{y}}_r$ that corresponds to ${\overline{ϵ}}_r$ is found. To save computational time, here, we employed the empirical relation between ${\stackrel{\sim }{y}}_r$ and ${\overline{ϵ}}_r$ that is given in Bagi et al. (2025). The relation is accurate to within 5% in the relevant parameter region.

D.2. Scattering time, adiabatic cooling term, and evolution of ζΔτ_r

As mentioned in Section B.1, the values for the common downstream zone density and bulk Lorentz factor were taken as the mass averaged value of each quantity between the two shock fronts. With ρ_* and Γ_* obtained, φ_* and t_*, sc were calculated using Equations (8) and (9), respectively.

The scattering time in the RMS zone(s) was slightly more complicated, since the plasma density and bulk Lorentz factor changes significantly across the RMS transition region. In this paper, we used

$$\begin{array}{c}\hfill t_{r,\mathrm{sc}}^{-1}=\frac{t_{d,\mathrm{sc}}^{-1}+t_{u,\mathrm{sc}}^{-1}}{2},\end{array}$$(D.3)

that is, we set the RMS zone scattering frequency to be the average of the upstream and immediate downstream scattering frequencies. In the adiabatic cooling term, we used the immediate downstream value for the density.

To solve Equation (C.5), we needed d(ζΔτ_r)/dR. For the photon-to-proton ratio, we used the value from Equation (B.4) evaluated in the upstream region. The optical depth across the shock was evaluated in the non-relativistic limit as $\mathrm{\Delta }{\tau }_r=\sqrt{{\mathcal{N}}_{\mathrm{sc}}}=\sqrt{{\stackrel{\sim }{y}}_r/4{\stackrel{\sim }{\theta }}_r}$.

D.3. Initial photon distributions

In this paper, the downstream zone was not empty at the beginning of the simulation. This is because the initial plasma is rather hot, and the reverse and the forward shocks did not form at the same mass coordinate. Let us denote the mass coordinate where the reverse shock formed by m_rs, and the mass coordinate where the forward shock formed by m_fs, where m_fs > m_rs. The number of photons in the downstream at the start of the simulation was then given by

$$\begin{array}{c}\hfill N_{\ast ,\mathrm{init}}={\int }_{m_{\mathrm{rs}}}^{m_{\mathrm{fs}}}\zeta (m)\frac{\mathit{dm}}{m_p}.\end{array}$$(D.4)

These photons were put in a thermal Wien distribution at temperature θ_* = w_cdm_p/4m_eζ_cd (see Equation (B.6)). Here, w_cd and ζ_cd are the dimensionless enthalpy and the photon-to-proton ratio at the contact discontinuity, respectively.

As mentioned in Section C.4, the RMS zone(s) were initiated with a thermal Wien distribution at a temperature of ${\overline{ϵ}}_d/3$. The number of photons in this initial distribution was given by (see Equation (A.2))

$$\begin{array}{c}\hfill N_{r,\mathrm{init}}=\frac{4\pi R^2\zeta \mathrm{\Delta }{\tau }_r}{{\sigma }_{\mathrm{T}}},\end{array}$$(D.5)

where the values of R, ζ, and Δτ_r were taken to be the respective values at the radius where the shock forms. In the hydro simulation, the forward shock and the reverse shock did not form at the same time. The KRA simulation was set to start once the first of the two shocks had formed. The second RMS zone was included once the KRA simulation had reached the shock formation radius of the second RMS.

The total number of photons that would pass through each shock was obtained by integrating the photon number in between the mass coordinate where the shock started (m_rs or m_fs) and the mass coordinate where the shock would end (m_ll or m_ul). The initial total photon number in the upstream was then given by

$$\begin{array}{c}\hfill N_{u,\mathrm{tot}}={\int }_{m_1}^{m_2}\zeta (m)\frac{\mathit{dm}}{m_p}-N_{r,\mathrm{init}},\end{array}$$(D.6)

where the mass coordinate limits m₁ and m₂ equal m_ll (m_fs) and m_rs (m_ul) for the reverse shock upstream (forward shock upstream), respectively.

All Tables

All Figures

Fig. 1. Flowchart showing the methodology of the paper. Each box shows one step of the chain, with the treated physics given in light blue and the simulation code used given in red. The final result is the time-resolved signal in the observer frame. The hydrodynamics are explained in Section 2, the RMS modelling in Section 3, and the photon decoupling in Section 4. The observed time-resolved signal is shown in Section 5.

In the text

Fig. 2. Top panel: Sketch of the Lorentz factor profile across the ejecta with the positions of the five different regions mentioned in subsection 2.4 marked. The dashed vertical lines indicate the position of each region relative to the KRA zones given in the bottom panel. Bottom panel: Geometry of the KRA with both reverse and forward shocks included. Green, red, and purple colours indicate the upstream zones, the RMS zones, and the common downstream zone, respectively. The zones are coupled via source terms as indicated in the figure. The RMS and downstream zones are evolved using the Kompaneets equation (Equation (C.2)) while the upstream zones are assumed to be in a thermal Wien distribution. The relevant temperature in each zone is given. This panel is adapted from Samuelsson et al. (2022), which shows the geometry of the KRA with three zones. In between the two panels, the subscripts used to refer to the different regions and zones are shown for clarity.

In the text

Fig. 3. Lorentz factor (top) and the comoving density (bottom) at different optical depths as obtained from the hydrodynamical simulation using the vacuum boundary conditions. The dashed black line gives the corresponding profile at the start of the simulation. The comoving photon distributions in the five different zones of Komrad dynamical at the corresponding times are given in Figure 4. The simulation parameters are given in Table 1.

In the text

Fig. 4. Comoving photon distribution in the five different zones as obtained by Komrad dynamical. The line coding is shown by the legend in the top-left panel, where RS and FS stand for reverse shock and forward shock, respectively. The four different panels show the distributions at different optical depths, corresponding to the optical depths in Figure 3. The observed photon distribution consists of a superposition of comoving spectra, leading to a smoothening and broadening effect on the observed spectrum. The simulation parameters are given in Table 1.

In the text

Fig. 5. Observed spectral evolution (top left), light curve (top right), count spectrum of a Band fit to the time-integrated spectrum (bottom left), and the evolution of the best-fit Band parameters (bottom right). The dashed grey line in the top-left panel shows the time-integrated spectrum taken as the average over the first second of observations, and purple shading shows the energy range of the GBM, with the darker colour indicating the region with the highest effective area. The dot-dashed red line in the light curve shows the light curve in the BGO band scaled to the peak of the NaI light curve to highlight the difference between the two. The red, blue, cyan, and green data points in the count spectrum show the three NaI detectors and the one BGO detector used for the fit, respectively. The best-fit peak energy is given in units Ep, 2 = Ep/100 keV, and the mock dataset for the fit had a S/N = 50. The best-fit parameters with errors for the time-integrated spectrum were α = −0.96 ± 0.06, β = −2.23 ± 0.09, and Ep = 156 ± 16 keV. All errors given represent 1σ. The simulation parameters are given in Table 1.

In the text

Fig. 6. Comparison of a Band function fit (top) and a 2SBPL fit (bottom) to the time-integrated spectrum with S/N = 300. The red, blue, cyan, and green data points in the count spectrum show the three NaI detectors and the one BGO detector used for the fit, respectively. It is clear from the uneven residuals in the top panel that the Band function did not produce a satisfactory fit. The 2SBPL function performed much better in this region of high S/N, being statistically preferred with ΔAIC = 114. The best-fit Band parameters were α = −0.92 ± 0.02, β = −2.33 ± 0.02, and Ep = 152 ± 4 keV, and the best-fit 2SBPL were α1 = −0.69 ± 0.06, α2 = −1.45 ± 0.04, β = −2.49 ± 0.04, Ep = 178 ± 7 keV, Eb = 27 ± 3 keV.

In the text