© The Authors 2024

1. Introduction

The early steep decay (ESD) is a common feature in gamma-ray burst (GRB) light curves. It is observed in X-rays during the transition between the GRB prompt phase and the subsequent GRB afterglow phase. At the end of the violent prompt phase, the observed luminosity in X-rays drops rapidly by several orders of magnitude with a temporal index of −3 to −5. Typical durations for the ESD are ∼10² − 10³ s (Nousek 2006). This behaviour seems to be a robust feature in GRBs: if an X-ray telescope – such as the X-Ray Telescope (XRT) on board Swift – manages to observe a GRB early enough, this peculiarity is almost always seen.

A natural explanation is that the ESD is the consequence of high-latitude emission (HLE) from the last flashing shell in the optically thin regime (Fenimore et al. 1996; Kumar & Panaitescu 2000; Genet & Granot 2009). In this case, the bolometric luminosity, L_bol, and the peak energy of the spectrum, E_p, are expected to decrease with time approximately as L_bol ∝ t⁻³ and E_p ∝ t⁻¹. This theoretical expectation was recently observationally corroborated by Tak et al. (2023; see also Uhm et al. 2022). By performing time-resolved fits of broad GRB prompt pulses, the authors found a parameter evolution consistent with that expected from HLE in a majority of the pulses examined. Given that HLE is observed in the decay phases of prompt pulses, as found by Tak et al. (2023), one could argue that it is plausibly the origin for the ESD as well, since there is often a smooth transition between the last prompt pulse and the ESD.

By contrast, Ronchini et al. (2021) performed time-resolved analysis during the ESD using data from XRT and found that the spectral evolution does not match that predicted for HLE. During the observations, the peak energy seems to cross the whole XRT band, indicating a stronger spectral evolution as E_p ∝ t⁻² to t^−2.5. By fitting a power-law spectrum to the XRT data, they found a relation between the fitted spectral index, α(t), and the ratio of the maximum flux at the onset of the decay to the flux at time t, F_max/F(t). Ronchini et al. (2021) referred to this correlation as the α − F relation, which they interpreted as being due to adiabatic cooling of the emitting particles in the context of a proton synchrotron origin for the GRB emission.

In photospheric models of GRBs, the prompt radiation is emitted when the initially opaque jet transitions to the optically thin regime (Paczyński 1986; Goodman 1986). This transition usually occurs once the acceleration of the outflow is complete, at a characteristic radius, R_ph, given by (e.g. Hascoët et al. 2013)

$$\begin{array}{c}\hfill R_{\mathrm{ph}}=\frac{f{\kappa }_{\mathrm{T}}\dot{E}}{8\pi c^3(1+\sigma ){\mathrm{\Gamma }}^3}=5.9\times {10}^{12}\frac{f{\dot{E}}_{52}}{(1+\sigma ){\mathrm{\Gamma }}_2^3}\mathrm{cm}.\end{array}$$(1)

Here, Ė is the isotropic equivalent total injected power in the outflow, Γ is the bulk Lorentz factor, σ is the magnetization at large distances, f is the number of leptons per baryon, and κ_T = 0.4 cm² g⁻¹ is the Thomson opacity. In the following, we assume a negligible magnetization at large distances (σ ≪ 1). The total power, Ė, is equal to the kinetic power, Ė_K, in this case. The corresponding geometrical timescale is given by

$$\begin{array}{c}\hfill {\tau }_{\mathrm{geo}}=\frac{R_{\mathrm{ph}}}{2c{\mathrm{\Gamma }}^2}=\frac{f{\kappa }_{\mathrm{T}}{\dot{E}}_{\mathrm{K}}}{16\pi c^4{\mathrm{\Gamma }}^5}=9.8\times {10}^{-3}\frac{f{\dot{E}}_{\mathrm{K},52}}{{\mathrm{\Gamma }}_2^5}\mathrm{s},\end{array}$$(2)

which is very short unless Γ is small or Ė_K or f is huge (see e.g. Dereli-Bégué et al. 2022; Samuelsson & Ryde 2023).

As discussed by Hascoët et al. (2012), a short geometrical timescale would not be in line with a HLE interpretation for the ESD in photospheric models (or the decay observed at the end of prompt pulses), since the total duration of the ESD should be of the order of ∼τ_geo (however, see Pe’er et al. 2006 for a longer-lasting ESD originating from photons diffusing in the surrounding cocoon). Furthermore, τ_geo is very sensitive to Γ as evident from Eq. (2), while the ESD appears to be a robust feature, nearly always present at the start of the early afterglow of GRBs with a generic temporal decay index of −3 to −5 (Nousek 2006).

