# Numerical simulations of internal shocks in spherical geometry: hydrodynamics and prompt emission

A Charlet, J Granot, P Beniamini

► To cite this version:

A Charlet, J Granot, P Beniamini. Numerical simulations of internal shocks in spherical geometry: hydrodynamics and prompt emission. Astronomy & Astrophysics - A&A, 2025, 699, pp.A139. 10.1051/0004-6361/202453351. [hal-04891889](https://hal.science/hal-04891889v1)

HAL Id: hal-04891889

https://hal.science/hal-04891889v1

Submitted on 16 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Numerical simulations of internal shocks in spherical geometry:
## Hydrodynamics and prompt emission

A. Charlet<fnref n="1" />,<fnref n="2" />,*, J. Granot<fnref n="1" />,<fnref n="3" />,<fnref n="4" />, and P. Beniamini<fnref n="1" />,<fnref n="3" />,<fnref n="4" />

Received 9 December 2024 / Accepted 20 May 2025

ABSTRACT

*Context.* Among the models used to explain the prompt emission of gamma-ray bursts (GRBs), internal shocks is a leading one. Its most basic ingredient is a collision between two cold shells of different Lorentz factors in an ultrarelativistic outflow, which forms a pair of shock fronts that accelerate electrons in their wake. In this model, key features of GRB prompt emission such as the doubly broken power-law spectral shape arise naturally from the optically thin synchrotron emission at both shock fronts.

**Aims.** We investigate the internal shocks model as a mechanism for prompt emission based on a full hydrodynamical analytic derivation in planar geometry, extending this approach to spherical geometry using hydrodynamic simulations.

*Methods.* We used the moving mesh relativistic hydrodynamics code GAMMA to study the collision of two ultrarelativistic cold shells of equal kinetic energy (and power). Using the built-in shock detection, we calculated the corresponding synchrotron emission by the relativistic electrons accelerated into a power-law energy distribution behind the shock in the fast-cooling regime.

*Results.* During the first dynamical time after the collision, the spherical effects cause the shock strength to decrease with radius. The observed peak frequency decreases faster than expected by other models in the rising part of the pulse and the peak flux is saturated even for moderately short pulses. This is likely caused by the very sharp edges of the shells in our model, while smoother edges would probably mitigate this effect. Our model traces the evolution of the peak frequency back to the source activity time scales.

**Key words.** hydrodynamics – radiation mechanisms: non-thermal – relativistic processes – methods: numerical – gamma-ray burst: general

**1. Introduction**

Relativistic outflows are common among astrophysical phenomena featuring accretion and/or explosions. Such outflows are observed as sources of bright non-thermal emission, indicating conversion of their kinetic energy into radiation. Internal shocks, collisions between parts of the outflow with different velocities, are one of the proposed dissipation mechanism in many astrophysical contexts. They were first introduced by Rees (1978) to explain “knots,” which are resolved inhomogeneities in the jet of galaxy M87 and were subsequently invoked to explain emission of radio-loud quasars (Spada et al. 2001), microquasars (Kaiser et al. 2000; Malzac 2014), and GRBs (Rees & Mészáros 1994; Sari & Piran 1997; Kobayashi et al. 1997; Daigne & Mochkovitch 1998). Internal shocks arise naturally when assuming inhomogeneities in the (proper) velocity of the outflow. In this scenario, a faster shell catches up with a slower one ejected at an earlier time by the central engine. The shells collide, and under the right conditions (see e.g., Pe’er 2014, for a review) two shock fronts will form: a forward shock (FS) propagating into the slow shell and a reverse shock (RS) propagating into the faster one. Particle acceleration is expected

at these shock fronts, which consequently powers the emission from the outflow.

A common approach used to calculate the dissipation efficiency and calculate the emissions from internal shocks assumes plastic collisions: shells merge inelastically and continue propagating as a single shell (Daigne & Mochkovitch 1998; Beloborodov 2000; Spada et al. 2000; Guetta et al. 2001; Bošnjak et al. 2009; Malzac 2014; Bustamante et al. 2017; Rudolph et al. 2022, 2023). The energy dissipation is associated to a FS if the Lorentz factor (LF) of the fused shell is closer to the fast shell, or to a RS otherwise. While providing a useful approximation that allowed us to reproduce a number of features associated with internal shocks, this approach crudely approximates the location of the shock front. Several studies have focused on the shock physics (e.g., Kino et al. 2004; Pe’er et al. 2017), but without considering the temporal evolution for a generic parameter space. This type of work was done by Rahaman et al. (2024a) (hereafter R24a), who provided an analytical framework for internal shocks between two cold, homogeneous, unmagnetized shells of arbitrary proper velocities in planar geometry.

This hydrodynamical framework is then completed with the parametrization from Genet & Granot (2009) for the optically thin synchrotron emission of a single shock front propagating over a range of radii. They found an analytical solution for the observed flux using integration over the equal arrival time

<footnote n="1">¹ Astrophysics Research Center of the Open University (ARCO), The Open University of Israel, 1 University Road, PO Box 808, Raanana 4353701, Israel</footnote>

<footnote n="2">² Centre de Recherche Astrophysique de Lyon (CRAL), ENS de Lyon, UMR 5574, Université Claude Bernard Lyon 1, CNRS, Lyon 69007, France</footnote>

<footnote n="3">³ Department of Natural Sciences, The Open University of Israel, PO Box 808, Raanana 4353701, Israel</footnote>

<footnote n="4">⁴ Department of Physics, The George Washington University, 725 21st Street NW, Washington, DC 20052, USA</footnote>

<footnote>* Corresponding author: arthur.charlet@ens-lyon.fr</footnote>
surface (EATS) for a Band function spectral shape. Rahaman et al. (2024b) (hereafter, R24b) built their model on this basis, adding refinements such as the ratio between shock front LF and downstream fluid LF, while also considering the contribution from both shock fronts using the estimates obtained in R24a. This allowed them to self-consistently calculate the flux of a single collision at any observed frequency, which can then be used as a building block for full light curves and/or spectra of relativistic outflows. In this framework, the relative positions of the shock fronts and its dependency on the initial conditions are an important factor for the shape of the observed light curve, a feature that tends to get washed away by the ballistic approach. While the choice of initial conditions in most papers exploring the internal shock model for the GRB prompt emission have led to a very short-lived FS of negligible contribution to the total emission, the choice made in R24b led to the contributions from the two fronts covering similar timescales. This allowed the FS to convert enough energy to have a sizable contribution to the emission, producing a two-component time-resolved spectrum. This physical scenario for prompt GRB emission may provide an explanation for the significant deviations from the pure Band function fit (Band et al. 1993) in the time-resolved analysis of the brightest bins of observed bursts (Preece et al. 2014; Vianello et al. 2018).

The hydrodynamical solution is significantly simplified in the planar case, which is a decent approximation as long as the dissipation radius varies by a factor $\le 2$ during the shock propagation. In particular, the shock strengths remain fixed with time or radius in that case. Because of this R24a determined all quantities in planar geometry, and R24b introduced spherical effects on the peak frequencies and luminosities in an approximated manner by neglecting variation of the shock fronts LF and shock strength with time. Considering that such effects may modify the emission signature significantly, a fully consistent spherical approach is necessary to properly quantify the applicability of R24b results.

The present study focuses on the construction of the thin-shell synchrotron emission corresponding to the thin cooling layer behind the shock in the fast-cooling regime of spherical colliding shells through numerical simulations. In Sect. 2, we introduce the physical framework of this study, beginning with a summary of the main results from R24a (Sect. 2.1), followed by a presentation of the numerical code (Sect. 2.2), the model used to derive observed flux from our simulation results (Sect. 2.3), and finishing with the numerical setups for this study (Sect. 2.4). We present the results in Sect. 3, starting with the comparison of our fiducial spherical run with the planar case (Sect. 3.1) before exploring the fully spherical regime over the doubling radius (Sect. 3.2). Our calculated observed flux is presented in Sect. 3.3, then we present a study of the behavior of its peaks (Sect. 3.4) and try to derive information on the source activity from a few selected GRBs using those results (Sect. 3.5). Finally, we offer our first insights into the cooling effects on the observed spectra through an approximated marginally fast-cooling regime (Sect. 3.6). We present our conclusions and discuss our results in Sect. 4.

## 2. Physical scenario and numerical methods

### 2.1. Hydrodynamical framework

In R24a, the authors showed that the collision of two homogeneous cold relativistic shells is determined by seven basic parameters, which are presented Table 1: the time, $t_{\text{off}}$, between the

<table><caption>**Table 1.** Physical parameters of the setup.</caption><thead><tr><th>Symbol</th><th>Definition</th></tr></thead><tbody><tr><td>$E_{k0,i}$</td><td>Available kinetic energy in shell i just before collision</td></tr><tr><td>$\Delta_{0,i}$</td><td>Radial width of shell i just before collision</td></tr><tr><td>$u_i$</td><td>Proper speed of shell i</td></tr><tr><td>$t_{\text{off}}$</td><td>Time between ejection of the shells</td></tr><tr><td>$L_i$</td><td>Source power during the ejection of shell i</td></tr><tr><td>$t_{\text{on,i}}$</td><td>Activity time of the source during ejection of shell i</td></tr></tbody></table>

ejection of the two shells, their proper speeds ($u_1, u_4$), where $u \equiv \Gamma\beta$, along with initial radial widths ($\Delta_{0,1}, \Delta_{0,4}$) and initial kinetic energies ($E_{k0,1}, E_{k0,4}$). Here and elsewhere, the subscript 0 denotes the initial values of properties, namely, at ejection or at the collision radius. Alternatively, the activity time ($t_{\text{on},1}, t_{\text{on},4}$) and the power ($L_1, L_4$) of the source during the emission of the shells may be used in place of the shells width and kinetic energy for a set of parameters focused on the source activity. The two are easily related through $\Delta_{0,i} = \beta_i c t_{\text{on},i}$ and $E_{k0,i} = L_i t_{\text{on},i}$, assuming constant jet power and outflow velocity within each ejection interval. In particular, the assumption of no velocity spread within a shell means that all fluid elements move at the same velocity and its width does not change during propagation. As shown in Table 2, those seven parameters can be combined into four derived parameters required to describe the post-collision shock hydrodynamics: the collision radius, $R_0$, the shells’ radial widths ratio, $\chi$, the proper velocities ratio, $a_u$, and the proper density ratio, $f$. We also introduced the collision time $t_0 = R_0/\beta_4 c$ in the source frame, where we implicitly chose $t=0$ at the ejection of the front of shell 4. The shells were assumed to be cold, $p_1 = p_4 = 0$, and the proper density of shell i was obtained from

$$ \rho'_{i} = \frac{1}{\Gamma_{i}(\Gamma_{i} - 1)} \frac{E_{k0,i}}{4\pi R_{0}^{2} \Delta_{0,i} c^{2}}, \qquad (1) $$

where primes denote quantities measured in the rest frames of the corresponding fluids. The collision of the two shells produces a pair of shocks: a reverse shock (RS) propagating into shell 4 and a forward shock (FS) propagating into shell 1. The two shocked regions, region 3 (shocked shell 4) and region 2 (shocked shell 1), are separated by a contact discontinuity (CD). At the collision radius, the use of the shock jump conditions together with pressure equality across the CD and the equation of state (EoS) in the shocked regions allow for an analytical derivation of all relevant hydrodynamical quantities. Values at the collision radius are noted with the subscript 0.

In planar geometry, quantities are constant across a shocked shell and with propagation. From the constant shock front velocities one easily derives shock crossing times and radii, determining the emission timing properties. In Appendix A, we present the main results from R24a. Those do not hold anymore in spherical geometry, prompting the present numerical study.

### 2.2. Numerical method

This study was performed with the code GAMMA (Ayache et al. 2022) with the aim to solve the equations of special relativistic hydrodynamics (SRHD) in one dimension, using a finite-volume Godunov scheme, the HLLC (Mignone & Bodo 2005) solver for relativistic hydrodynamics, piecewise linear spatial reconstruction, and following an arbitrary Lagrangian-Eulerian approach.
<table><caption>**Table 2.** List of derived parameters.</caption><thead><tr><th>Symbol</th><th>Definition</th><th>Expression</th></tr></thead><tbody><tr><td>R<sub>0</sub></td><td>Collision radius</td><td><sup>β<sub>1</sub>β<sub>4</sub>ct<sub>off</sub></sup>/<sub>β<sub>4</sub> − β<sub>1</sub></sub></td></tr><tr><td>χ</td><td>Radial width ratio</td><td><sup>Δ<sub>0,1</sub></sup>/<sub>Δ<sub>0,4</sub></sub></td></tr><tr><td>a<sub>u</sub></td><td>Proper speed ratio</td><td><sup>u<sub>4</sub></sup>/<sub>u<sub>1</sub></sub> &gt; 1</td></tr><tr><td>f</td><td>Proper density ratio</td><td><sup>n'<sub>4</sub></sup>/<sub>n'<sub>1</sub></sub> = χ <sup>E<sub>k0,4</sub>Γ<sub>1</sub>(Γ<sub>1</sub> − 1)</sup>/<sub>E<sub>k0,1</sub>Γ<sub>4</sub>(Γ<sub>4</sub> − 1)</sub></td></tr></tbody></table>

This means **GAMMA** can compute fluxes for any interface velocity. The HLLC solver adds a calculation of the contact discontinuity (CD) wave speed to the two-wave HLL solver (Harten et al. 1983), the default behavior of **GAMMA** sets the interface velocity to that of the CD. We used this default setting throughout the entirety of this work, as a mesh moving with the flow's velocity evolves propagating shells over a wide range of scales using bigger time steps: in such a mesh, the limiting velocity becomes the sound speed; whereas in fixed-mesh approaches, it would be the flow speed instead. Such a moving mesh also offers the added benefit of Lagrangian behavior of the cells like natural refinement of zones with high gradients, as the cell size will follow the compression of the fluid.

**GAMMA** offers the choice between several time integration methods, among which we chose the third-order Runge-Kutta. The time step is adaptive, based on a Courant-Friedrich-Lewy (CFL) condition (Courant et al. 1928). To be consistent with the derivation of R24a, we implemented the Taub-Mathews equation of state (Mathews 1971) following Mignone & McKinney (2007). We used none of the adaptive mesh refinement methods present in **GAMMA** because we wanted to be able to properly identify cells from one time step to the next and compare their properties. Finally, we use the shock detection algorithm introduced in Zanotti et al. (2010) already implemented in **GAMMA**. We found setting the shock detecting threshold to 0.15 produces satisfactory detections of the two fronts across the simulation. **GAMMA** also contains a method to inject an electron distribution in shocked cells and let them evolve using a reformulation of the cooling equation into an advection equation, which we did not use in this work given our choice of assumptions for the radiative efficiency (see discussion in Rahaman et al. 2024b, and references therein). We leave the exploration of different cooling regimes to a future work.

## 2.3. Flux calculation from cells

In this study, we derived the observed flux from our simulations by applying the method for an infinitely thin shell, as described in Genet & Granot (2009), to cells right downstream of each shock front. Here, the infinitely thin shell approximation is validated as $\Delta r/R \sim 10^{-9}$ for shocked cells in the simulations run for this work. The contributing cells are chosen right downstream of the shocks, counting only one contribution of each cell crossed by the shock, waiting the few time steps necessary for the shock to cross to another cell to add a new contribution. We also assumed the accelerated electrons follow a power-law energy distribution of index $p$ ($dN_e/d\gamma_e \propto \gamma_e^{-p}$ for $\gamma_e \ge \gamma_m$) and all the energy given to the electrons of a cell by the shock passage is radiated in less than a numerical time step. This requires

the emission to be deep in the “fast-cooling” regime. After identifying the contributing cell right downstream of the shock, we derived the minimum Lorentz factor of the post-shock electron distribution by normalizing the total available energy over the electron population:

$$ \gamma_m = \frac{p-2}{p-1} \frac{m_p}{m_e} \frac{\epsilon_e}{\xi_e} \frac{e'_{\text{int}}}{\rho' c^2}, \quad (2) $$

for $p > 2$, where $e'_{\text{int}}$ is the comoving (or proper) internal energy density; $\epsilon_e$ is the fraction of internal energy transferred to the fraction $\xi_e$ of total electrons. For this study, we assumed an equipartition of energy $\epsilon_e = \epsilon_B = 1/3$ and chose values of $p = 2.5$, $\xi_e = 10^{-2}$ for the accelerated electron distribution power-law index and fraction of accelerated electrons respectively. To this Lorentz factor corresponds the comoving peak frequency

$$ v'_{\text{m}} = v'_{\text{syn}}(\gamma_{\text{m}}) = \frac{q_e B'}{2\pi m_e c} \gamma_{\text{m}}^2. \quad (3) $$

The contribution of a cell at radius $r$ and time $t$ traveling with dimensionless velocity $\beta$ (corresponding to the Lorentz factor $\Gamma$) to the flux at an observed frequency, $\nu$, and observed time, $T$, is:

$$ F_{\nu}(\tau \ge 1) = \frac{1+z}{4\pi d_{\text{L}}^2} \tilde{L}_{\nu_{\text{m}}} \tau^{-2} S\left(\tau \frac{\nu}{\nu_{\text{m}}}\right), \quad (4) $$

where $S$ is a normalized function verifying $S(x) = xS(x) = 1$ for $x = 1$. The normalized time, $\tau$, is defined as

$$ \tau = \frac{T - T_{\text{ej}}}{T_{\theta}}, \quad T_{\text{ej}} = (1 + z) \left( t - \frac{r}{\beta c} \right), \quad T_{\theta} = \frac{(1 + z)r}{2\Gamma^2 c}, \quad (5) $$

and the peak luminosity is

$$ \tilde{L}_{\nu_{\text{m}}} = \frac{\Gamma \Delta V^{(3)}}{T_{\theta}} \frac{\epsilon_e e'_{\text{int}}}{W(p) \nu'_{\text{m}}}. \quad (6) $$

Here, $\Delta V^{(3)}$ is the (three-)volume of the cell in the rest frame of the central source and $W(p) = 2(p-1)/(p-2)$. We obtained Eq. (6) by comparing the isotropic energy from a single pulse in the formulation of GG09 and the formulation of De Colle et al. (2012). The derivation is detailed in Appendix B.2. The quantity $\tilde{L}_{\nu_{\text{m}}}$ is a numerical equivalent luminosity derived from synchrotron power and is not to be confused with the usual luminosity that appears in similar equations for the flux in GG09 or R24b. From their Appendix D, we can express the comoving luminosity behind a shock front:<fnref n="1" />

$$ L'_{\nu'_{\text{m}}} = \frac{L'_{\text{bol}}}{W(p)\nu'_{\text{m}}} = \frac{4}{3}\pi R^2 \beta_{\text{ud}} c \frac{\epsilon_e e'_{\text{int}}}{W(p)\nu'_{\text{m}}}, \quad (7) $$

using $e'_{\text{int}} = (\Gamma_{\text{ud}} - 1)\rho_{\text{d}}c^2$ to obtain $\beta_{\text{ud}}$. Finally we identified:

$$ \tilde{L}_{\nu_{\text{m}}} = \frac{3\Gamma\Delta r}{\beta_{\text{ud}}cT_{\theta}} L'_{\nu_{\text{m}}}. \quad (8) $$

In this work, the spectral shape, $S$, will either be the synchrotron broken power law with low- and high-frequency spectral slopes, $b_1$ and $b_2$, respectively (syn-BPL), or the Band function:

$$ S(x) = e^{1+b1} \begin{cases} x^{b_1} e^{-(1+b_1)x} & x \le x_b, \\ x^{b_2} x_b^{b_1-b_2} e^{-(b_1-b_2)x} & x \ge x_b, \end{cases} \quad (9) $$

<footnote n="1">¹ Since the electrons are assumed by construction to be deep in the fast-cooling regime, it follows that the bolometric luminosity is directly related to $\nu_m$ and the spectral luminosity at $\nu_m$ and is independent of the cooling frequency.</footnote>
where $x_b = (b_1 - b_2)/(1 + b_1)$. Following R24b, the choice of setting both fronts deep in the fast-cooling regimes sets values of $b_1 = -1/2$ and $b_2 = -p/2 = -1.25$ for spectral slopes. This choice is further discussed in the conclusion. The effect of setting the weaker FS in the marginally fast-cooling regime is discussed in Sect. 3.6. Summing all contributions from the run gives the total observed flux as a function of observed time and frequency.

In the following, quantities are plotted in units of the normalization frequency, $v_0$, and flux, $F_0$, defined as the peak observed frequency and flux radiated at the RS at collision in R24b:

$$v_0 = 2\Gamma_0 v'_{m,RS}(R_0), \qquad (10)$$

$$F_0 = \frac{(1+z)}{12\pi d_L^2} 2\Gamma_0 L'_{v'm,RS}(R_0). \qquad (11)$$