In light of the above argumentation, we investigate in this paper an alternative origin for the ESD in photospheric models, namely that it results from the intrinsic evolution of the dying central engine (as suggested by Hascoët et al. 2012). Depending on the last stages in the life of the central source, the photospheric emission may indeed generate an ESD mimicking that expected from HLE or a stronger spectral evolution, such as the one found by Ronchini et al. (2021). Assuming a power-law decline in the emitted power and Lorentz factor of the flow, we estimate the observed light curve and spectrum under different scenarios and compare the results to observations.

The paper is structured as follows. In Sect. 2 we introduce our model for the last stages of central activity. To predict the observed signal, we study two benchmark scenarios: a non-dissipative model in Sect. 2.1 and a dissipative model in Sect. 2.2. In Sect. 3 we present our results and discuss them, with a specific emphasis on parameter dependence and underlying assumptions, in Sect. 4. We conclude in Sect. 5. We employ the notation Q_x = Q/10^x throughout the text.

2. The ESD in photospheric models

We modelled the last stages of the central engine with power-law decays for the injected power and Lorentz factor of the relativistic wind as

$$\begin{array}{c}\hfill {\dot{E}}_{\mathrm{K}}={\dot{E}}_{\mathrm{b}}{\left(\frac{t}{t_{\mathrm{b}}}\right)}^{-\lambda },\mathrm{\Gamma }={\mathrm{\Gamma }}_{\mathrm{b}}{\left(\frac{t}{t_{\mathrm{b}}}\right)}^{-\gamma },t>t_{\mathrm{b}},\end{array}$$(3)

where t is the time in the source rest frame, t_b is the break time when the ESD starts, and Ė_b and Γ_b are, respectively, the kinetic power and Lorentz factor at t_b. The choice of power-law decays for Ė_K and Γ is discussed in Sect. 4.2.

If the emitted bolometric luminosity follows the evolution of Ė_K given in the equation above – that is, if the radiative efficiency does not vary too much along the ESD – then one expects λ ∼ 3 to account for ESD light curves. If, in addition, f does not vary too much during the ESD either, one finds from Eq. (1) that the photospheric radius evolves as R_ph ∝ t^{3(γ − 1)}. Thus, the photosphere can either shrink or inflate with time depending on the value of (γ − 1).

2.1. Non-dissipative model

In a non-dissipative photospheric model, the radiation below the photosphere is at all times kept in thermodynamic equilibrium with the plasma. Under this assumption, the temperature of the observed radiation and the observed isotropic equivalent bolometric luminosity at the photosphere are given by (e.g. Piran 1999; Daigne & Mochkovitch 2002; Pe’er et al. 2007)

$$\begin{array}{c}\hfill T_{\mathrm{obs}}=T_0{\left(\frac{R_{\mathrm{ph}}}{R_{\mathrm{s}}}\right)}^{-2/3},L_{\mathrm{bol}}={\dot{E}}_{\mathrm{th}}{\left(\frac{R_{\mathrm{ph}}}{R_{\mathrm{s}}}\right)}^{-2/3}.\end{array}$$(4)

Here, T₀, Ė_th, and R_s = ΓR₀ are, respectively, the initial temperature, the injected power in thermal form at the base of the jet, and the saturation radius, and R₀ is the distance from the central engine at the base of the jet. The initial temperature is given by

$$\begin{array}{c}\hfill k{T}_0\approx k{\left(\frac{{\dot{E}}_{\mathrm{th}}}{4\pi R_0^{2}ac}\right)}^{1/4}=1.2{\dot{E}}_{\mathrm{th},52}^{1/4}{R}_{0,7}^{-1/2}\mathrm{MeV}.\end{array}$$(5)

We assumed that thermal energy is efficiently converted into kinetic energy below the photosphere, implying that Ė_th ≈ Ė_K.

The efficiency of the non-dissipative model can be evaluated from Eqs. (1) and (4):

$$\begin{array}{c}\hfill {ϵ}_{\gamma }=\frac{L_{\mathrm{bol}}}{{\dot{E}}_{\mathrm{th}}}={\left(\frac{R_{\mathrm{ph}}}{R_{\mathrm{s}}}\right)}^{-2/3}=3\times {10}^{-3}{\left(\frac{f{\dot{E}}_{\mathrm{K},52}}{R_{0,7}{\mathrm{\Gamma }}_2^4}\right)}^{-2/3},\end{array}$$(6)

which is quite small for standard values of the parameters. Similarly, the observed temperature, equal to T₀ϵ_γ, is expected in the keV range.

Due to HLE and contributions from different optical depths to the released emission, the observed spectrum consists of a superposition of black bodies at different temperatures. This leads to a modified, broadened spectrum compared to a Planck function in the observer frame. This spectrum can be approximated by a cutoff power-law function with a spectral index α = 0.4 and a peak energy of $E_{\mathrm{p}}^{\mathrm{th}}=3.9k{T}_{\mathrm{obs}}$ (Beloborodov 2010). Assuming λ = 3 and that R₀ and f do not vary during the ESD, a HLE-like decline of the peak energy with time as t⁻¹ requires γ = 27/32 ∼ 0.8.

2.2. Dissipative model

Energy injection below the photosphere provides a way to increase the efficiency, raise the peak energy of the emitted spectrum, and transform its shape via Comptonization. In this section, we considered an unspecified dissipation mechanism that continuously injects energy into the photon distribution. Such a mechanism could be, for instance, magnetic dissipation (Drenkhahn & Spruit 2002; Giannios 2008), long-lasting dissipation via turbulence or multiple shocks (Rees & Mészáros 2005; Zrake et al. 2019), or collisional heating between neutrons and protons (Beloborodov 2010). We chose not to detail the processes involved and rather adopted some very simplified assumptions to estimate the spectral evolution and the outflow parameters during the ESD.

The bolometric luminosity was obtained by fixing the efficiency such that

$$\begin{array}{c}\hfill L_{\mathrm{bol}}={ϵ}_{\gamma }{\dot{E}}_{\mathrm{K}},\end{array}$$(7)

with ϵ_γ = 0.1 − 0.5. The efficiency can be large if dissipation takes place close to the photospheric radius (Beloborodov 2010; Gottlieb et al. 2019; Samuelsson & Ryde 2023). Although the efficiency could very well vary during the ESD, here we kept it constant for simplicity. This is somewhat in line with the numerical simulation presented in Gottlieb et al. (2019), where ϵ_γ is found to fluctuate around a central value of ϵ_γ ∼ 0.5. In this paper, we focus on the global behaviour and neglect possible small-scale variations in Ė_K, Γ, and/or ϵ_γ, which we argue can be at least partially suppressed (see Sect. 4.3). If there exists a global time dependence for ϵ_γ, the results presented in Sect. 3 can still be obtained by adjusting λ and γ accordingly.

If the bolometric radiation came from a blackbody¹, the peak energy would be related to the effective temperature, T_eff, as (Beloborodov 2013)

$$\begin{array}{cc}& E_{\mathrm{p}}^{\mathrm{th}}\approx 4\mathrm{\Gamma }k{T}_{\mathrm{eff}}=44L_{\mathrm{bol},52}^{1/4}{\left(\frac{{\mathrm{\Gamma }}_2}{R_{\mathrm{ph},13}}\right)}^{1/2}\mathrm{keV}\hfill \\ \hfill & =57L_{\mathrm{bol},52}^{-1/4}{\mathrm{\Gamma }}_2^2{\left(\frac{{ϵ}_{\gamma }}{f}\right)}^{1/2}\mathrm{keV},\hfill \end{array}$$(8)

where the third equality employed Eq. (7). We assumed that the Comptonized spectrum satisfies $E_{\mathrm{p}}\sim \mathrm{a}\mathrm{few}{E}_{\mathrm{p}}^{\mathrm{th}}$, as argued in Beloborodov (2013). If, in addition, one adopts the simplifying assumptions that ϵ_γ, f, and the ratio $E_{\mathrm{p}}/E_{\mathrm{p}}^{\mathrm{th}}$ do not vary too much during the ESD, one finds from Eq. (8) that the peak energy follows a power law:

$$\begin{array}{c}\hfill E_{\mathrm{p}}\propto {\dot{E}}_{\mathrm{K}}^{-1/4}{\mathrm{\Gamma }}^2\propto t^{\lambda /4-2\gamma }.\end{array}$$(9)

Adopting λ = 3, a HLE-like behaviour (i.e. E_p ∝ t⁻¹) is obtained for γ = 7/8, while the strongest spectral evolution found by Ronchini et al. (2021) requires γ = 1.4–1.6.