Here, $\Gamma_0$ is the Lorentz factor of the shocked material at $R_0$, which is at this radius equal downstream of both shock fronts. We refer to their Appendix G for the full derivation. Here, we write $\tilde{v} = v/v_0$ as the normalized frequency and $\tilde{T} = 1 + \tilde{T} = (T - T_{ej,0})/T_{\theta,0}$ as the normalized time – as defined in Eq. (5) with the relevant values at $t_0$ and $R_0$.

## 2.4. Numerical setup

To expand on the theoretical framework of R24b, we explored the collision of two ultrarelativistic shells of equal energy and width with a moderate proper velocity ratio of $a_u = 2$. The kinetic energy available in the shells is taken to be $10^{52}$ erg, a typical order of magnitude value for the isotropic equivalent value corresponding to a single GRB pulse. The choice of equal activity time and equal shell width is doubly motivated by the GRB spectrum observations and to ensure a high radiative efficiency. From R24b, the peaks' relative prominence is tied to the shells' sizes: if $t_{on,1} \ge 2t_{on,4}$, the FS peak is too prominent relative to the RS and vice versa (see e.g., their Fig. 4). They also showed that the rarefaction wave traveling back towards the center after shock crossing may catch up with and suppress the second shock before it finishes crossing its shell. This is valid for values of $\chi$ that are not close to 1, limiting the radiative efficiency of the process in such a case. The off time between pulses, $t_{off}$, is set to the same value as the pulses themselves, as observations suggests a correlation between pulse width and interval between pulses, outside of quiescent periods (see e.g., Nakar & Piran 2002). We set both activity times, $t_{on}$, and the off time to a typical order of magnitude for the activity timescale of 0.1 s. The shells are ultrarelativistic with Lorentz factors of $u_1 = 100$ and $u_4 = 200$, which are minimum values to ensure the assumption of optically thin emission.