Lastly, we assumed that the dissipation transforms the observed spectrum to be similar to that of typical GRB observations. Specifically, we considered the observed spectrum to be a Band function (Band et al. 1993), with a soft, low-energy power-law index α = −1 and a high-energy power-law index β = −2.3. We discuss the influence of this choice on the results in Sect. 4.1.

3. Results

3.1. Spectral evolution during the ESD

In Fig. 1, the evolution of the spectral shape during the ESD is shown, in the non-dissipative (left panel) and in the dissipative scenario (right panel). Both scenarios are shown for γ = 7/8 (dashed lines, intended to mimic a HLE-like evolution), and γ = 1.5 (solid lines). As expected, when γ = 1.5, the decrease in E_peak is much faster compared to when γ = 7/8.

Fig. 1. Snapshot spectra showing the evolution of the spectral luminosity during the ESD in the non-dissipative scenario (left) and the dissipative scenario (right). Dashed lines have γ = 7/8, while solid lines have γ = 1.5, as also indicated in the figure. The spectra are plotted at even intervals in time from the onset of the ESD at t/tb = 1 (dark colour) until t/tb = 10 (bright colour). Note that the two spectra for γ = 7/8 and γ = 1.5 overlap at t/tb = 1. The purple shading indicates the energy sensitivity of XRT. Parameter values used are given in Table 1.

In Fig. 2 we show the obtained α − F relation. The maximum flux, F_max, is calculated as the integrated spectral flux within the XRT energy band (0.5 − 10 keV) at the onset of the ESD (first spectrum in Fig. 1), and F(t) is the flux in the XRT band at time t. The spectral index is estimated as in Ronchini et al. (2021):

$$\begin{array}{c}\hfill {\alpha }_{\mathrm{XRT}}=-\frac{log[N_{\nu =10\mathrm{keV}/h}(t)/N_{\nu =0.5\mathrm{keV}/h}(t)]}{log(10\mathrm{keV}/0.5\mathrm{keV})},\end{array}$$(10)

Fig. 2. Photon index as a function of Fmax/F along the ESD branch. The grey points are the best-fit values from time-resolved analyses of eight GRBs, as obtained in Ronchini et al. (2021). Parameter values used are given in Table 1.

where N_ν = E/h is the spectral number density, evaluated at energy E.

The line-coding in Fig. 2 is the same as for the corresponding scenarios in Fig. 1. The grey points show the best-fit values obtained in Ronchini et al. (2021) in their time-resolved spectral fits during the ESD of eight GRBs. It is clear from the figure that unless there is some dissipation occurring below the photosphere, the spectrum is initially too hard and at late times too soft to account for the observations. The Ronchini et al. (2021) results are well reproduced in the dissipative case with γ = 1.5 (full red line). The peak energy then crosses the entire XRT band, starting at E_p ≳ 100 keV and reaching 0.5 keV, the lower limit of the XRT energy band considered here, at F_max/F ∼ 150 (t/t_b ≳ 10; see Fig. 1). Conversely, γ = 7/8, which mimics the spectral evolution during HLE, fails (by far) to reproduce the Ronchini et al. (2021) results (dashed green line). It should however be noted that the Ronchini et al. (2021) results rely on the delicate correction for absorption, which becomes quite large below 1 keV (see the discussion in Sect. 4.4).

3.2. Outflow parameters during the ESD

The evolution of the Lorentz factor Γ, peak energy E_p, and photospheric radius R_ph in the dissipative scenario are shown in Fig. 3. The quantities are plotted against observer time:

$$\begin{array}{c}\hfill t_{\mathrm{obs}}=(1+z)[t+{\tau }_{\mathrm{geo}}],\end{array}$$(11)

Fig. 3. Evolution of the Lorentz factor, peak energy, photospheric radius, and normalized XRT-flux during the ESD for the dissipative model. The XRT flux in the bottom-right panel is shown for the low-energy index, α, equal to −2/3 (dark colour), −1 (intermediate colour), and −1.5 (bright colour). Parameter values are given in Table 1.

where τ_geo appears to account for the propagation time of the plasma to the photosphere: t_prop = R_ph/2cΓ² = τ_geo, assuming emission on the line of sight. For γ = 7/8, τ_geo is always negligible and Γ, E_p and R_ph just follow power laws of slopes −7/8, −1, and −3/8, respectively. For γ = 1.5, the quantities start decaying as power laws of slopes −1.5, −2.25, and +1.5, from which they deviate once τ_geo becomes comparable to t.

One can estimate when this transition occurs by setting t = τ_geo. This gives

$$\begin{array}{c}\hfill t=t_{\mathrm{HLE}}\equiv {\left(t_{\mathrm{b}}^{5\gamma -\lambda }\frac{100{\mathrm{\Gamma }}_{\mathrm{b},2}^5}{f{\dot{E}}_{\mathrm{b},52}}\right)}^{1/(5\gamma -\lambda -1)},\end{array}$$(12)

with the above equation valid only when 5γ − λ > 1. Adopting λ = 3, γ = 1.5, f = 1, and t_b = 20 s, we obtain t_HLE = 175 s, which, with z = 1, gives t_obs = 700 s. This is the time when HLE starts to dominate. From the figure, it is clear that the deviations become noticeable earlier than this once t_obs ∼ 200 s, which corresponds to τ_geo/t ∼ 0.1.

When t > t_HLE, HLE can become important and our results should be evaluated with caution since we only considered emission on the line of sight. HLE from the optically thick fireball results in the bolometric flux decreasing as t⁻² with a slowly evolving peak-energy (Pe’er & Ryde 2011). We note that when t > t_HLE, the observer would still see HLE from the last emitted regions even if the central engine activity ceased.

Finally, the ESD light curve (in the XRT spectral band 0.5 − 10 keV) is shown in the bottom-right panel of Fig. 3. We represent the light curves for three values of the low-energy spectral index: α = −2/3, −1, and −1.5. When γ = 7/8, the flux immediately starts to decrease rapidly after t_b due to the drop in L_bol. For γ = 1.5 however, the decay is more gradual. This is because the drop in luminosity is partially prevented by the peak energy moving into the XRT band, as also evident from Fig. 1. Both these types of behaviours have been observed. For instance, GRB 090618 and GRB 230420A show very steep decays right after the prompt emission, while GRB 081221 and GRB 210305A show smoother transitions to the ESD phase, similar to the solid black line in Fig. 3².

4. Discussion

4.1. Parameter dependence

There are two observables that we can compare our predictions against during the ESD: the typical XRT light curve and the spectral evolution. The XRT light curve commonly drops ∼2 order of magnitude in ∼100 s, seemingly as a power law. As evident from the bottom-right panel in Fig. 3, the flux decrease predicted by our model depends on the value of γ as well as on the spectral shape. However, it depends most strongly on the choice of λ; a larger λ leads to a more rapid decrease. The drop in flux can be circumvented to an extent by a different choice of γ and/or a different spectral shape, but this requires fine-tuning once λ > 5. The same is true for λ < 3. Thus, we deduce that the isotropic equivalent of the kinetic power must drop as a power law in time with index ∼3 − 5 to account for the observed XRT light curve. If this scenario is correct, this provides an important constraint on the behaviour of the central engine at the end of the relativistic jet launching (see the discussion in Sect. 4.2).

The results regarding the spectral evolution presented in Fig. 2 are most sensitive to the assumed shape of the emitted spectrum. The general trend in the α − F plane (Fig. 2) is that at early times, when F_max/F ≳ 1, as long as E_p is above the XRT band, then α_XRT ∼ |α|. In the case of γ = 1.5 in the dissipative model, E_p rapidly crosses the entire XRT band, which means that α_XRT ∼ |β| at late times. How quickly this transition occurs depends on the initial peak energy: if the peak energy is high, one probes the low-energy part for longer, and the transition occurs later.

According to the Ronchini et al. (2021) analysis, α_XRT seems to saturate around ∼2.3 once F_max/F ≳ 100. If E_p ≪ 0.5 keV at late times, this requires the existence of a persistent high-energy power law with slope β ∼ −2.3. This feature is non-trivial to maintain in a photospheric framework, since if dissipation halts deep below the photosphere, the high-energy photons lose their energy due to Compton scattering. Inelastic scatterings between neutrons and protons may provide such a signature; however, this requires a highly relativistic jet and may, therefore, not be efficient at late times in the current framework (Beloborodov 2010). We note that the α − F relation can be obtained if γ is time-dependent in such a way that the decrease in E_p halts just as α_XRT ∼ 2.3. However, unless one can find a physical motivation for why this would occur, this scenario seems unlikely. Therefore, we conclude that dissipation should continue all the way to the photosphere to accommodate the late-time spectrum (however, see the discussion in Sect. 4.4).