To complete this setup, we defined the external medium density behind shell 4, $\rho'_{ext,b}$ and in front of shell 1, $\rho'_{ext,f}$ from a density contrast parameter of $\tilde{\rho}'_{ext} = \rho'_4/\rho'_{ext,b} = \rho'_1/\rho'_{ext,f}$. Both regions of the external medium were set with the same velocity as the shell they are in contact with. Such an “external” medium is required for the rarefaction wave to set after shock crossing and hardly affects the results. The external pressure is defined by introducing the relativistic temperature, $\Theta_0$, as a parameter, we set the pressure to be constant everywhere, namely, as $p_0 = \min(\rho'_1, \rho'_4)\Theta_0 c^2$. Setting the pressure equally over the whole simulation box avoids any unwanted effects from pdV working between the external medium and the shells. All these parameters are summarized in Table 3. The simulation was run in 1D over 6100 cells, of which there are $N_{sh} = 3000$ for each shell and 50 on each side for the external medium. This choice of the resolution was meant to ensure that the shock fronts have

<table><caption>**Table 3.** Run parameters in CGS units.</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>E<sub>k,0</sub></td><td>10<sup>52</sup> erg</td></tr><tr><td>u<sub>1</sub></td><td>100</td></tr><tr><td>u<sub>4</sub></td><td>200</td></tr><tr><td>t<sub>on,i</sub></td><td>0.1 (0.5) s</td></tr><tr><td>t<sub>off</sub></td><td>0.1 s</td></tr><tr><td>N<sub>sh</sub></td><td>3000 (15 000)</td></tr><tr><td>ρ'<sub>1</sub></td><td>4.7 × 10<sup>-12</sup> g cm<sup>-3</sup></td></tr><tr><td>ρ'<sub>4</sub></td><td>1.2 × 10<sup>-12</sup> g cm<sup>-3</sup></td></tr><tr><td>Θ<sub>0</sub></td><td>5 × 10<sup>-5</sup></td></tr><tr><td>ρ̃'<sub>ext</sub></td><td>5 × 10<sup>-2</sup></td></tr><tr><td>Δ<sub>0,i</sub></td><td>3 × 10<sup>9</sup>(1.5 × 10<sup>10</sup>) cm</td></tr><tr><td>t<sub>0</sub></td><td>2.7 × 10<sup>3</sup> s</td></tr><tr><td>R<sub>0</sub></td><td>8 × 10<sup>13</sup> cm</td></tr><tr><td>a<sub>u</sub></td><td>2</td></tr><tr><td>χ</td><td>1</td></tr><tr><td>f</td><td>0.25</td></tr></tbody></table>

not propagated more than $2 \times 10^{-3} R_0$ before the downstream values are properly established. This can be estimated by applying the formula for the crossing time (Eq. (A.3)) to the width of five numerical cells, greater than the typical numerical size of the detected shock front with our choice of threshold parameter. Finally, we set an additional passive tracer with different values in each shell to help separate them easily during post-processing.

To fully explore the effects of spherical geometry on shock hydrodynamics by going well above the doubling radius for the crossing radius of each front, we ran another simulation in spherical geometry with larger activity times (and initial shells width) by a factor of 5. Meanwhile, we kept the off time constant. We increased the number of numerical cells per shell by the same factor. The values corresponding to this run are given in brackets in Table 3.

## 3. Effects of spherical geometry

Having calibrated our simulation against the analytical results in planar geometry in Appendix A, we explored the spherical effects through two simulations. In the first, all parameters are equal for the sake of comparison and in the second, the activity times $t_{on,i}$ (and, thus, the shell width $\Delta_{0,i}$) have been multiplied by 5 to explore hydrodynamics over a larger range of radii.

### 3.1. Spherical effects on shells structure and shock dynamics

A first result of performing the same run in spherical geometry is the variation in crossing times (radii): for the same values of parameters, $\Delta R_{RS}/R_0$ increases from 1.29 to 1.52 when switching from planar to spherical geometry, while $\Delta R_{FS}/R_0$ only grows from 1.62 to 1.67. Figure 1 displays snapshots of the comoving density, proper velocity, and pressure (similarly to Fig. A.1) at times $t_0$, 2 $t_0$, and 2.6 $t_0$. In the analytical expectations from planar geometry, the density and pressure are rescaled by $(R/R_0)^{-2}$ to account for propagation, then plotted to better highlight the differences with R24a. Spherical effects modify the structure of the shocked regions: the density profile shows compression from the shocks to the contact discontinuity, a proper velocity decrease with radius, and a slightly increasing pressure profile from the RS to the FS. As the shells and shocks propagate,
<figure><img src="image_2.png" /><figcaption>Fig. 1. Snapshots of the run in spherical geometry, with all parameters similar to the fiducial run, at $t = t_0$ (left), before first (reverse) shock crossing (center), and after shock crossing (right). Dashed black lines show the analytical expectations from the results in planar geometry. The legend is the same for all panels.</figcaption></figure>

the density in regions 1 and 4 (i.e., the unshocked portions of both shells) will decrease with $r^{-2}$. This is due to mass conservation of a fluid element in spherical geometry, with our assumption of little to no velocity spread between the two interfaces. At the zeroth order, where we consider constant shock strength, the downstream density, $\rho_d$, and pressure will also evolve following a $r^{-2}$ law. From pressure continuity (i.e., the pressure gradient is rather small in the shocked regions as they are in causal contact), the pressure at the CD will follow the same law instead of the $r^{-2\hat{\gamma}}$ from adiabatic expansion, with $\hat{\gamma}$ being the adiabatic index. Specifically, $p_d \approx p_0(R/R_0)^{-2} \approx p_{CD} \approx (\rho_{CD}/\rho_0)^{\hat{\gamma}}$ and $\rho_d = \rho_0(R/R_0)^{-2}$, implying $\rho_{CD}/\rho_d = (R/R_0)^{2(1-\hat{\gamma})}$. This results in a compression of the fluid towards the center with the ratio between densities behind the shock and at the CD following a $-2(1-\hat{\gamma})$ law. This implies a density scaling $\rho_{CD} \propto r^{-2/\hat{\gamma}}$ between $r^{-1.5}$ and $r^{-1.2}$ ($\rho_{CD}/\rho_d \propto r^{2(1-\hat{\gamma})}$ between $r^{0.5}$ and $r^{0.8}$) for an adiabatic index varying between 4/3 and 5/3.