If α_XRT ∼ |α| at early times and α_XRT ∼ |β| at late times, then the dispersion observed in α_XRT should reflect the dispersions found for α and β in GRB catalogues. Specifically, the spread in α_XRT at late times should increase, since the values obtained for β in time-resolved analyses of the prompt emission ranges from −4 to −2 (Poolakkil et al. 2021). Thus, the perceived saturation around ∼2.3 is likely due to the small number of GRBs in the current sample. If the behaviour persists as more GRBs are added to the sample, an alternative explanation for the power-law index at late times may be needed.

4.2. Insights into the progenitor systems

To account for the observations using the model presented herein, a power-law decay of the kinetic power with an index ∼3 − 5 is required. In this section, we briefly discuss what this implies for the dying central engine.

Many detailed numerical simulations of compact binary mergers have been published in the wake of GW170817. Although these numerical simulations cannot yet probe the long timescales discussed here (several 100 s), they can still provide valuable insights. One such insight is that the mass accretion rate onto the black hole after collapse seems to decrease as ∼t⁻² (Christie et al. 2019; Metzger & Fernández 2021; Hayashi et al. 2023; Gottlieb et al. 2023). Accounting for a conversion efficiency between the accretion rate and the jet kinetic power, as well as the uncertainty in ϵ_γ, a decay of the kinetic power as Ė_K ∝ t⁻³ seems possible. This also motivates our choice of a power-law decay for the jet kinetic power in Eq. (3).

For collapsars, the picture is less clear. The mass accretion rate depends on the surrounding stellar density profile, ρ(r)∝r^−a. Assuming free-fall accretion, the black hole mass accretion rate goes as ∼t^{1 − 2a/3}, with the stellar material at an initial radius R reaching the central black hole of mass M_BH at a characteristic time $t=R^{3/2}/\sqrt{2G{M}_{\mathrm{BH}}}$ (e.g. Gottlieb et al. 2022). Assuming perfect conversion between accreted mass rate and jet kinetic power, a power-law decrease of Ė_K ∝ t⁻³ after t = 20 s requires a ∼ 6 at R > 7.5 × 10⁹ cm.

A density decrease with a ∼ 6 is much steeper than what is expected within a stellar envelope (e.g. Woosley & Heger 2006). Indeed, it was recently found that the density profile after core collapse is quite shallow in the inner regions with a ∼ 1.5 (Halevi et al. 2023). However, in the same work, the density was found to decrease much more rapidly close to the stellar edge, at R ≳ 3 × 10⁹ cm. Such a combined stellar density distribution would lead to a jet power that is initially constant, followed by a strong decay. Furthermore, it is sufficient that the diffuse material extends up until ≳3 × 10¹⁰ cm. This would generate an accretion rate ∝t⁻³ over a timescale of ∼200 s in the central engine frame, which is further stretched by a factor of (1 + z) in the observer frame. Thus, we conclude that the envisioned scenario predicts a diffuse density profile with a ∼ 6 near the stellar edge in collapsars, which should extend to a few times 10¹⁰ cm. Since the ESD is a robust feature, these progenitor properties should be quasi-universal, which, therefore, constitutes a powerful test for the model.

The kinetic power is related to the Lorentz factor as Ė_K = Γṁc², where ṁ is the observed mass ejection rate in the jet. If the comoving density in the jet is constant, then ṁ ∝ Γ and we naturally obtain Γ ∝ t^−1.5. In the second considered case, where Γ ∝ t^−7/8, we required a comoving jet density that decreases with time. This could possibly be due to a cleared jet funnel resulting in less mixing at later times.

4.3. Photospheric blending

The ESD is often smooth in time, even if the earlier prompt emission has been highly irregular and chaotic. This may seem to contradict the scenario discussed in this paper: if the ESD probes the dying central engine, surely some variability is expected. However, a variable central engine will inevitably generate emission periods with higher optical depths, leading to larger photospheric radii. These dense shells are going to re-trap emission from inner layers that may already be optically thin. The radiation from these different layers will blend and be emitted together at a later time. This effect ensures that the central engine variability is smoothed out in the observer frame.

Imagine a dense shell emitted at a time t₁. A second shell emitted at a later time t₂ ≡ t₁ + δt will be affected by photospheric blending if the radiation released from the photosphere of the second shell reaches the first shell while it is still optically thick. This gives a condition on δt as

$$\begin{array}{c}\hfill \delta t<\frac{R_{\mathrm{ph},1}}{{\beta }_{1}c}(1-{\beta }_1)-\frac{R_{\mathrm{ph},2}}{{\beta }_{2}c}(1-{\beta }_2)\approx {\tau }_{\mathrm{geo},1}-{\tau }_{\mathrm{geo},2}.\end{array}$$(13)