Both of these effects are seen Fig. 2 (top panel) which displays densities rescaled with $(r/R_0)^2$ to take out spherical effects: the values taken right downstream of the shocks show very little to no variation in radius, while those on both sides of the CD increase with the same $r^{0.7}$ (or $\rho' \propto r^{-1.3}$ without rescaling) fitting with our observed regime of intermediate shock strengths. This results in intermediate adiabatic index values according to the Kumar & Granot (2003) formulation used in R24a:

$$ \hat{\gamma} = \frac{4\Gamma_{ud} + 1}{3\Gamma_{ud}}, \qquad (12) $$

which those authors demonstrated is equivalent to the Taub-Mathews EoS for a cold upstream medium. To a highest order of precision, the variation of the shock strength $\Gamma_{ud} - 1$ with propagation changes this behavior; as seen in the bottom panel of Fig. 2, the shock strengths decrease with propagation at a stronger rate for the RS than the FS. This translates to a slight decrease in the post-RS rescaled density and stronger decrease in the post-RS pressure than that seen downstream of the FS (as shown in the 3rd panel of Fig. 2). This explains the increasing pressure profile from RS to FS seen in the snapshots. Similarly, we observed a higher post-RS flow velocity than that seen downstream of the FS, where the first increases with propagation, while the second decreases with it, creating the profile observed Fig. 1. In terms of power-law scalings, this approach has proven useful to explain behavior below the doubling radius

(i.e., $R \lesssim 2R_0$); however, when exploring a broader range of radii, a different behavior is observed above it.

## 3.2. Asymptotical hydrodynamics in the spherical regime

We present the results of the simulation run with both activity times multiplied by a factor of 5, while keeping $t_{off}$ and all other quantities constant. This means the shells are now wider than their separation by this same factor, allowing us to explore the evolution of the various relevant quantities for the flux calculation over a range of radii ~10 $R_0$, almost an order of magnitude greater than the situation presented in Sect. 3.1. In Fig. 3, we show the fluid and shock fronts Lorentz factors in the top panel and the shock strengths in the bottom panel. The data were smoothed by a rolling average window to eliminate the strong noise caused by the propagation of numerical errors in the interface velocity over the large number of code iterations for this simulation. Such oscillations could be reduced by the use of an higher spatial order reconstruction algorithm instead of the first-order piecewise linear algorithm currently present in GAMMA. The downstream and shock Lorentz factors evolve at the same rate, keeping the ratio $g = \Gamma/\Gamma_{sh}$ introduced in R24a close to constant. At the collision radius, the values are in accordance to the analytical expectations in planar geometry derived in R24a, before decreasing in a power-law up to the doubling radius and smoothly connecting to a constant value after a few $R_0$. Values relative to the RS change by a greater amount in this asymptotical regime compared to the planar case than values relating to the FS: the shock and downstream Lorentz factor grow by ~10% at the RS and go down by ~3% at the FS, while the shock strength diminishes by ~3% and ~0.7%, respectively. Thus, we modeled any hydrodynamical or derived quantity $X$ with the law:

$$ X(R) = \left[ (X_{0*} \tilde{R}^n)^s + (X_{\text{sph}} \tilde{R}^h)^s \right]^{1/s} = X_0 f_X(R), \qquad (13) $$

with $\tilde{R} = R/R_0$ and $f_X$ a function verifying $f_X(R_0) = 1$. The index $h$ corresponds to the expected scaling from propagation effects at constant shock strength, which can be 0. Then, $\tilde{X}_{\text{sph}} \tilde{R}^h$ is the asymptotical value at large radii and $X_{0*}$ is defined by the value at collision radius $X_0 = X(R_0) = [X_{0*}^s + X_{\text{sph}}^s]^{1/s}$. We note that $s$ values need to be of the opposite sign of $n$ (and of $h$ when nonzero). The values obtained for the shock strength and downstream Lorentz factors by the fitting procedure are given Table 4.
<figure><img src="image_3.png" /></figure>

**Fig. 2.** Top to bottom: comoving densities, Lorentz factors of the fluid and shock fronts, pressure, and shock strengths. The hydrodynamical quantities are measured downstream of shocks ('d,3' downstream of the RS, 'd,2' for the FS), at the CD by averaging over a few cells crossing the interface, or on each side of it (CD$_3$, CD$_2$) for the density. Legend is the same for all panels. The derived Lorentz factor of the shock fronts is displayed with the measured ones. Densities and pressures are in units of the values behind the RS in the planar case and rescaled to $(r/R_0)^2$. We smoothed the numerical oscillations in data by a rolling average. Horizontal dashed lines show the corresponding analytical values from R24a.

While the change in those quantities and thus the related quantities such as luminosity and peak frequency are not very significant in absolute values, it is their rate of change that has the most impact on the observable data at very early times.

### 3.3. Emission from two spherical colliding shells

In Figs. 4 and 5, we present the resulting light curves and time-integrated spectra from our simulations, comparing the two geometries in the same panel for both the Band and the broken power-law spectral shapes. To highlight the difference between the assumptions for propagation chosen in R24b and the fully spherical results, we compare our results to the hybrid approach instead of the fully planar case detailed in Appendix B.3. The total light curve is weakly affected by the geometry, with slightly increased peak times corresponding to our estimation in Sect. 3.1 and a smoothed shape compared to the planar case. In comparison, the individual contributions of each shock front vary significantly, with the RS peak luminosity increasing by a factor of 1.25 and that of the FS decreasing by a factor of 0.7.

<figure><img src="image_4.png" /></figure>

**Fig. 3.** Top: Lorentz factors of the shock fronts, the fluid downstream of the shocks and at the CD with radius. Bottom: evolution of shock strengths with radius. Data were smoothed by rolling average. The quantities converge smoothly to a constant value after $r \approx 4R_0$, shown by the dotted horizontal line for the quantities relative to the RS and the FS. Legend is the same as in Fig. 2.

The time-integrated spectra in Fig. 5 show the peak frequencies of both shock fronts decrease by different amounts: the quantities relative to the RS decrease faster with radius than those relative to the FS. Thus, the observed range of the frequencies ratio between the two spectral components can be attained with a wider range of $a_u$ in this fully spherical approach. From the hybrid approach, we retained the doubly broken power-law spectrum with indices $1 + b_1 = 0.5$ and $1 + b_2 = -0.25$ at low and high energies, directly related to our choice of synchrotron spectral slopes, and an intermediate part $1 + b_{\text{mid}} \approx 0.1$. The effective middle slope $b_{\text{mid}}$ depends on the flux normalization for each front and is thus very variable with the choice of hydrodynamical conditions: an analytical estimate using R24b show $1 + b_{\text{mid}} \in [-0.05, 0.45]$ when exploring the parameter space $(\chi, a_u) \in [-0.5, 1.5] \times [1.01, 5]$, under the assumptions of ultrarelativistic shells and constant central engine power. These slopes are the same for both the Band and syn-BPL spectral shapes. The differences between the results of this present study and those from R24b are highlighted by redrawing in dotted line the Band + hybrid light curve (respectively, the time-integrated spectrum) in the syn-BPL panel. The main change lies in the total received flux (respectively, the fluence), especially at lower frequencies, but the essential characteristics of the light curve (respectively, the time-integrated spectrum) are similar.

### 3.4. Peak frequency and flux

The individual variations of peak frequency and flux associated to each shock front are displayed in Fig. 6 against normalized time, where the spherical effects systematically lowers both the peak frequency and flux for both shocks. The spherical geometry also introduces a difference in the time evolution of the peak frequency between the rising and the high-latitude emission (HLE) part. While the flux rises, the peak frequency decreases faster than $v_{\text{pk}} \propto \dot{T}^{-1}$ expected from purely geometric considerations.
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Light curves in the spherical (full lines) and hybrid (planar geometry and spherical assumptions for peak frequency and luminosity, dashed-dotted lines) cases at fixed frequency $\tilde{\nu} = 0.1$. Top panel: light curves with the Band's spectral shape. Bottom panel: broken power-law spectral shape. The total light curve in the hybrid + Band case is redrawn in the bottom panel (dotted line) for the sake of comparison. The individual contributions from each shock change but the resulting light curve do not, apart from the slightly increased breaks from the increase in crossing times. The choice of spectral shape decreases total emission at this frequency, while keeping the same time evolution.</figcaption></figure>

This is especially the case for the RS, with $v_{pk,RS} \propto \tilde{T}^{-1.9}$. In the following, we present a way to estimate those quantities from the hydrodynamics with fewer calculations than the full flux calculation (detailed in Sect. 2.3). We built this approach on the asymptotical regimes identified for the hydrodynamical quantities in Sect. 3.2.

Assuming the hydrodynamical variables can be modeled according to the law given by Eq. (13), after extracting the relevant variables for flux calculation from the simulation results, we fit them with the same type of law:

$$ \Gamma_d(R) = \left[ (\Gamma_{0*} \tilde{R}^{-m/2})^s + (\Gamma_{sph})^s \right]^{1/s} \equiv \Gamma_{f,0} f_{\Gamma}(R), \quad (14) $$

$$ v'_{m}(R) = \left[ (v'_{0*} \tilde{R}^{d})^{s} + (v'_{sph} \tilde{R}^{-1})^{s} \right]^{1/s} \equiv v'_{f,0} f_{v'}(R), \quad (15) $$

$$ L'_{v'm}(R) = \left[ (L'_{0*} \tilde{R}^{a})^{s} + (L'_{sph} \tilde{R}^{s})^{s} \right]^{1/s} \equiv L'_{f,0} f_{L'}(R). \quad (16) $$

We performed the fit with the *curve_fit* function from SciPy, which is based on a non-linear least square method. We give in Table 4 the obtained values for the fits. The numerical values at $R_0$ differ slightly (≤2%) from the analytical expectations, we note them with the subscript “f, 0” to avoid confusion.

Then, to obtain an estimate for the peak frequency and flux in the rising part of the prompt emission, we assume we can define an effective radius on the EATS at which the peak contribution can be calculated $R_{\text{eff}}(T) = y_{\text{eff}}R_{\text{L}}(T)$, with $R_{\text{L}}(T)$ as the maximal radius of the EATS (found along the line of sight). At a given observer time, $T$, this radius can also be identified by its angle with the line of sight – or, more conveniently, by the quantity $\xi_{\text{eff}} = (\Gamma\theta)_{\text{eff}}^2$. A schematic of this method is given Fig. C.1.

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Time-integrated spectra following the same structure as in Fig. 4. The spherical effects decrease the peak frequencies for both shocks, with a stronger effect on the RS. The peaks from both contributions are closer than in the planar case, resulting in a three-part spectrum. This effect is more visible with the broken power-law spectral shape where the two peaks are distinct from each other in the hybrid approach. The choice for the spectral shape also diminishes the flux at low frequencies.</figcaption></figure>

<table><caption>**Table 4.** Fitting parameters for emission-related quantities.</caption><thead><tr><th></th><th>X<sub>0*</sub></th><th>X<sub>sph</sub></th><th>Index</th><th>|s|</th></tr></thead><tbody><tr><td>(Γ<sub>ud</sub> − 1)<sub>RS</sub></td><td>0.103</td><td>0.0735</td><td>n = −0.64</td><td>5.0</td></tr><tr><td>(Γ<sub>d</sub>/Γ<sub>0</sub>)<sub>RS</sub></td><td>1.01</td><td>1.08</td><td>m = −0.29</td><td>24.0</td></tr><tr><td>(v'<sub>m</sub>/v'<sub>0</sub>)<sub>RS</sub></td><td>0.72</td><td>0.40</td><td>d = −3.46</td><td>1.0</td></tr><tr><td>(L'<sub>v'm</sub>/L'<sub>0</sub>)<sub>RS</sub></td><td>1.03</td><td>1.43</td><td>a = 1.84</td><td>3.1</td></tr><tr><td>(Γ<sub>ud</sub> − 1)<sub>FS</sub></td><td>0.0278</td><td>0.0236</td><td>n = −0.21</td><td>18.5</td></tr><tr><td>(Γ<sub>d</sub>)<sub>FS</sub></td><td>0.999</td><td>0.982</td><td>m = 0.06</td><td>102.2</td></tr><tr><td>(v'<sub>m</sub>/v'<sub>0</sub>)<sub>FS</sub></td><td>0.40</td><td>0.65</td><td>d = −2.67</td><td>1.0</td></tr><tr><td>(L'<sub>v'm</sub>/L'<sub>0</sub>)<sub>FS</sub></td><td>1.05</td><td>1.16</td><td>a = 1.39</td><td>5.7</td></tr></tbody></table>