Here, R_ph, i, β_i, and τ_geo, i are, respectively, the photospheric radius, velocity, and geometrical timescale of the shell emitted at t_i, and the last approximation holds if Γ₁, Γ₂ ≫ 1. If the central engine is variable at the end of its life, the ESD would deviate from a single power law in time. However, small-scale variations of timescales δt would be suppressed in the light curve.

Rapid variability is not observed during the ESD but X-ray flares, with typical timescales of δt/t_obs ∼ 0.1 to 1, are quite common (Burrows et al. 2005; Nousek 2006). Photospheric blending could be an interesting avenue to explore with regards to X-ray flares, since δt ∼ τ_geo. In the dissipative case when γ = 1.5, one has τ_geo ≲ t_obs after ∼400 s. Therefore, at late times we have δt ≲ t_obs. Thus, photospheric blending may potentially act as a filter, suppressing small timescale variability while leaving large timescale variability unaffected. We leave a detailed investigation for a future time.

The effect of photospheric blending is general to photospheric models. Specifically, it should be present during the prompt emission as well. However, it is likely negligible during this phase as the geometrical timescale is expected to be small for standard GRB parameters. It becomes significant in cases when Γ is small, for example when t ≫ t_b in the current framework, or in cases when the kinetic power is exceptionally high.

4.4. Assumptions in previous works

In this section, we mention some underlying assumptions in the works of Ronchini et al. (2021) and Tak et al. (2023), which may influence our results.

In Tak et al. (2023), the peak energy is found to decrease with time as E_p ∝ t⁻¹ during the decay phases of a large fraction of the GRB pulses studied. This decline is in agreement with the theoretical expectation from HLE. However, in their study, only data from the Fermi Gamma-ray Burst Monitor (GBM) are used. The low-energy threshold of GBM is 8 keV, with full effective area above ∼20 keV (Meegan et al. 2009). Therefore, E_p remains clearly visible in the GBM only a short time after the peak of the pulse. Thus, the performed analysis cannot probe very deep into the HLE-regime, with the peak energy in some of the studied pulses being tracked for ≲10 s. Since the behaviour E_p ∝ t⁻¹ is expected only once the line-of-sight contribution has faded, a robust conclusion about the decay rate of the peak energy may require additional observations in lower-energy bands. This may be possible for GRBs detected by the future space mission SVOM to be launched in 2024 (Wei et al. 2016). Its two gamma-ray instruments, ECLAIRs and GRM, offer a spectral coverage of the prompt emission from 4 keV to 5 MeV (Bernardini et al. 2017).

In Fig. 2, the evolution of the model spectra in the XRT band during the ESD is shown in comparison with data points obtained by Ronchini et al. (2021). However, such a comparison is not straightforward. The data points are generated by fitting a power-law function to the XRT spectral data during the ESD, and, thus, they include instrumental effects and background noise. The evolution of the photon index for the models on the other hand, is calculated using Eq. (10), following the definition of α_XRT in Ronchini et al. (2021). Equation (10) estimates the photon index by approximating the model spectra between 0.5 keV and 10 keV by a single power law (i.e. a straight line across the purple shaded region in Fig. 1). Such a prescription neglects the shape of the model spectrum within the XRT band. A fairer comparison against the data points would be achieved if one instead fitted a power-law function to mock data, where the mock data are generated by folding the model spectra through the XRT response matrix. It is plausible that the spectral index obtained using this method would differ from that obtained using Eq. (10), especially in regimes of low signal-to-noise.

Lastly, the photon index obtained when fitting the XRT data is highly sensitive on the modelling of the X-ray absorption below ≲1 keV. Absorption of X-rays occurs in the Milky Way, in the host galaxy of the GRB, and possibly also in the intergalactic medium (Wilms et al. 2000; Behar et al. 2011). X-ray spectra from distant sources are commonly fitted with an absorbed power law, where the absorption is modelled with a galactic plus an intrinsic hydrogen column density, N_H (Starling et al. 2013). The absorbed power-law model often gives good fits to the data. However, the value of the photon index obtained is sometimes highly degenerate with the best-fit value for the hydrogen column density (e.g. Valan et al. 2023). This point is indeed raised and discussed in Ronchini et al. (2021), who argue that their general results are robust against such a degeneracy. However, it highlights the importance of a correctly modelled absorption. Because of this, and the argumentation in the paragraph above, one should interpret the observational results presented in Fig. 2 with some caution.

5. Conclusion