Both quantities are joined by $y = [1 + g^{-2}\xi]^{-1}$ and we take $R_L = R_0 \tilde{T}$, where we implicitly choose a constant $\Gamma_{\text{sh}}$ to approximate the EATS at the zeroth order. Then, the peak frequency and flux are given by:

$$ (1+z)v_{\text{pk}}(T) = k_v v_{f,0} \frac{f_{\Gamma}(R_{\text{eff}}) f_{v'}(R_{\text{eff}})}{1 + \xi_{\text{eff}}}, \quad (17) $$

$$ F_{v_{\text{pk}}}(T) = k_F F_{f,0} 3g^2 \xi_{\text{max}} \left( \frac{f_{\Gamma}(R_{\text{eff}})}{1 + \xi_{\text{eff}}} \right)^3 \frac{f_{L'}(R_{\text{eff}})}{f_{\Gamma}(R_{\text{L}})^2}. \quad (18) $$

Here, $v_{f,0} = 2\Gamma_{f,0}v'_{f,0}$ and $F_{f,0} = ((1+z)/(4\pi d_L^2))2\Gamma_{f,0}L'_{f,0}$. $k_v$ and $k_F$ are two free parameters of order unity introduced to compensate for the difference between analytical and numerical values at $R_0$ such that $k_v v_{f,0} = v_0$ and $k_F F_{f,0} = F_0$. The observed peak frequency is obtained by calculating the Doppler factor and comoving peak frequency at $R_{\text{eff}}$, and the flux is obtained by assuming the flux integral (Eq. (B.7)) can be approximated by the value of the integrand at $\xi_{\text{eff}}$ times the extension of the EATS in this variable, $\xi_{\text{max}} = g^2(\tilde{T} - 1)$. We then looked for the function
<table><caption>**Table 5.** Fitting parameters for the peak frequency and flux.</caption><thead><tr><th>Front</th><th>k<sub>ξ</sub></th><th>ξ<sub>sat</sub></th><th>s</th><th>k<sub>ν</sub></th><th>k<sub>F</sub></th><th>d<sub>ν</sub></th><th>d<sub>F</sub></th></tr></thead><tbody><tr><td>RS</td><td>0.29</td><td>6.98</td><td>1.0</td><td>0.98</td><td>0.54</td><td>-1.31</td><td>-2.76</td></tr><tr><td>FS</td><td>0.51</td><td>2.89</td><td>1.0</td><td>0.97</td><td>1.28</td><td>-1.02</td><td>-2.47</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Normalized time evolution of peak frequency and flux for both shock fronts in the spherical (full line) and hybrid (dash-dotted) cases, calculated from the flux obtained with the Band spectral shape. In spherical geometry, the peak frequency decrease at a faster rate during the rising phase and the flux reaches saturation before the break to the high-latitude emission.</figcaption></figure>

$ξ_{eff}(T)$ fitting the data. We note that contributions from regions of the EATS beyond the relativistic beaming angle $1/Γ$ become increasingly suppressed, thus, we defined a saturation angle $ξ_{sat}$ on the order of unity and we chose the following fitting function for $ξ_{eff}$:

$$ξ_{eff}(T) = k_ξ (ξ_{max}^{-s} + ξ_{sat}^{-s})^{-1/s} = k_ξ ([g^2(\tilde{T} - 1)]^{-s} + ξ_{sat}^{-s})^{-1/s}, \quad (19)$$

where $k_ξ$, $ξ_{sat}$, and $s$ are free parameters. The high-latitude emission part is modeled as the emission from the final effective radius $R_{eff,f} = y_{eff}(T_f)R_f$, with $T_f$ the observed time corresponding to the shock crossing:

$$v_{pk,HLE}(T) = v_{pk}(T_f)τ_f^{d_v}, \quad (20)$$

$$F_{v_{pk,HLE}}(T) = F_{v_{pk}}(T_f)τ_f^{d_F}, \quad (21)$$

and $τ_f$ is the normalized time defined at $R_{eff,f}$ following Eq. (5). We give in Table 5 the values found to fit the peak frequency and flux (from Fig. 6) with an accuracy within a few percent.

### 3.5. Maximal frequency and crossing radius

Among the spectral signatures, the behavior of the peak flux versus peak frequency shown in the left panel of Fig. 7 presents the most striking change between planar and spherical dynamics. In this figure, we compare the purely spherical case to the hybrid

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Left: $νF_ν$ at peak frequency vs peak frequency for the reverse shock in the spherical and hybrid case. The fully spherical case shows peak flux over a larger range of frequencies. Right: ratio between break frequency $ν_{bk}$ and frequency at half flux $ν_{1/2}$. Legend is the same for both panels.</figcaption></figure>

approach of R24b, which features planar dynamics along with spherical geometry for calculating the radiation. The spherical case transitions to the anticipated correlation of high-latitude emission ($νF_ν$)$_{pk} \propto ν^3$ at a lower frequency but presents a high flux over a wider range. We define $ν_{bk}$ as the peak frequency at this break. Comparing to Fig. 11 in R24a, the combined effects from spherical geometry on the peak flux and frequency is greater than a pure increase of $ΔR/R_0$ in the hybrid approach; simply increasing the crossing radius in the R24 formalism is not enough to obtain the curves in the spherical case.

In the right panel of Fig. 7, we display the evolution with a normalized crossed radius, $ΔR/R_0$. of the ratio between $ν_{bk}$ and the frequency $ν_{1/2}$, defined as the point where the peak flux in the rising part is half that at $ν_{bk}$. This allows us to add our modifications to the results presented in Table 11 from R24b: from the ratio in frequencies for peak flux obtained from the data presented in Yan et al. (2024) and assuming the RS is the main contribution to the observed flux, they inferred the radius crossed by the RS during the burst $ΔR/R_0$ and, thus, the ratio of activity time to off time $t_{on,4}/t_{off}$. We present in Table 6 our updated values for the crossed radius in light of spherical effects and compare them to the values inferred in R24b. On average, the obtained crossed radius for our fiducial spherical case is half of the value obtained by R24b and about a third for the obtained activity over the off time ratio.

### 3.6. Dissipation efficiency and marginally fast-cooling regime

The total efficiency, $ε_{tot} = ε_{int}ε_eε_{rad}$, of the internal shock process is easily derived in this framework. When it is deep in the fast-cooling regime, it implies a radiative efficiency of $ε_{rad} = 1$ behind both shock fronts. Here, $ε_{int}$ is the conversion rate of kinetic energy into internal energy, commonly called the dissipation efficiency (or “thermal efficiency,” noted $ε_{th}$, in R24a,b). We calculate it by summing the internal energy from each contributing cell and comparing this quantity to the initial available
<table><caption>**Table 6.** Estimation of $\Delta R/R_0$ and corresponding $t_{on}/t_{off}$ for a choice of GRBs.</caption><thead><tr><th>GRB</th><th>ν<sub>bk</sub>/ν<sub>1/2</sub></th><th>(ΔR/R<sub>0</sub>)<sub>R24</sub></th><th>(ΔR/R<sub>0</sub>)<sub>sph</sub></th><th>(t<sub>on</sub>/t<sub>off</sub>)<sub>sph</sub></th></tr></thead><tbody><tr><td>140606B</td><td>0.45</td><td>1.66</td><td>0.81</td><td>0.57</td></tr><tr><td>131011A</td><td>0.59</td><td>1.04</td><td>0.51</td><td>0.37</td></tr><tr><td>170607A</td><td>0.36</td><td>2.36</td><td>1.09</td><td>0.74</td></tr><tr><td>151027A</td><td>0.65</td><td>0.84</td><td>0.43</td><td>0.31</td></tr><tr><td>150514A</td><td>0.60</td><td>0.99</td><td>0.50</td><td>0.36</td></tr><tr><td>120326A</td><td>0.70</td><td>0.71</td><td>0.35</td><td>0.26</td></tr><tr><td>190829A</td><td>0.50</td><td>1.40</td><td>0.69</td><td>0.49</td></tr></tbody></table>

kinetic energy. We obtained a dissipation efficiency of 8.5% in the planar case (equal to the analytical value using Eq. (B5) in R24b) and of 7.3% in the spherical case. The total efficiency is thus $\epsilon_{tot} = 2.4\%$ in the spherical case. In R24a,b, the dissipation efficiency is obtained by comparing the total final internal energy present in the shells after shock crossing to the initial available kinetic energy. Both our “local” and their “global” approaches give the same result in the planar case, but the quantities are different in the spherical case due to the shock strengths evolving as they cross their respective shells. Using the definition R24a,b would give a dissipation efficiency of 3.5% in the spherical case.

The spectrum we obtained using this thin-shell assumption displays a low-energy $\nu F_\nu$ slope of $1 + b_1 = 0.5$, indicative from fast-cooling, which does not seem to fit observed bursts. In Ravasio et al. (2019), for instance, their fits using a doubly broken power-law model give values of the low-energy slopes from several bursts that are more consistent with a slow-cooling regime. While a model where cooling is derived consistently is beyond the scope of the current paper, we began to look into this issue by recalculating the contributions from the FS in an approximated marginally fast-cooling regime. For this, we took $b_{1,FS} = 1/3$ instead of $-1/2$ and we applied an efficiency factor of $\epsilon_{rad,FS} = 1/2$ to the computed FS flux as in R24b. This value of $1/2$ is an arbitrary middle ground between $\epsilon_{rad,FC} \sim 1$ in the fast-cooling case and $\epsilon_{rad,SC} \ll 1$ in the slow-cooling case. With this choice, we obtained an overall efficiency of the prompt emission $\epsilon_{tot} = 2\%$ for the spherical simulation. We compare in Fig. 8 the instantaneous spectra at $T = T_0$ ($\bar{T} = 1$) obtained in this marginally fast-cooling FS case with the fast-cooling regime. Both were calculated from the simulation in spherical geometry. For this marginally fast-cooling FS, we did not obtain the slow-cooling low-energy slope of $1+b_{1,FS} = 4/3$, but rather $\sim 0.8$ from the cumulative contribution of the two shock fronts.

Setting the FS in the marginally fast-cooling regime is not enough to obtain low energy slopes more in line with observations from our set of simulations. However, we need to mention that marginally fast-cooling regime has three issues here. First, it requires a very specific set of physical conditions that is contrary to our general approach. Second, if the FS is marginally fast-cooling, the RS is not very deep in the fast-cooling regime, and should show a cooling break near the FS peak photon energy, below which $1 + b_1 = 4/3$ is expected in the total spectrum. Third, we used an approximated way of calculating the radiation in this regime, we expected the thin shell assumption to break and, thus, the effects of a finite size emission region to appear, which are bound to modify emission signatures. Such effects will be the focus of a follow-up paper dedicated to the exploration of the cooling regimes.

## 4. Conclusions

We present a fully spherical, self-consistent numerical approach to the internal shocks model for the prompt GRB emission,

<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Comparison of instantaneous spectra between two fast-cooling shock fronts (full lines) versus a fast-cooling RS and a marginally fast-cooling FS (dash-dotted lines) at $\bar{T} = 1$. The low energy slopes are steeper and the overall flux is lower with a marginally fast-cooling FS.</figcaption></figure>

showing the impact of the spherical hydrodynamics effects on emission, even for crossing radii lower than twice the initial collision radius. Our basic two-shell collision model describes a single spike in the prompt GRB light curve, which usually consists of multiple pulses that in this picture correspond to multiple collisions. After calibrating our methods against the analytical results obtained in Rahaman et al. (2024a,b) using a planar approach for hydrodynamics and approximate spherical effects for the emission, we numerically extend their conclusions to the fully spherical case, while keeping the assumption of an infinitely thin emitting shell. We present the structural and dynamical differences of shocked shells in spherical geometry and highlight the non-negligible effects of the evolution in Lorentz factors and shock strengths on the emission over the shock crossing time. In particular, the quantities associated with each shock front evolve at a different rate due to the differences in shock strengths. We then expect that following the spectral properties of a burst over time will allow us to distinguish this two-zones emission model from a one-zone model. While we assumed equal values of the microphysical parameters ($\epsilon_e, \epsilon_B, \xi_e, p$) between the two shock fronts, any diversity in those parameters will most likely add diversity to the obtained emission between the two shock fronts.

We produced light curves and (time-integrated) spectra with two different spectral shapes, the smooth Band function and the synchrotron broken power law. While the Band function is a phenomenological model of the emission spectrum, the synchrotron broken power law is physically motivated. Calculating flux with the latter results in less flux at lower frequencies, but conserves the general time evolution of the light curve. In both choices of spectral shape, our framework naturally obtains a doubly broken power-law shape for the spectrum, shape that has been successfully used to fit prompt GRB data (see e.g., Burgess et al. 2014; Oganesyan et al. 2017; Ravasio et al. 2019). The validity of the
results obtained with the rough spectral shape compared to using
the empiric smooth Band function motivates the use of this sim-
pler emission function to explore more complicated flows as well
as other cooling regimes with different spectral shapes in fur-
ther works while retaining physical accuracy. The value of the
spectral slope between the two breaks varies significantly with
the choice of initial conditions. While it is not totally consis-
tent with the range of values obtained in Ravasio et al. (2019),
this is not a contradictory result. Those slopes are more depen-
dent on the choices made for the fitting functions than the low-
and high-energy part and a better test of our model would be
to perform fits of observed data from calculated flux across an
expected parameter space.

While the effects of spherical geometry on the shock hydro-
dynamics are not striking when comparing light curves and time-
integrated spectra, they cannot be ignored when studying the
evolution of spectral shape within a single spike in the prompt
GRB light curve. In particular we obtain shapes that are more
similar to the observed ones for a choice of GRBs presented
in Yan et al. (2024) and our inferred crossing radii and source
activity time over off time a peak flux and frequency evolution
are smaller by a factor of 2 and 3, respectively, when compared
to the estimations done by Rahaman et al. (2024b).

This work relies on many highly idealized assumptions,
especially the choice of homogeneous shells in both density
and Lorentz factor. This assumption may be lifted to obtain
more realistic emission signatures of internal shocks by simu-
lating more physically motivated flows. A striking example is
the sharp edge in the peak flux versus peak frequency plot that
originates from the sharp edge of the shells, emission stopping
at once when the shock crosses. Such a sharp edge is not seen
in observed GRBs (Yan et al. 2024). The exploration of internal
shocks in higher dimensions is also expected to introduce addi-
tional effects. Studies on the external shock of GRB outflows
(i.e., the one that form due to the interaction with the external
medium) in 2D have shown the growth of Rayleigh-Taylor insta-
bilities (RTI) significantly modifies the dynamics of the external
reverse shock and causes its emission to peak at a later time (e.g.,
Duffell & MacFadyen 2013; Ayache et al. 2022). The shear flow
and/or turbulence near RTI fingers at the contact discontinuity
may also be source of particle acceleration.

Additionally, the synchrotron spectra used in this work is
obtained in the limiting case where all the non-thermal electron
energy is radiated away in less than a numerical time step. This
limits the emission region to a single numerical cell represent-
ing an infinitely thin shell, behind each shock front. This case
assumes being very deep in the fast-cooling regime behind both
shock fronts, producing spectral slopes that do not agree to the
low energy part of observed spectra fitted with a doubly broken
power law, such as in Ravasio et al. (2019). Their method obtains
low-energy slopes indicative of a slow-cooling regime, which is
naturally not reached in our assumption. A simple calculation
by changing the FS emission to an approximate marginally fast-
cooling regime was not enough to obtain slopes in this range.
However, exploring the impact of other cooling regimes and
extended emission zones in a consistent way using the capacity
of GAMMA to evolve a non-thermal electron distribution could be
a step towards solving this issue and will be the focus of a future
work. Changing the microphysical parameters is one of the ways
to explore other cooling regimes and try to obtain the observed
photon indices. In Bošnjak & Daigne (2014), they find the best
agreement to observations using varying microphysical parame-
ters with shock conditions. Numerically, this could be achieved

by locally computing the value of these parameters following
recipes. GAMMA already proposes one such recipe to compute *p*,
exploring its effects and implementing local evolution of other
parameters may be pursued in a future work.

The implementation for the local electron distribution in
GAMMA only features synchrotron and adiabatic cooling, but
inverse Compton scatterings with the photons produced at both
fronts are expected to modify the distribution in the context of
GRB emission (see, e.g., Nakar et al. 2009 for a comprehen-
sive reference, Jacovich et al. 2021; McCarthy & Laskar 2024;
Pellouin & Daigne 2024 for recent numerical implementations
in the context of GRB afterglows). Implementing proper inverse
Compton and synchrotron self-Compton is a challenging task
that may be tackled in future works.

***Acknowledgements.*** We thank Frédéric Daigne for the insightful discussions and the CRAL for its welcome. PB is supported by a grant (no. 2020747) from the United States-Israel Binational Science Foundation (BSF), Jerusalem, Israel, by a grant (no. 1649/23) from the Israel Science Foundation and by a grant (no. 80NSSC 24K0770) from the NASA astrophysics theory program.

References

Ayache, E. H., van Eerten, H. J., & Eardley, R. W. 2022, MNRAS, 510, 1315

Band, D., Matteson, J., Ford, L., et al. 1993, ApJ, 413, 281

Beloborodov, A. M. 2000, ApJ, 539, L25

Bošnjak, Ž., & Daigne, F. 2014, A&A, 568, A45

Bošnjak, Ž., Daigne, F., & Dubus, G. 2009, A&A, 498, 677

Burgess, J. M., Preece, R. D., Connaughton, V., et al. 2014, ApJ, 784, 17

Bustamante, M., Heinze, J., Murase, K., & Winter, W. 2017, ApJ, 837, 33

Courant, R., Friedrichs, K., & Lewy, H. 1928, Mathematische Annalen, 100, 32

Daigne, F., & Mochkovitch, R. 1998, MNRAS, 296, 275

De Colle, F., Granot, J., López-Cámara, D., & Ramirez-Ruiz, E. 2012, ApJ, 746, 122

Duffell, P. C., & MacFadyen, A. I. 2013, ApJ, 775, 87

Genet, F., & Granot, J. 2009, MNRAS, 399, 1328

Granot, J. 2005, ApJ, 631, 1022

Granot, J., Cohen-Tanugi, J., & Silva, E. C. 2008, ApJ, 677, 92

Guetta, D., Spada, M., & Waxman, E. 2001, ApJ, 557, 399

Harten, A., Lax, P. D., Leer, B., & v., 1983, SIAM Review, 25, 35

Jacovich, T. E., Beniamini, P., & van der Horst, A. J. 2021, MNRAS, 504, 528

Kaiser, C. R., Sunyaev, R., & Spruit, H. C. 2000, A&A, 356, 975

Kino, M., Mizuta, A., & Yamada, S. 2004, ApJ, 611, 1021

Kobayashi, S., Piran, T., & Sari, R. 1997, ApJ, 490, 92

Kumar, P., & Granot, J. 2003, ApJ, 591, 1075

Malzac, J. 2014, MNRAS, 443, 299

Mathews, W. G. 1971, ApJ, 165, 147

McCarthy, G. A., & Laskar, T. 2024, ApJ, 970, 135

Mignone, A., & Bodo, G. 2005, MNRAS, 364, 126

Mignone, A., & McKinney, J. C. 2007, MNRAS, 378, 1118

Nakar, E., & Piran, T. 2002, MNRAS, 331, 40

Nakar, E., Ando, S., & Sari, R. 2009, ApJ, 703, 675

Oganesyan, G., Nava, L., Ghirlanda, G., & Celotti, A. 2017, ApJ, 846, 137

Pe'er, A. 2014, Space Sci. Rev., 183, 371

Pe'er, A., Long, K., & Casella, P. 2017, ApJ, 846, 54

Pellouin, C., & Daigne, F. 2024, A&A, 690, A281

Preece, R., Burgess, J. M., Von Kienlin, A., et al. 2014, Science, 343, 51

Rahaman, S. M., Granot, J., & Beniamini, P. 2024a, MNRAS, 528, 160

Rahaman, S. M., Granot, J., & Beniamini, P. 2024b, MNRAS, 528, L45

Ravasio, M., Ghirlanda, G., Nava, L., & Ghisellini, G. 2019, A&A, 625, A60

Rees, M. 1978, MNRAS, 184, 61P

Rees, M. J., & Mészáros, P. 1994, ApJ, 430, L93

Rudolph, A., Bošnjak, Ž., Palladino, A., Sadeh, I., & Winter, W. 2022, MNRAS, 511, 5823

Rudolph, A., Petropoulou, M., Bošnjak, Ž., & Winter, W. 2023, ApJ, 950, 28

Sari, R., & Piran, T. 1997, MNRAS, 287, 110

Spada, M., Ghisellini, G., Lazzati, D., & Celotti, A. 2001, MNRAS, 325, 1559

Spada, M., Panaitescu, A., & Meszaros, P. 2000, ApJ, 537, 824

Vianello, G., Gill, R., Granot, J., et al. 2018, ApJ, 864, 163

Yan, Z.-Y., Yang, J., Zhao, X.-H., Meng, Y.-Z., & Zhang, B.-B. 2024, ApJ, 962, 85

Zanotti, O., Rezzolla, L., Del Zanna, L., & Palenzuela, C. 2010, A&A, 523, A8
## Appendix A: Hydrodynamics in planar geometry

In planar geometry, the quantities are constant across the two shocked regions. In particular, the shocked material's proper speed, $u_3 = u_2 \equiv u_0$, can be conveniently determined in the rest frame of shell 1 (see Appendix B of R24a):

$$u_{21} = u_{41} \left( \frac{2f^{3/2}\Gamma_{41} - f(1+f)}{2f(u_{41}^2 + \Gamma_{41}^2) - (1+f^2)} \right)^{1/2}, \quad (A.1)$$

where $u_{41}$ and $\Gamma_{41}$ are the proper velocity and Lorentz factor of shell 4 in the rest frame of shell 1. Eq. (A.1) is transformed to the source frame by velocity transformation via $u_0 = \Gamma_{21}\Gamma_1(\beta_{21} + \beta_1)$. The density and pressure in the shocked shells are determined from the relative Lorentz factor between up- and downstream, $\Gamma_{34}$ for the RS and $\Gamma_{21}$ for the FS, using shock jump conditions. $\Gamma_{21}$ is derived from Eq. (A.1), and $\Gamma_{34}$ is easily obtained by either velocity transformation or using the result $(\Gamma_{21}^2 - 1) = f(\Gamma_{34}^2 - 1)$ that applies to the collision of cold shells (see R24a). The shock front velocities $\beta_{RS}$ and $\beta_{FS}$ can be derived from these relative Lorentz factors (see Appendix C of R24a for the full derivation) and, in turn, we would obtain the crossing times of both shocks:

$$\beta_{RS} = \frac{\beta_4 - 4\Gamma_{34}\left(\frac{u_0}{\Gamma_4}\right)}{1 - 4\Gamma_{34}\left(\frac{\Gamma_0}{\Gamma_4}\right)}, \quad \beta_{FS} = \frac{\frac{1}{4\Gamma_{21}}\left(\frac{u_1}{\Gamma_0}\right) - \beta_2}{\frac{1}{4\Gamma_{21}}\left(\frac{\Gamma_1}{\Gamma_0}\right) - 1}, \quad (A.2)$$

$$t_{RS} = \frac{\beta_4 t_{on,4}}{\beta_4 - \beta_{RS}}, \quad t_{FS} = \frac{\beta_1 t_{on,1}}{\beta_1 - \beta_{FS}}. \quad (A.3)$$

To make sure our numerical methods fit the analytical framework of Rahaman et al. (2024a), we ran the simulation described in Sect. 2.4 up to the point where each of the shocks have crossed their respective shells. Fig. A.1 shows three snapshots of our fiducial simulation: (a) the initial conditions, (b) $t_0$ after the collision, and (c) 1.5 $t_0$ after the collision. The dashed black lines show the expected values using the analytical results from R24a for the same set of parameters. The zones are determined with a combination of the passive tracers, the position of the shocks identified by GAMMA, and an interface detection algorithm based on peak detections in the gradient of the hydrodynamical variables to identify the rarefaction wave front when no shock is present in a given shell. We see a very good agreement for both the hydrodynamical variables and the position of the interfaces of both shocks and the rarefaction wave after the first shock crossing. This interface analysis is completed with a method to recover the shock front velocities by detecting when the front enters a cell for the first time and assuming the front is at the interface between this cell and its already shocked neighbor at the half time step.

The values across the shocked region remain constant in time with very good accuracy compared to the analytical results: the proper speed, $u$, matches its expected value with a precision of the order $10^{-5}$, and derived quantities such as the shock strength $\Gamma_{ud} - 1$ offers a match to within a few times $10^{-4}$. We also recover the analytical values for the shock crossing times $t_{RS}$ and $t_{FS}$ with an accuracy of $10^{-3}$, consistent with a systematic error of a few time steps due to the finite width of the numerical shock.

## Appendix B: Flux calculation in the power-law approximation

### B.1. Analytical flux

We present here the formalism for the analytical calculation of the flux received at observer time $T$ and frequency $\nu$ from a range

of emitting radii. From Eq. (4) we see that a single pulse emitted at a single time and radius in the source frame can be seen over a range of observed times. This means that several source times and radii contribute to a single observer time, defining a surface called equal arrival time surface (EATS) as the locus of all emission points from which photons will arrive at the same time to the observer. In R24b, itself expanding on previous works (e.g., Granot 2005; Granot et al. 2008; Genet & Granot 2009), the authors derive the observed flux coming from a shock front propagating with Lorentz factor $\Gamma_{sh} = \Gamma_{sh,0}(R/R_0)^{-m/2}$, emitting between radii $R_0$ and $R_f = (1 + \Delta R/R_0)R_0$ at a peak frequency $v'_{p} = v'_0(R/R_0)^d$ and with peak luminosity $L'_{v_p} = L'_0(R/R_0)^a$. As quantities introduced in this subsection are relative to a single shock front, we have dropped the subscript RS (respectively FS) for the sake of readability. Additionally in R24, compared to GG09, it is the shocked material traveling with a Lorentz factor of $\Gamma = g\Gamma_{sh}$ that is considered to be the source of emission. This means that while the shock front itself serves to determine the location of the emitting surface, it is the emitting material that will be considered for the Doppler factor.

The flux is calculated against normalized time $\tilde{T}$ defined as:

$$\tilde{T} = \frac{T - T_{\text{ej},0}}{T_0}, \quad (B.1)$$

where $T_{\text{ej},0}$ and $T_0$ are the effective (observed) ejection time and angular time as defined Eq. (5) for the shock front at $R_0$ the initial collision radius:

$$T_0 = (1+z) \left( \frac{1-\beta_{\text{sh},0}}{\beta_{\text{sh},0}} \right) \frac{R_0}{c} \approx (1+z) \frac{R_0}{2c\Gamma_{\text{sh},0}^2}, \quad (B.2)$$

$$T_{\text{ej},0} = (1+z) \left( t_0 - \frac{R_0}{\beta_{\text{sh},0}c} \right). \quad (B.3)$$

We do not provide here the details of the EATS and intermediate variables definitions, and instead refer again to appendix F of R24b for full details of the calculation. We introduce the normalized radius $y(\tilde{T}) \equiv R/R_L(\tilde{T})$ and obtain the flux by integrating from $y_{\text{min}}$ to $y_{\text{max}}$, defined as:

$$y_{\text{min}}(\tilde{T}) = \begin{cases} 1 & \text{if } \tilde{T} \le 1, \\ \tilde{T}^{-1/(m+1)} & \text{if } \tilde{T} \ge 1, \end{cases} \quad (B.4)$$

$$y_{\text{max}}(\tilde{T}) = \begin{cases} 1 & \text{if } \tilde{T}_f \le 1, \\ \tilde{T}_f^{-1/(m+1)} & \text{if } \tilde{T}_f \ge 1, \end{cases} \quad (B.5)$$

where $\tilde{T}_f = (1+\Delta R/R_0)^{m+1}$. The Doppler factor is also expressed as a function of $y$:

$$\delta = (1+z) \frac{\nu}{\nu'} = \frac{1+z}{\Gamma(1-\beta\mu)} \approx (1+z) \frac{2(m+1)g\Gamma_{\text{sh},0}y^{\frac{m}{2}}\tilde{T}^{\frac{-m}{2(m+1)}}}{g^2y^{-1} + (m+1-g^2)y^m}. \quad (B.6)$$

The received flux at normalized observer time $\tilde{T}$ for any given comoving luminosity $L'_{v'}(y)$ can be expressed as:

$$F_{\nu}(\tilde{T}) = \frac{1+z}{8\pi d_L^2} \int_{y_{\text{min}}}^{y_{\text{max}}} \left| \frac{d\mu}{dy} \right| \delta^3 L'_{v'}(y) dy. \quad (B.7)$$

With the change of variable:

$$\frac{d\mu}{dy} = \frac{y^{-2}(1 + my^{m+1})}{2(m+1)\Gamma_{\text{sh},0}^2} \tilde{T}^{m/(m+1)}, \quad (B.8)$$
<figure><img src="image_10.png" /><figcaption>Fig. A.1. Snapshots of the fiducial planar run ($a_u = 2$, $E_{k0,4} = E_{k0,1}$, $t_{on,1}$: $t_{on,4}$: $t_{off} = 1:1:1$.) at $t = t_0$ (left), before first shock crossing (center), and after shock crossing (right). Shells are color-coded according to the identified zones, the dashed lines show theoretical values. $t_0 = 2.7 \times 10^3$ s and $R_0 = 8 \times 10^{13}$ cm for this set of parameters.</figcaption></figure>

and using Eq. (B.6) in B.7:

$$F_v(\tilde{T}) = \frac{(1+z)g^3\Gamma_{\text{sh},0}\tilde{T}^{\frac{-m}{2(m+1)}}}{2\pi d_L^2} \int_{y_{\text{min}}}^{y_{\text{max}}} dy L'_{v'}(y)y^{-1-\frac{m}{2}} \times \left[1+g^2\frac{y^{-(m+1)}-1}{m+1}\right]^{-2} \left(\frac{1+my^{m+1}}{g^2+(m+1-g^2)y^{m+1}}\right). \quad (B.9)$$

In the general case, we chose power-law scalings for the luminosity and peak frequency with radius:

$$L'_{v'} = L'_0 \left(\frac{R}{R_0}\right)^a S \left(\frac{v'}{v'_p}\right) = L'_0 \tilde{T}^a y^a S \left(\frac{v'}{v'_p}\right), \quad (B.10)$$

$$\frac{v'}{v'_p} = \tilde{T}^{-d} y^{m/2-d} \left[ 1 + g^2 \frac{(y^{-(m+1)} - 1)}{m+1} \right] \frac{v}{v_0}. \quad (B.11)$$

## B.2. Numerical flux normalization

To obtain the appropriate flux from a numerical cell corresponding to the infinitely thin shell approximation used in R24b, we compare the isotropic energy obtained by integrating Eq. (4) with the equivalent integration of the flux formula from De Colle et al. (2012) for the contribution from a single computational cell, which for an isotropic comoving emissivity per unit volume and frequency, $P'_{v',jk}$, reads:

$$\begin{gathered} \Delta F_{v,ijk} = \frac{(1+z)^2}{4\pi d_L^2} \frac{\Delta V_{ijk}^{(4)}}{\Delta T_{\text{obs,i}} \Gamma_{jk}^2 (1-\beta_{jk} \cos \theta_{jk})^2}, \\ \text{for } \left| T_{\text{obs,z,i}} - t_{jk} + \frac{r_{jk} \cos \theta_{jk}}{c} \right| < \frac{\Delta T_{\text{obs,z,i}}}{2}. \quad (B.12) \end{gathered}$$

Here $T_{\text{obs,z}} = T_{\text{obs}}/(1+z)$ and $v_z = (1+z)v$ are measured in the source's cosmological frame, as is the (Lorentz invariant) 4-volume $\Delta V^{(4)}$. In this section, subscript i refer to observer time bins over which flux is calculated, subscript j to the simulation time step, and subscript k to the numerical cell within data file j. This approach attributes all of the contribution from any given numerical 4D cell jk to a single observer time interval i. This approximation relies on sufficiently fine spatial (and temporal) resolution and on the Doppler factor not significantly varying

within the cell, due to the variation of its velocity direction with respect to the direction to the observer. However, our simulation is spatially 1D: the Doppler factor significantly varies within each numerical cell of lab-frame volume $\Delta V_{ijk}^{(4)} = 4\pi r^2 \Delta r$ and contributes to different observed time bins. We thus divide 1D cells further into subcells along the $\theta$ coordinate in which the Doppler factor can be considered constant.

The observed flux depends only on $\mu = \cos\theta$ through the Doppler factor and the photon arrival time but not on the $\phi$ coordinate as we calculate flux along the outflow axis. Thus each subdivision of volume $\Delta V_{ijk}^{(4)} = 2\pi r_{jk}^2 \Delta r_{jk} \Delta \mu_i \Delta t_j$ contributes to the time bin centered on $T_{\text{obs,i}} = (1+z)(t_j - r_{jk}\mu_i/c)$ of width $\Delta T_{\text{obs,i}} = (1+z)(r_{jk}/c)\Delta \mu_i$. Recalling $E_{\text{iso}} = \int dv_z \int dT_{\text{obs,z}} L_{v_z,\text{iso}} = (4\pi d_L^2/(1+z)) \int dT_{\text{obs}} \int dV F_v(T_{\text{obs}})$, the isotropic energy for a spherical cell reads:

$$\begin{aligned} E_{\text{iso,jk}}^{(\text{num})} &= \frac{4\pi d_L^2}{1+z} \int dv \sum_i \Delta T_{\text{obs,i}} \Delta F_{v,ijk} \\ &= \int dv_z \sum_i \Delta V_{ijk}^{(4)} \frac{P'_{v',jk}}{\Gamma_{jk}^2 (1 - \beta_{jk} \mu_i)^2} \\ &= 2\pi r_{jk}^2 \Delta r_{jk} \Delta t_j \frac{P'_{v',m,jk}}{\Gamma_{jk}^2} \int_0^\infty dv_z \int_{-1}^1 d\mu (1 - \beta\mu)^{-2} S \left( \frac{v'}{v'_{m,jk}} \right) \\ &= \Delta V_{jk}^{(4)} (1 - \beta_{jk}) \Gamma_{jk}^2 \nu_{m,z,jk} P'_{v',m,jk} \quad \mathcal{A} \cong \frac{1}{2} \Delta V_{jk}^{(4)} \nu_{m,z,jk} P'_{v',m,jk} \quad \mathcal{A}, \end{aligned} \quad (B.13)$$