In this paper we have investigated the end phase of the prompt emission in photospheric models of GRBs. Due to the short geometrical timescale expected in these models for typical GRB parameter values, the interpretation of the ESD as HLE is challenging. Instead, we interpret the ESD as an emission signature from the dying central engine.

We modelled the fading central engine by prescribing a power-law decay for the injected power, Ė_K, and the Lorentz factor, Γ, and constructed simple non-dissipative and dissipative frameworks to obtain the observed spectrum as a function of time (Fig. 1). In the dissipative case, we find that the photospheric emission from the dying central engine can mimic the spectral evolution predicted from HLE, if the kinetic power and the Lorentz factor decrease as t⁻³ and t^−7/8, respectively. If the Lorentz factor decreases more rapidly, as t^−1.5, the dissipative model can reproduce the α − F relation obtained by Ronchini et al. (2021), as shown in our Fig. 2. This requires the existence of a persistent, high-energy power law to account for the spectral shape at late times, indicating that dissipation should take place near the photosphere. In both cases, we find that if the kinetic power decreases more quickly than t⁻⁵, some fine-tuning of the other model parameters is necessary to be consistent with the observations. These results rely on a series of simplifying assumptions. A more detailed approach could include a time-varying efficiency and a physically motivated calculation of the spectral shape.

We used the deduced late-time behaviour of the central engine to gain insights into the progenitor systems (Sect. 4.2). We argue that the jet kinetic power decreasing as t⁻³ is in rough agreement with current state-of-the-art numerical simulations of compact binary mergers. If the comoving density in the jet remains constant, this in turn implies that Γ ∝ t^−1.5. For collapsars, the evolution Ė_K ∝ t⁻³ after the initial prompt phase requires a diffuse density profile near the stellar edge, where the mass density, ρ ∝ r^−a, decreases rapidly with a ∼ 6. Here we assumed free-fall accretion. The diffuse structure should extend from R ≲ 10¹⁰ cm to R ≳ 3 × 10¹⁰ cm to account for the observed duration of the ESD.

Lastly, we discussed photospheric blending: the fact that outer regions of the jet with high optical depths may obscure emission from inner, optically thin layers, leading to a blending of the radiation from the different regions. Although a generic feature in photospheric models of GRBs, this effect is likely negligible during the prompt phase for standard GRB parameter values, since its typical duration scales with the geometrical timescale of the dense regions (Eq. (13)). However, it should help suppress short timescale variability during the ESD in the current context, since the geometrical timescale can become significant at late times.

To conclude, spectral and temporal observations of the ESD can be reproduced by late-time photospheric emission from a dying central engine. To account for the observed light curve, the model predicts a decline in the kinetic power with a temporal index between ∼ − 3 and −5. A decline in the Lorentz factor with a temporal index ∼ − 7/8 reproduces a spectral behaviour similar to HLE, while a temporal index of ∼ − 1.5 reproduces the α − F relation. The behaviour should be quasi-universal to GRB central engines during their last stages, and the model predictions can be tested against long-lasting numerical simulations in the future.

Acknowledgments

We thank the anonymous referee for many good suggestions and the XRT team for providing the data and light curves. F.A. is supported by the Swedish Research Council (Vetenskapsrådet; 2022-00347). F.D. and R.M. acknowledge the Centre National d’Études Spatiales (CNES) for financial support in this research project.

References

All Tables

All Figures

Fig. 1. Snapshot spectra showing the evolution of the spectral luminosity during the ESD in the non-dissipative scenario (left) and the dissipative scenario (right). Dashed lines have γ = 7/8, while solid lines have γ = 1.5, as also indicated in the figure. The spectra are plotted at even intervals in time from the onset of the ESD at t/tb = 1 (dark colour) until t/tb = 10 (bright colour). Note that the two spectra for γ = 7/8 and γ = 1.5 overlap at t/tb = 1. The purple shading indicates the energy sensitivity of XRT. Parameter values used are given in Table 1.

In the text

	Fig. 2. Photon index as a function of Fmax/F along the ESD branch. The grey points are the best-fit values from time-resolved analyses of eight GRBs, as obtained in Ronchini et al. (2021). Parameter values used are given in Table 1.
In the text

	Fig. 3. Evolution of the Lorentz factor, peak energy, photospheric radius, and normalized XRT-flux during the ESD for the dissipative model. The XRT flux in the bottom-right panel is shown for the low-energy index, α, equal to −2/3 (dark colour), −1 (intermediate colour), and −1.5 (bright colour). Parameter values are given in Table 1.
In the text