having changed variables to $x \equiv v'/v'_m = v(1-\beta\mu)/[\nu_m(1-\beta)] = v_z(1-\beta\mu)/[\nu_{m,z}(1-\beta)]$ and denoting $\mathcal{A} = \int S(x)dx$. In De Colle et al. (2012) the flux at a fixed observed time bin was calculated by summing over j and k. To obtain the isotropic energy we performed the opposite, namely, fixing j and k and summing over i (and integrated over frequency). Performing the equivalent integration on a single shell in the infinitely thin shell approximation (Genet & Granot 2009) is expressed as:

$$\begin{aligned} E_{iso}^{(th)} &= \frac{4\pi d_L^2}{1+z} \int dT_{\text{obs}} \int dv F_v \\ &= \bar{L}_{v_m} T_{\theta} \int \frac{d\tilde{T}}{\tilde{T}^2} \int dv S \left(\frac{v}{v_m}\tilde{T}\right) = \frac{1}{2}\nu_{m,z}\bar{L}_{v_{m,z}} T_{\theta,z}\mathcal{A}, \end{aligned} \quad (B.14)$$
using change of variable $x \equiv (\nu/\nu_m)\tilde{T}$. Finally, we obtain the equivalent numerical luminosity for simulation cell k of data file j in the infinitely thin shell approximation:

$$ \tilde{L}_{\nu_{m,z},jk} = \frac{\Delta V_{jk}^{(4)} P'_{\nu_{m}jk}}{T_{\theta,z,jk}} = \frac{\Delta E'_{\nu_{m}jk}}{T_{\theta,z,jk}}. \qquad (B.15) $$

Here $\Delta E'_{\nu_{m}jk}$ is the energy emitted per unit frequency in the comoving frame by the four-dimensional cell jk. Going further with the infinitely thin shell approximation, we can assume all the energy given to the accelerated electrons is radiated over a time smaller than the numerical time step, and we write:

$$ \tilde{L}_{\nu_{m,z},jk} = \frac{\Gamma_{jk} \Delta V_{jk}^{(3)} \epsilon_e e'_{jk}}{W(p) \nu'_{m,jk} T_{\theta,z,jk}}. \qquad (B.16) $$

### B.3. Flux calibration

While the general analytical formula for any scaling indices $a, d, m$ was not obtained and the integral in Eq. (B.9) must be calculated numerically, two special cases with a simpler integral form are of interest to us. We use the first case $a, d, m = 0, 0, 0$ to calibrate the post-processing methods in a consistent way from our fully planar simulation before extending them to the physically-relevant spherical geometry. With these simplifying assumptions, the observed flux is:

$$ F_{\nu}(\tilde{T} \ge 1) = \frac{(1+z)g^3\Gamma_{\text{sh}}L_0'}{2\pi d_L^2} \int_{y_{\text{min}}}^{y_{\text{max}}} \frac{dy \, y^{-2} S(x)}{[1+g^2(y^{-1}-1)]^3}. \qquad (B.17) $$

With

$$ x = \frac{\nu'}{\nu_0'} = [1 + g^2(y^{-1} - 1)] \tilde{v}, \qquad (B.18) $$

$$ \tilde{v} = \nu/\nu_0, \qquad (B.19) $$

$$ \tilde{T} = \frac{T - T_{\text{ej}}(R_0, \beta_{\text{sh}})}{T_{\theta}(R_0, \Gamma_{\text{sh}})} \equiv 1 + \tilde{T}, \qquad (B.20) $$

and the integration limits $y_{\text{min}}, y_{\text{max}}$ are defined in Eqs. (B.4)-(B.5). Choosing $x$ as the new integration variable, we obtain:

$$ F_{\nu} = \frac{(1+z)2\Gamma_0 L_0'}{4\pi d_L^2} \left(\frac{\nu}{\nu_0}\right)^2 \int_{x_1}^{x_u} dx \, x^{-3} S(x). \qquad (B.21) $$

The upper and lower integration limits are, respectively, $x_u = x(y_{\text{min}})$ and $x_l = x(y_{\text{max}})$ and the integral is performed numerically using the quad method from scipy odeint<fnref n="2" />. We present in Fig. B.1 light curves at $\tilde{v} = 10^{-1}$ and instantaneous spectra at $\tilde{T} = 1$ obtained from our fiducial run in planar geometry, compared to the analytical flux obtained analytically with Eq. (B.21). The figure presents the observed flux using either the smooth Band function (top row) or a broken power-law (bottom row). Depending on the local spectral shape, the relative contributions of both shock to the light curve between their respective peak frequencies may vary, also causing a change in the peak flux at this intermediary frequency range, but the overall shape of the light curve obtained by the sum of both contributions do not. This difference between peak frequencies is better shown in the instantaneous spectrum where the shape of the curve between the peaks shows more variation with spectral shape. Our numerical fluxes agree with their analytical counterparts to a good degree of precision.

<figure><img src="image_11.png" /><figcaption>Fig. B.1. Light curves (left column) and instantaneous spectra (right column) obtained from the fiducial run in planar geometry (full lines) compared to the analytical expectations (dash-dotted lines). Normalizations of $\tilde{T}$ and $\tilde{v}$ are defined at the RS. The top row displays the result for a Band spectral shape, and in the bottom the synchrotron broken power-law. The respective contributions of each shock to the total observed flux (in black) are given in red for the RS and blue for the FS.</figcaption></figure>

## Appendix C: The effective radius approximation

<figure><img src="image_12.png" /><figcaption>Fig. C.1. Schematic representation of the effective radius approximation used to derive peak frequency and flux from the hydrodynamical data depending on the asymptotical regime. The EATS is shown in red.</figcaption></figure>

<footnote n="2">² An analytic solution of the integral exist for both spectral shapes, but in the case of the Band function this form features the incomplete Gamma function, which is not implemented in Python for negative values of the first parameter. See Appendix F in R24b.</footnote>