GRB 260924A: Extreme Energetics, Bright GeV Emission, and Exceptional X-ray and Optical Luminosity

Introduction

Gamma-ray bursts (GRBs) are high-energy transients whose prompt emission exhibits diverse temporal structures, with rapid variability and durations ranging from milliseconds to minutes von Kienlin et al. 2020. Their emission can continue well beyond the brightest gamma-ray pulses, motivating observations across multiple energy bands and timescales Kumar & Zhang 2015.

Relativistic outflows provide a framework for interpreting GRB emission. The prompt emission may arise from internal dissipation within the outflow Rees & Meszaros 1994 and radiation released near its photosphere Mészáros & Rees 2000. As the ejecta interact with the surrounding medium, a forward shock propagates into it and a reverse shock can cross the ejecta Mészáros & Rees 1997; Kobayashi 2000. Electrons accelerated at these shocks may produce synchrotron emission Sari et al. 1998, while inverse-Compton scattering can contribute at higher photon energies Sari & Esin 2001. The relative contributions of these processes depend on the outflow properties and environment and can be constrained by the observed temporal and spectral evolution.

The combination of Fermi (Meegan et al. 2009; Atwood et al. 2009) and Swift Gehrels et al. 2004 provides a broad view of GRB emission from the prompt phase to the developing afterglow. The Fermi Gamma-ray Burst Monitor (GBM) and the Swift Burst Alert Telescope (BAT) detect and characterize the temporal and spectral structure of the prompt gamma-ray emission, while the Fermi Large Area Telescope (LAT) extends this view to GeV energies. LAT observations have revealed that this high-energy emission can show a delayed onset and persist well beyond the brightest lower-energy gamma-ray pulses Ajello et al. 2019. Following the trigger, the rapid response of Swift provides a complementary view at lower energies, with the X-Ray Telescope (XRT; Burrows et al. 2005) tracing the early X-ray evolution and the Ultraviolet/Optical Telescope (UVOT; Roming et al. [76]) identifying and monitoring the ultraviolet (UV) and optical counterpart. Ground-based optical and near-infrared (NIR) observations extend this multiwavelength coverage.

Multiwavelength observations provide complementary constraints on the radiation mechanisms and emitting regions respon- [superscript star] E-mail: [email protected] sible for the temporal and spectral structure of GRBs (Fraija 2015; Fraija et al. 2016; Ajello et al. 2019; Fraija et al. 2019b; Becerra et al. 2023; Dainotti et al. 2024). The combined afterglow evolution across different bands probes the properties of the relativistic ejecta and the density structure of the circumburst medium (Fraija et al. 2017). Departures from a simple decay can also reveal the need to consider additional emission components, continued energy supply (e.g., see Laskar et al. 2015; Zhao et al. 2021; Becerra et al. 2019; Pereyra et al. 2022; Fraija et al. 2021), or evolving shock microphysics (e.g., see Rossi & Rees 2003; Lemoine 2013; Fraija et al. 2024b). Broad spectral and temporal coverage therefore allows afterglow models to be tested against several observational constraints within a common physical framework.

GRB 260924A was detected by both Fermi (Fermi GBM Team 2026) and Swift (D’Elia et al. 2026). Its GBM and BAT light curves show a bright, structured initial episode followed by a second prominent gamma-ray pulse beginning near 334 s after the GBM trigger (Gupta et al. 2026a). Between these episodes, BAT detects weaker emission, while Fermi/LAT reveals high-energy activity extending to tens of GeV (Holzmann Airasca et al. 2026). XRT observed the decline following the initial activity and the rise associated with the second pulse, although its early coverage ended before that episode was complete. Swift/UVOT and ground-based observations provided optical and NIR follow-up, complementing the XRT coverage and enabling multiwavelength monitoring over the following days. Optical spectroscopy revealed both the Ly$\alpha$ break and the Lyman limit, yielding a redshift of $z = 3.51 \pm0.10$ (Castro-Tirado et al. 2026). Using $z = 3.51$ and the time-integrated GBM spectrum, Gupta et al. (2026b) estimated an isotropic-equivalent gamma-ray energy of $E_{\gamma,\mathrm{iso}} \simeq9.74 \times10^{54}$ erg in the 1–10,000 keV rest-frame band, placing GRB 260924A among the most energetic GRBs detected and at an energy scale comparable to that of GRB 221009A (see e.g. Burns et al. 2023; O’Connor et al. 2023; Lesage et al. 2023). The combination of a large prompt energy release, extended GeV activity, and extensive afterglow coverage makes GRB 260924A a useful case for examining the relationship between the prompt emission and the subsequent broadband evolution.

In this paper, we present a multiwavelength analysis of GRB 260924A, combining Fermi/GBM and LAT observations with Swift/BAT, XRT, and UVOT data and ground-based optical and NIR follow-up. We characterize the temporal and spectral evolution of the two bright gamma-ray episodes, the intervening activity, and the subsequent afterglow. We use the light curves, spectral measurements, and LAT photon energies to assess possible radiation mechanisms and emitting regions. We model the afterglow with forward-shock synchrotron and synchrotron self-Compton (SSC) emission, an early reverse-shock SSC contribution, and energy injection with evolving shock microphysics. We also place the event in the context of GRB comparison samples using rest-frame X-ray luminosities, optical light curves transformed to a common redshift, and observed LAT photon fluxes.

The paper is organized as follows. Section 2 describes the observations, data analysis, and common temporal reference. § 3 presents the temporal evolution across the observed energy bands. § 4 presents the BAT, GBM, and LAT spectral analyses and selected LAT high-energy events. § 5 presents our external-shock model of the multiwavelength afterglow and the inferred physical parameters. § 6 discusses the physical implications of these results, and § 7 summarizes our conclusions.

Observations and Data Analysis

This section describes the observations and analysis procedures used to construct the multiwavelength data set of GRB 260924A, including the processing of the Fermi and Swift data and the compilation of ground-based optical and NIR photometry. Unless otherwise stated, times are given in the observer frame relative to the Fermi/GBM trigger, adopted as $T_0$.

Fermi/GBM

Light curves

GRB 260924A was detected by the Fermi/GBM instrument at 2026 September 24 01:13:23 UTC ($T_0$) (Fermi GBM Team 2026). We constructed the GBM count-rate light curves directly from time-tagged event (TTE) data in 2-s bins relative to $T_0$ (Figure 1). We selected PHA channels whose energy boundaries lay entirely within the nominal bands of 8–50 and 50–300 keV for NaI and 300–5000 keV for BGO. Count rates were corrected for detector deadtime. For each detector and band, a linear background was fitted over [−125, −35] and [250, 320] s relative to $T_0$ and subtracted. The resulting rates from NaI detectors n3, n4, and n5 were summed within each band, while BGO b0 was displayed separately. Uncertainties include counting statistics and the covariance of the fitted background parameters.

Although automated methods such as Bayesian blocks can be used to define time intervals, we selected the boundaries by visual inspection of the GBM light curves to delimit the main emission episodes and facilitate comparison with BAT and LAT. The adopted boundaries are marked in Figure 1.

Spectral analysis

For the first bright episode, we analysed the spectrum integrated over −6.1–164.9 s and the three subintervals −6.1–50, 50–115, and 115–164.9 s relative to $T_0$. We used NaI detectors n3, n4, and n5 and BGO detector b0, with nominal energy ranges of 8–900 keV and 0.3–40 MeV, respectively.

For each detector, we constructed an effective response by averaging the time-dependent RSP2 response over the corresponding interval, weighting by the total observed counts.

The background uncertainties, including the covariance of the fitted coefficients, were propagated to the integrated background counts. Further details of the response construction, background estimation, and exposure calculation are given in Appendix A.

Spectral fitting used a Poisson–Gaussian profile likelihood, combining Poisson-distributed observed counts with Gaussian constraints on the estimated background counts. We report the fit statistic as PG = 2S, where S is the negative log-likelihood up to data-dependent constants. We modelled the continuum using either the Band function (Band et al. 1993), which smoothly joins a low-energy power law with an exponential cutoff to a high-energy power law, or the Comptonized model, defined here as a power law with an exponential cutoff. We also tested whether adding a blackbody (BB) improved the Band or Comptonized fits (Ryde 2004; Ryde & Pe’er 2009). The BB prescription and parameter bounds are given in Appendix A5.

To examine the earliest spectral evolution, we divided the −6.1 to 50 s interval at 12.29 s, selecting this boundary from the 0.512-s GBM light curve before fitting the spectra. The second bright episode was analysed over 332.8–390 s, with 334–390 s used as an alternative extraction interval. Spectral fits were performed in detector count space; unfolded $\nu F_{\nu}$ spectra are shown for visualization. Tests of sensitivity to the response treatment, background model, and channel selection are described in Appendix A.

The spectra displayed in Figure 2 and Figure 3 and the parameters listed in Table 1 and Table 2 use linear temporal backgrounds; quadratic backgrounds were examined as a sensitivity test.

Interval (s)	Model	$\alpha$	$\beta$	$E_{\mathrm{peak}}$ (keV)	$k_{\mathrm{B}}T$ (keV)	PG	$N_{\mathrm{ch}}$	$N_{\mathrm{src}}$
−6.1–164.9	Band	$-0.85^{+0.01}_{-0.01}$	$-3.52^{+0.35}_{-0.64}$	$596^{+12}_{-12}$	–	687.25	485	4
	Band+BB	$-0.91^{+0.02}_{-0.02}$	$-3.79^{+0.46}_{-1.59}$	$681^{+27}_{-26}$	$42.7^{+3.5}_{-3.4}$	666.31	485	6
−6.1–50$^{a}$	Band	$-0.65$	$-12^{b}$	350.7	–	524.22	485	4
	Band+BB	$-0.77$	$-12^{b}$	453.9	33.10	507.30	485	6
50–115	Band	$-0.87^{+0.01}_{-0.01}$	$-3.16^{+0.19}_{-0.27}$	$764^{+18}_{-17}$	–	715.46	485	4
	Band+BB	$-0.91^{+0.02}_{-0.02}$	$-3.27^{+0.23}_{-0.35}$	$836^{+37}_{-34}$	$60.2^{+9.9}_{-8.7}$	708.59	485	6
115–164.9	Band	$-0.89^{+0.02}_{-0.01}$	$-4.15^{c}$	$504^{+16}_{-17}$	–	501.52	485	4
332.8–390	Comptonized	$-1.05^{+0.03}_{-0.03}$	–	$714^{+64}_{-55}$	–	355.92	391	3
	Band	$-1.03^{+0.03}_{-0.03}$	$-2.15^{d}$	$617^{+73}_{-72}$	–	349.15	391	4

Interval (s)	Model	$A/(10^{-2}\ \mathrm{ph\ cm^{-2}\ s^{-1}\ keV^{-1}})$	$A_6$	$\Delta PG$
−6.1–164.9	Band	$1.810^{+0.017}_{-0.016}$	–	–
	Band+BB	$1.624^{+0.043}_{-0.042}$	$1.33^{+0.39}_{-0.34}$	20.94
−6.1–50$^{a}$	Band	1.19	–	–
	Band+BB	0.85	3.8	16.93
50–115	Band	$2.42^{+0.02}_{-0.02}$	–	–
	Band+BB	$2.26^{+0.06}_{-0.06}$	$0.47^{+0.26}_{-0.21}$	6.88
115–164.9	Band	$1.88^{+0.03}_{-0.03}$	–	–
332.8–390	Comptonized	$0.82^{+0.02}_{-0.02}$	–	–
	Band	$0.85^{+0.03}_{-0.03}$	–	6.77

Fermi/LAT

We analysed Fermi/LAT events within a 15° region of interest (ROI) centred on the Swift/XRT position, (RA, Dec) = (133.76°, 12.04°), over the 0.1–1000 GeV energy range. We selected event class 8 and event type 3, applied a maximum zenith angle of 100°, and used the P8R3_TRANSIENT020E_V3 instrument response functions. We used an unbinned maximum-likelihood analysis with Fermitools version 2.5.3 and fermtools-data version 0.18.

We analysed five non-overlapping intervals spanning $T_0 + 70.04$ to $T_0 + 1100$ s (Table 3), covering the initial episode, inter-episode activity, the second pulse, and the subsequent high-energy tail. Alternative overlapping intervals were used to test sensitivity to the temporal selection and were not included as independent light-curve measurements (§ 3.2).

Interval (s)	$\Gamma_{\mathrm{LAT}}$	$F_{\mathrm{ph}}$ ($10^{-5}\ \mathrm{ph\ cm^{-2}\ s^{-1}}$)	$F_E$ ($10^{-8}\ \mathrm{erg\ cm^{-2}\ s^{-1}}$)	TS
70.04–164.9	$2.21 \pm0.08$	$60.25 \pm4.38$	$47.43 \pm8.95$	1058.36
164.9–332.8	$2.24 \pm0.14$	$10.01 \pm1.23$	$7.40 \pm2.35$	282.41
332.8–390	$2.23 \pm0.34$	$5.06 \pm1.48$	$3.78 \pm2.96$	63.55
390–800.53	$2.28 \pm0.15$	$3.20 \pm0.43$	$2.17 \pm0.67$	309.28
800.53–1100	$2.52 \pm0.45$	$2.00 \pm0.72$	$0.93 \pm0.52$	26.95

The background model included the Galactic diffuse template gll_iem_v07.fits, the isotropic component iso_P8R3_TRANSIENT020E_V3_v1.txt, and 243 sources from the 4FGL-DR4 catalogue within 25° of the Swift/XRT position. The same catalogue sources were included in all five principal intervals and were modelled as point sources with their positions and spectral parameters fixed. The GRB spectrum was described by a power law, $\mathrm{d}N/\mathrm{d}E \propto E^{-\Gamma_{\mathrm{LAT}}}$, where $\Gamma_{\mathrm{LAT}}$ is the photon index. Under the source hypothesis, the free parameters were the GRB normalization and photon index and the normalizations of the Galactic and isotropic diffuse components. The GRB position and diffuse spectral shapes were held fixed. All events satisfying the selection criteria were included in the likelihood analysis without an additional cut on their fitted source-association probabilities.

For each interval, we evaluated the likelihood-ratio test statistic (TS) by comparing models with and without the GRB, as described in Appendix B. Both diffuse normalizations were reoptimized under each hypothesis within the bounds [0, 30]. The quoted photon-flux and photon-index errors are nominal $1\sigma$ statistical uncertainties derived from the local covariance matrix.

Energy fluxes were obtained by integrating the energy-weighted fitted spectra over 0.1–1000 GeV. Their uncertainties were propagated from the fitted covariance matrix, including the covariance between normalization and photon index. For the light-curve display in Figure 4, we integrated the same fitted models over 0.1–300 GeV; the likelihood fits, Table 3, and temporal analysis use 0.1–1000 GeV.

We also examined the energies and arrival times of individual LAT events. Their model-dependent source-association probabilities, $P_{\mathrm{GRB}}$, were obtained from a fit over 0–1100 s after the GBM trigger (hereafter the full-window fit). The bottom-left panel of Figure 4 displays all 404 selected events, with greyscale indicating these probabilities. Table 4 lists the ten highest-energy events with $P_{\mathrm{GRB}} \ge0.9$. Their properties are discussed in § 4.2.

$E_{\mathrm{obs}}$ (GeV)	$t-t_{\mathrm{GBM},0}$ (s)	Separation (deg)	$P_{\mathrm{GRB}}$
43.31	82.24	0.070	0.999992
24.95	122.92	0.41	0.998767
14.24	304.15	0.08	0.999993
11.59	141.64	0.18	0.999957
5.30	665.55	0.05	0.999996
3.95	98.99	0.28	0.999857
3.52	88.12	0.61	0.999316
3.11	587.17	0.11	0.999985
2.63	136.79	0.24	0.999920
2.36	147.86	0.23	0.999937

Further implementation details and sensitivity tests are provided in Appendix B.

Swift/BAT

Swift/BAT triggered on and located GRB 260924A ((D’Elia et al. 2026)). Throughout this work, Swift times were referred to the Fermi/GBM reference time using $t_{\mathrm{GBM}} = t_{\mathrm{BAT}} + 4.71$ s.

The BAT count-rate light curve was constructed from mask-weighted rates in the 15–350 keV band using 1-s bins and is shown in the upper panel of Figure 5. It was used to compare the temporal structure with that observed by GBM and XRT (§ 3.3). Spectral fits were performed over 15–150 keV.

We extracted BAT spectra over $T_0 + 70$ to $T_0 + 165$ s, $T_0 + 250$ to $T_0 + 320$ s, and $T_0 + 334$ to $T_0 + 390$ s. These intervals sample the initial bright episode, the weaker inter-episode activity, and the second bright pulse, respectively. We fitted each spectrum with a power-law model and examined whether adding an exponential cutoff improved the fit. The spectral responses accounted for changes in spacecraft pointing within each interval. Table 5 summarizes the power-law fits.

Interval (s)	$\Gamma_{\mathrm{BAT}}$	$F_{15-150}$ ($10^{-7}\,\mathrm{erg\,cm^{-2}\,s^{-1}}$)	$\chi^2/\mathrm{dof}$
70–165	$1.08^{+0.02}_{-0.02}$	$3.14 \pm0.04$	36.30/56
250–320	$1.55^{+0.28}_{-0.28}$	$0.08^{+0.01}_{-0.01}$	46.40/56
334–390	$1.16^{+0.03}_{-0.04}$	$1.48 \pm0.03$	53.49/56

To examine spectral evolution within the two bright episodes, we also extracted spectra in shorter subintervals and fitted them with power-law models. The resulting photon indices, $\Gamma_{\mathrm{BAT}}$, are shown in the lower panels of Figure 5 and discussed in § 4.1.

The spectra were fitted using the $\chi^2$ statistic, including channel-dependent systematic uncertainties. We report 90 per cent confidence intervals for one parameter of interest, corresponding to $\Delta\chi^2 = 2.706$. Energy-flux confidence intervals were obtained by profiling over the photon index, accounting for the normalization–index covariance. Details of the treatment and verification of the BAT systematic uncertainties are provided in Appendix § 12.

Swift/XRT

Swift/XRT began its pointed observations of GRB 260924A approximately 110.6 s after the BAT trigger, detecting a bright, previously uncatalogued X-ray source within the BAT localization ((D’Elia et al. 2026)). Subsequent astrometric refinement yielded a position of RA = 133.76° and Dec = +12.04° (J2000), with a 90 per cent confidence radius of 2.1″ ((Evans et al. 2026)). The refined XRT analysis included 6.4 ks of observations spanning approximately 99 s to 43.5 ks after the BAT trigger, comprising 365 s in Windowed Timing (WT) mode, including the initial slew data, followed by observations in Photon Counting (PC) mode ((Dicharia et al. 2026)).

We obtained the Swift/XRT 0.3–10 keV count-rate and flux light curves, together with the hardness-ratio products, from the UK Swift Science Data Centre (UKSSDC) GRB repository ((Evans et al. 2009, 2010))¹. Figure 6 shows the count-rate evolution and the photon index $\Gamma_{\mathrm{XRT}}$ inferred from the hardness ratio. Table 6 reproduces the automated temporal fit reported by the catalogue, retaining its original BAT time reference. Times plotted in Figure 6 were converted to the GBM reference using the offset specified in § 2.3.

Segment	$\alpha_i$	Break after segment, $t_{b,i}$ (s; BAT $T_0$)
1	$0.3^{+0.5}_{-0.7}$	$130^{+8}_{-6}$
2	$2.85^{+0.12}_{-0.12}$	$256.4^{+15.7}_{-11.1}$
3	$1.316^{+0.025}_{-0.029}$	$4.67^{+0.22}_{-0.44}\times10^3$
4	$-1.50^{a}$	$6.30^{+4.87}_{-0.26}\times10^3$
5	$1.61^{+0.05}_{-0.05}$	—

The X-ray temporal analysis and comparison with the optical/NIR afterglow are presented in § 3.4, with the late-time break test described in Appendix § 13.1.

To characterize the early X-ray spectral evolution outside the flare interval, we analysed three spectra accumulated over 263.5–314.0 s (WT), 3652.5–3757.0 s (WT), and 3763.4–4656.0 s (PC) after the GBM trigger. We fitted these spectra jointly over 0.3–10 keV with an absorbed power-law model. The Galactic absorption column was fixed, and the intrinsic column at $z = 3.51$ was tied across the three spectra, while their photon indices and normalizations varied independently.

The fits used the W statistic, with spectra grouped to at least one source count per bin. Parameter uncertainties correspond to 90 per cent confidence for one parameter of interest ($\Delta W = 2.71$), with the remaining parameters reoptimized. Table 7 summarizes the fitted parameters and absorption-corrected flux densities. Figure 7 shows the spectra, the models folded through their instrumental responses, and the corresponding residuals. Additional regrouping was applied solely for visualization and did not affect the spectral fits.

Mode	Time since GBM trigger (s)	Exposure (s)	$\Gamma_{\mathrm{XRT}}$	$F_{\nu,\mathrm{unabs}}$ (1 keV) ($\mu$Jy)
WT	263.5–314.0	50.45	$1.78^{+0.07}_{-0.07}$	$631^{+38}_{-36}$
WT	3652.5–3757.0	104.36	$2.25^{+0.15}_{-0.14}$	$26.0^{+2.4}_{-2.2}$
PC	3763.4–4656.0	889.03	$2.01^{+0.10}_{-0.10}$	$18.4^{+1.4}_{-1.3}$

UV, optical, and near-infrared observations

We compiled the available UV, optical, and NIR photometry of GRB 260924A reported by multiple follow-up facilities. The data set includes observations obtained with the Ultraviolet/Optical Telescope (UVOT) on board Swift (Klingler & D’Elia 2026); the 1.5-m telescope at Sierra Nevada Observatory (OSN) (Hu et al. 2026); the Nordic Optical Telescope (NOT) with ALFOSC (Izzo et al. 2026); the Lulin One-meter Telescope (LOT) and Super Light Telescope (SLT) (Aryan et al. 2026); the GROWTH-India Telescope (Vijaykumar et al. 2026); the SYSU 80-cm telescope ([61, 58]); the Ondřejov D50 telescope (Novotny et al. 2026); the Calar Alto 2.2-m telescope with CAFOS (Castro-Tirado et al. 2026); the 1.3-m COLIBRÍ telescope with the two-channel DDRAGO imager (Fortin et al. 2026); the MDM 1.3-m telescope with MIRAGE (Mo et al. 2026); and AbAO (Pankov et al. 2026).

We placed the photometry on the AB magnitude system and corrected it for Galactic foreground extinction adopting $E(B - V) = 0.041$ mag (Schlafly & Finkbeiner 2011) and $R_V = 3.1$ (Fitzpatrick 1999). For visualization, measurements in related passbands were displayed together as $g/g^\prime$ and $r/r^\prime/R$. The temporal fit reported in § 3.4 used only the six $r/r^\prime$ detections, while the original filter information was retained for the spectral energy distribution analysis. Upper limits were retained at the confidence level reported by the corresponding observations.

The upper panel of Figure 8 shows the compiled UV, optical, and NIR light curves together with the XRT flux-density light curve at 1 keV. The temporal evolution and optical-to-X-ray spectral energy distributions are discussed in § 3.4.

Multiwavelength temporal evolution

We characterize the multiwavelength evolution of GRB 260924A across the two bright gamma-ray episodes, the intervening activity, and the subsequent afterglow. The temporal comparisons and optical-to-X-ray spectral energy distributions (SEDs) provide the observational basis for the modelling presented in § 5.

Fermi/GBM and Swift/BAT light curves

Figure 1 and Figure 5 show the GBM and BAT count-rate light curves. The initial bright, structured GBM episode extends to approximately $T_0 + 160$–165 s, followed by a distinct second pulse beginning near $T_0 + 334$ s that is also evident in BAT. Between these episodes, a candidate GBM excess at approximately $T_0 + 195$–235 s is most apparent in the NaI detectors, while BAT also traces weaker emission. Count-rate amplitudes are not directly comparable between instruments because their energy ranges and responses differ.

Because $T_0 + 250$–320 s was used to fit the GBM background, the background-subtracted light curve does not provide an independent test for source emission during that interval.

Fermi/LAT light curve and inter-episode emission

Figure 4 shows the LAT light curve for the five intervals analysed in § 2.2. The corresponding photon and energy fluxes are listed in Table 3.

LAT emission was detected between the two bright gamma-ray episodes. During $T_0 + 250$ to $T_0 + 320$ s, the fitted source at the XRT position has TS = 87.16, a 0.1–1000 GeV photon flux of $(5.78 \pm1.41) \times10^{-5}$ ph cm$^{-2}$ s$^{-1}$, and $F_{\mathrm{LAT}} = 2.06 \pm0.24$.

Subdividing the first principal interval into $T_0 + 104$ s reveals a decrease in the average 0.1–1000 GeV photon flux from $(11.45 \pm1.01) \times10^{-4}$ over 70.04–104 s to $(3.25 \pm0.39) \times10^{-4}$ ph cm$^{-2}$ s$^{-1}$ over 104–164.9 s, relative to $T_0$. The lower average flux in the later subinterval indicates an early decline, although these two measurements do not resolve the peak time. Neither these subdivisions nor the overlapping 250–320 s interval were included as additional independent points in the five-bin temporal fit.

The four principal intervals after $T_0 + 164.9$ s are described by $F_{\mathrm{ph}}(t) = A(t/164.9\ \mathrm{s})^{-\alpha_{\mathrm{LAT}}}$, with $\alpha_{\mathrm{LAT}} = 1.30^{+0.20}_{-0.19}$, $A = (1.60^{+0.50}_{-0.27}) \times10^{-4}$ ph cm$^{-2}$ s$^{-1}$, and $\chi^2/\mathrm{dof} = 0.31/2$. The model was averaged over each annual interval using normalized exposure weights. Parameter uncertainties correspond to $\Delta\chi^2 = 1$, with the other parameter reoptimized.

The first interval, excluded from this fit, has a measured photon flux of $(6.03 \pm0.44) \times10^{-4}$ ph cm$^{-2}$ s$^{-1}$, compared with an extrapolated interval-averaged prediction of $(2.68 \pm0.67) \times10^{-4}$ ph cm$^{-2}$ s$^{-1}$. Their ratio, $2.25 \pm0.59$, quantifies the early excess relative to the subsequent decay. The fit includes the interval containing the second bright episode, and its uncertainties are conditional on the adopted spectral models and exposure weights. Further details and sensitivity tests are given in Appendix B5.

Temporal structure and hardness evolution of the second bright episode

A second prominent gamma-ray pulse begins near $T_0 + 334$ s in both the GBM and BAT light curves. XRT records a contemporaneous flare, with an observed maximum near $T_0 + 344.5$ s followed by a partial decline. Its early coverage ends near $T_0 + 359.4$ s, preventing us from following the complete X-ray decay.

Figure 9 compares the pulse profiles and hardness evolution across GBM, BAT, and XRT. The GBM pulse is broader in 8–50 keV than in 50–300 keV. The hardness measurements are listed in Table 8.

Table 8. Count-rate hardness ratios of the second bright episode shown in Figure 9. Each instrument retains its own time bins; entries sharing a row are not necessarily simultaneous measurements.

Fermi/GBM Interval (s)	$H_{\mathrm{GBM}}$	Swift/BAT Interval (s)	$H_{\mathrm{BAT}}$	Swift/XRT Interval (s)	$H_{\mathrm{XRT}}$
334–338	$1.32^{+0.30}_{-0.22}$	333.93–337.93	$0.96^{+0.16}_{-0.13}$	321.11–324.52	$0.63 \pm0.12$
338–342	$1.85^{+0.12}_{-0.11}$	337.93–341.93	$1.06^{+0.07}_{-0.06}$	324.52–328.53	$0.96 \pm0.18$
342–346	$1.75^{+0.04}_{-0.04}$	341.93–345.93	$1.05^{+0.03}_{-0.03}$	328.53–331.67	$0.99 \pm0.18$
346–350	$1.35^{+0.05}_{-0.05}$	345.93–349.93	$0.95^{+0.03}_{-0.03}$	331.67–334.39	$1.03 \pm0.18$
350–354	$1.28^{+0.07}_{-0.07}$	349.93–353.93	$0.93^{+0.04}_{-0.04}$	334.38–337.04	$0.94 \pm0.17$
354–358	$1.06^{+0.09}_{-0.08}$	353.93–357.93	$0.86^{+0.05}_{-0.05}$	337.04–338.76	$1.37 \pm0.26$
358–362	$0.77^{+0.09}_{-0.08}$	357.93–361.93	$0.81^{+0.08}_{-0.07}$	338.76–340.46	$1.69 \pm0.26$
362–366	$0.75^{+0.10}_{-0.09}$	361.93–365.93	$0.82^{+0.10}_{-0.09}$	340.46–342.54	$1.63 \pm0.27$
366–370	$0.93^{+0.22}_{-0.17}$	365.93–369.93	$0.96^{+0.16}_{-0.13}$	342.54–344.35	$1.90 \pm0.29$
370–374	$1.14^{+0.43}_{-0.28}$	369.93–373.93	$0.57^{+0.12}_{-0.11}$	344.35–345.57	$1.77 \pm0.29$
374–378	$0.65^{+0.43}_{-0.27}$	373.93–377.93	$1.22^{+0.46}_{-0.31}$	345.57–347.13	$1.77 \pm0.28$
378–382	$0.38^{+0.35}_{-0.25}$	377.93–381.93	$0.54^{+0.16}_{-0.14}$	347.13–349.11	$1.98 \pm0.30$
382–386	$0.84^{+0.69}_{-0.35}$	381.93–385.93	$1.17^{+0.41}_{-0.28}$	349.11–351.22	$1.66 \pm0.27$
386–390	$0.62^{+1.23}_{-0.45}$	385.93–389.93	$1.20^{+1.11}_{-0.50}$	351.22–354.01	$2.28 \pm0.35$
–	–	–	–	354.01–356.61	$1.78 \pm0.29$
–	–	–	–	356.61–360.10	$1.62 \pm0.22$

GBM and BAT show decreasing hardness during the well-sampled early decay following the pulse maximum. XRT shows hardening during the rise, but its limited coverage prevents characterization of the spectral evolution throughout the full decay.

We fitted the 19 original WT bins covering 344.77–359.37 s, excluding the peak bin, with a power-law flare component proportional to $(t - 332.8\ \mathrm{s})^{-\alpha_{\mathrm{XRT}}}$. The underlying emission was represented by $B_{300}(t/300\ \mathrm{s})^{-\beta}$, with its normalization and slope jointly constrained by 30 pre-flare bins spanning 263.50–314.04 s. The fit yields $\alpha_{\mathrm{XRT}} = 1.11 \pm0.17$, with uncertainties quoted at the $1\sigma$ level ($\Delta\chi^2 = 1$ for one parameter of interest). The power-law model gives $\chi^2/\mathrm{dof} = 28.73/45$, whereas a constant excess above the underlying emission gives $\chi^2/\mathrm{dof} = 70.57/46$. A reconstruction using the time-dependent spectral conversion gives a compatible flux decay. The temporal index depends strongly on the assumed time origin, and the XRT hardness is consistent with remaining constant during the observed decline.

UV, optical, NIR, and X-ray evolution

The XRT light curve shows an initial decline interrupted by the flare discussed in § 3.3. After the flare, the observations resume with a declining light curve, followed by a rebrightening at several kiloseconds (Figure 6). The automated catalogue fit describes this feature with a rising fourth segment, $\alpha_4 = -1.50$, between breaks near 4.7 and 6.3 ks (Table 6). Hereafter, we refer to the afterglow phases before and after this rebrightening as the early and late afterglow, respectively.

The upper panel of Figure 8 compares the UV, optical, and NIR light curves compiled in § 2.5 with the XRT evolution at 1 keV. The observations span from a few minutes to several days after the burst, allowing us to characterize the late-time decay and reconstruct the optical-to-X-ray spectral energy distribution at selected epochs.

A single-power-law fit to six $r/r'$ detections spanning 0.19–2.44 d gives $\alpha_{\rm opt} = 1.46 \pm0.01$ for $F_\nu\propto t^{-\alpha_{\rm opt}}$, with $\chi^2/\mathrm{dof} = 25.18/4$. The excess scatter indicates that this slope characterizes the overall decline but does not fully describe the optical evolution.

For comparison, a power-law fit to 69 late-time XRT PC-mode measurements with bin centres spanning 0.17–3.40 d yields $\alpha_X = 1.57 \pm0.04$ for $F_\nu\propto t^{-\alpha_X}$, with $\chi^2 = 61.83$ for 67 degrees of freedom. The XRT fluxes were converted to monochromatic flux densities at 1 keV. Details of the XRT flux-density conversion, temporal binning, and broadband SED reconstruction are provided in Appendix § 13.

To test for a late-time steepening, we used the native, absorption-corrected 0.3–10 keV XRT fluxes rather than the reconstructed 1-keV light curve. Their count-to-flux conversion tracks the observed hardness ratio. Excluding exposures flagged by the UKSSDC for unreliable centroiding, a single power law gives $\alpha_{X,\mathrm{band}} = 1.55^{+0.04}_{-0.04}$ and $\chi^2/\mathrm{dof} = 66.32/67$. A continuous broken power law constrained to steepen improves the statistic by only $\Delta\chi^2 = 0.17$ for two additional parameters. Including flagged exposures increases the improvement to 5.41, but the break time remains poorly constrained. These products therefore do not establish a robust late-time X-ray break. The selection, fitting procedure, and sensitivity checks are described in Appendix § 13.1.

The lower panel of Figure 8 shows the optical-to-X-ray SEDs reconstructed for $0.18 < t < 0.20$, $0.85 < t < 1.00$, and $2.30 < t < 2.70$ d after the GBM trigger. These intervals are marked by shaded regions in the upper panel. The optical/NIR measurements were shifted to a common representative epoch within each interval, at which the XRT flux density was evaluated from the late-time temporal fit (Appendix § 13). Photometric upper limits are shown for completeness but were not included in the spectral fits.

Table 9 summarizes the reconstructed SEDs. The optical/NIR-only spectral indices are flatter than those obtained when XRT is included. The joint single-power-law fits yield $\chi^2/\mathrm{dof} = 17.41/1$, 8.80/1, and 4.70/2 for the three epochs, respectively, providing a poor description at the first two epochs under the adopted corrections and uncertainties. The XRT flux densities are predictions from the same late-time temporal fit and are therefore correlated; the latter two epochs lie within XRT coverage gaps.

Table 9. Optical/NIR and optical-to-X-ray spectral indices for the three SED epochs.

Interval (d)	$\beta_{\mathrm{opt/NIR}}$	$\beta_{\mathrm{opt/NIR/X}}$	$\chi^2/\mathrm{dof}$
0.18–0.20	$0.57 \pm0.07$	$0.85 \pm0.01$	17.41/1
0.85–1.00	$0.71 \pm0.06$	$0.89 \pm0.01$	8.80/1
2.30–2.70	$0.68 \pm0.10$	$0.90 \pm0.01$	4.70/2

Time-resolved gamma-ray spectra and LAT high-energy events

We examine the gamma-ray spectra across the main emission episodes and the intervening activity, together with the energies and arrival times of LAT photons. These results complement the temporal analysis in § 3 and provide spectral constraints for the modelling in § 5.

BAT spectra across the early episodes

Table 5 summarizes the power-law fits to the 15–150 keV BAT spectra over 70–165, 250–320, and 334–390 s after the GBM trigger. The inter-episode emission has a lower energy flux and a softer best-fitting spectrum than either bright episode, although its photon index is less precisely constrained.

For all three intervals, the cutoff-power-law fits converge to the power-law limit. The BAT spectra therefore do not constrain a finite cutoff energy or provide a BAT-only measurement of $E_{\rm peak}$.

We also tested for spectral evolution within the two bright BAT intervals using shorter spectra and their corresponding responses (Figure 5). In the first episode, the fitted power-law index increased from $1.07^{+0.03}_{-0.03}$ at 70–90 s to $1.21^{+0.04}_{-0.04}$ at 140–165 s. During the second pulse, it increased from $1.00^{+0.05}_{-0.04}$ at 334–344 s to $1.50^{+0.30}_{-0.29}$ at 379–390 s. Uncertainties are quoted at 90 per cent confidence. Compared with a shared photon index within each episode, allowing independent indices in its subintervals improves the fit by $\Delta\chi^2 = 31.93$ for three additional parameters in the first episode and $\Delta\chi^2 = 47.61$ for four in the second. Each subinterval retains an independent normalization in both comparisons. These results are conditional on the adopted diagonal treatment of systematic uncertainties.

These results support spectral softening within both bright episodes, whereas Table 5 reports their time-averaged spectral properties.

LAT high-energy events and spectral evolution

The bottom-left panel of Figure 4 shows the energies and arrival times of the selected LAT events. Among those with $P_{\rm GRB} \geq0.9$, 28 exceed 1 GeV and four exceed 10 GeV. The earliest event satisfying $P_{\rm GRB} \geq0.9$ arrived at 70.39 s with an energy of 0.306 GeV. The highest-energy event, 43.31 GeV, arrived at 82.24 s. Other prominent events include 24.95 GeV at 122.9 s and 14.24 GeV at 304.15 s, the latter between the two bright gamma-ray episodes. Table 4 lists the ten highest-energy events and their source-association probabilities. Their physical interpretation is discussed in § 6.1.

The middle-left panel of Figure 4 shows the fitted power-law spectra for the five principal intervals, expressed as $\nu F_\nu\propto E^{2-\Gamma_{\rm LAT}}$ over 0.1–1000 GeV. The corresponding photon indices are shown in the middle-right panel and listed in Table 3. The model envelopes reflect the local normalization–index covariance.

A joint likelihood fit with a shared photon index across the five intervals gives $\Gamma_{\rm LAT} = 2.24 \pm0.06$, with nominal one-parameter profile uncertainties defined by $2\Delta\ln\mathcal{L} = 1$. Allowing independent indices improves the fit by only $2\Delta\ln\mathcal{L} = 0.63$ for four additional parameters. The spectra are therefore compatible with a shared photon index at the adopted temporal resolution. Sensitivity tests are described in Appendix B.

Over the same intervals, the 0.1–1000 GeV energy flux decreases from $(47.43 \pm8.95) \times10^{-8}$ to $(0.93 \pm0.52) \times10^{-8}\ \mathrm{erg\ cm^{-2}\ s^{-1}}$ (Table 3).

GBM spectra of the bright episodes

The left panel of Figure 2 shows the GBM spectrum integrated over the first bright episode, from $-6.1$ to $164.9$ s relative to the GBM trigger. With the adopted linear background, adding a blackbody to the Band model improves the fit by $\Delta\mathrm{PG}=20.94$ for two additional parameters, with $k_{\mathrm{B}}T=42.7^{+3.5}_{-3.4}$ keV (Table 1). The sensitivity to the background prescription is discussed in Appendix A.

For the early $-6.1$–$50$ s interval, adding a BB improves the fit by $\Delta\mathrm{PG}=16.93$ relative to either Band or Comptonized with the adopted linear background, with $k_{\mathrm{B}}T\simeq33$ keV for Band+BB. Here, $\Delta\mathrm{PG}$ denotes the decrease in PG after adding the BB component. The Band+BB high-energy index reaches its lower bound, effectively approaching a Comptonized+BB spectrum over the observed energy range. Table 1 and Table 2 report the spectral parameters and normalizations, and Figure 3 shows the fits.

Splitting this interval at $12.29$ s also yields fit improvements, but competing temperature solutions prevent a unique determination of the BB temperature evolution. In the subsequent $50$–$115$ and $115$–$164.9$ s intervals, adding BB to Band gives $\Delta\mathrm{PG}=6.88$ and $1.11$, respectively, with a linear background. We display the Band fit for the final interval, where the BB parameters are poorly constrained.

The fit improvements depend on the background prescription and the allowed BB normalization range. We therefore do not assign a calibrated detection significance to the additional component. These tests are described in Appendix A.

For the second bright episode, the $332.8$–$390$ s spectrum is described by a Comptonized model (Figure 2, right panel), with $\alpha_{\mathrm{GBM}}=-1.05\pm0.03$, $E_{\mathrm{peak}}=714^{+64}_{-55}$ keV, and $\mathrm{PG}=355.92$ for 391 fitted channels under the linear-background prescription. The parameter uncertainties are one-parameter profile intervals at $\Delta\mathrm{PG}=1$.

Replacing Comptonized with Band improves the statistic by $\Delta\mathrm{PG}=6.77$, but the high-energy index is not bounded on both sides within the tested domain, and the improvement depends on the background and detector selection. Band+BB fits also remain sensitive to the imposed normalization bound and do not provide a stable additional-component estimate. We therefore retain Comptonized as the reference description. The overlapping $334$–$390$ s extraction gives consistent parameters and serves as a sensitivity check (Appendix A).

Adding a BB improves some of the early GBM fits, but these comparisons do not uniquely identify the origin of the additional spectral curvature.

Description of the multiwavelength afterglow observations

We use external-shock models to interpret the early and late afterglow phases identified in § 3.4. We first examine the late optical and X-ray emission to guide the choice of circumburst density profile, and then consider the early X-ray and LAT emission. We connect these stages by modelling the X-ray rebrightening through energy injection and evolving shock microphysics.

Late afterglow

The late optical and X-ray observations provide temporal and spectral constraints on the afterglow interpretation. Using the convention $F_{\nu}\propto t^{-\alpha}\nu^{-\beta}$, the separate light-curve fits in § 3.4 give $\alpha_{\mathrm{opt}}=1.46\pm0.01$ and $\alpha_{\mathrm{X}}=1.57\pm0.04$. Complementary spectral information comes from the optical-to-X-ray SEDs shown in the lower panel of Figure 8. Single-power-law fits yield $\beta_{\mathrm{opt/ NIR/ X}}=0.85\pm0.01$, $0.89\pm0.01$, and $0.90\pm0.01$ at $0.18$–$0.20$, $0.85$–$1.00$, and $2.30$–$2.70$, respectively (Table 9).

We examine whether these temporal and spectral properties are compatible with synchrotron emission in which the optical and X-ray bands sample the same power-law electron segment. In the slow-cooling regime, this would correspond to either $\nu_{\mathrm{m}}^{\mathrm{syn}}<\nu_{\mathrm{opt}}<\nu_{\mathrm{X}}<\nu_{\mathrm{c}}^{\mathrm{syn}}$ or $\nu_{\mathrm{m}}^{\mathrm{syn}}<\nu_{\mathrm{c}}^{\mathrm{syn}}<\nu_{\mathrm{opt}}<\nu_{\mathrm{X}}$, where $\nu_{\mathrm{m}}^{\mathrm{syn}}$ and $\nu_{\mathrm{c}}^{\mathrm{syn}}$ denote the characteristic and cooling synchrotron frequencies, respectively.

For a power-law distribution of accelerated electrons with index $p$, the synchrotron spectral indices are $\beta=(p-1)/2$ below the cooling frequency and $\beta=p/2$ above it. If the broadband indices are identified within a single synchrotron segment, the former interpretation gives $p=2\beta+1\simeq2.7$–$2.8$, whereas the latter gives $p=2\beta\simeq1.7$–$1.8$. For $p\simeq2.7$–$2.8$, a uniform medium predicts $\alpha\simeq1.28$–$1.35$, whereas a steady stellar wind predicts $\alpha\simeq1.78$–$1.85$. Within the $p>2$ prescription considered here, this comparison motivates exploring $\nu_{\mathrm{m}}^{\mathrm{syn}}<\nu_{\mathrm{opt}}<\nu_{\mathrm{X}}<\nu_{\mathrm{c}}^{\mathrm{syn}}$. This identification is treated as a working hypothesis, given the departures from a single power law in the first two optical-to-X-ray SEDs (Table 9).

In this spectral regime, the temporal decay provides a diagnostic of the circumburst density profile, parametrized as $n(R)\propto R^{-k}$. The cases $k=0$ and $k=2$ correspond to a uniform medium and a steady stellar wind, respectively, while $0<k<2$ describes an intermediate density gradient. For an adiabatic forward shock before a jet break, the temporal and spectral indices satisfy

Figure 10 illustrates the dependence of the predicted temporal index on $p$ and $k$, with contours marking the representative optical and X-ray decay indices. Within the stated assumptions, this comparison provides a physical motivation for considering a stratified circumburst medium.

Combining the reference temporal index $\alpha_{\mathrm{ref}}=(\alpha_{\mathrm{opt}}+\alpha_{\mathrm{X}})/2$ with the observed broadband spectral indices, the closure relations (Equation (1)) support an intermediate density gradient, with $k\simeq1.0$–$1.3$. We use $k=1.2$, corresponding to $\beta\simeq0.87$ and $p\simeq2.73$, in the modelling below. Here, $\alpha_{\mathrm{ref}}$ represents the midpoint between the separately fitted temporal indices.

The observational constraints on a possible jet break and their implications for the jet geometry are discussed in § 6.6.

Early afterglow

The joint analysis of three XRT spectra extracted outside the flare interval, with a shared intrinsic absorption column, gives $\Gamma_{\mathrm{XRT}}=1.78^{+0.07}_{-0.07}$ during $263.5$–$314.0$ s (WT), $2.25^{+0.15}_{-0.14}$ during $3652.5$–$3757.0$ s (WT), and $2.01^{+0.10}_{-0.10}$ during $3763.4$–$4656.0$ s (PC), with times measured from the GBM trigger (Table 7 and Figure 7). The earliest WT spectrum is harder than the two subsequent spectra, whose central photon indices do not indicate monotonic softening. For forward-shock synchrotron emission at $\nu_{\mathrm{X}}>\max(\nu_{\mathrm{m}}^{\mathrm{syn}},\nu_{\mathrm{c}}^{\mathrm{syn}})$, the expected photon index is $\Gamma_{\mathrm{X}}=1+p/2$. The standard adiabatic closure relation before a jet break gives $\alpha_{\mathrm{X}}=(3p-2)/4$, independent of the circumburst density index $k$ in this spectral regime. For $p\simeq2.35$–$2.4$, these relations predict $\Gamma_{\mathrm{X}}\simeq2.18$–$2.20$ and $\alpha_{\mathrm{X}}\simeq1.26$–$1.30$. The spectral prediction is consistent with the resomewhat harder. The predicted temporal decay is close to that of the third segment of the catalogue XRT fit (Table 6).

The LAT observations provide an additional test of the forward-shock interpretation. The light curve over 164.9–1100 s decays with $\alpha_{\mathrm{LAT}} = 1.30^{+0.20}_{-0.19}$ (§ 3.2), while the five principal spectral intervals are compatible with a shared photon index $\Gamma_{\mathrm{LAT}} = \beta_{\mathrm{LAT}}+1 = 2.24^{+0.06}_{-0.06}$ (§ 4.2). We consider synchrotron and SSC emission above their respective characteristic and cooling frequencies, $\nu_{\mathrm{LAT}} > \max(\nu_m^i,\nu_c^i)$, where $i = \mathrm{syn}, \mathrm{SSC}$. In this spectral segment, both mechanisms predict $\beta^i = p/2$. For an adiabatic, decelerating forward shock before a jet break, the temporal indices are

The SSC relation assumes scattering in the Thomson regime. For $p \simeq2.35$–$2.4$ and $k = 1.2$, these relations predict $\alpha^{\mathrm{syn}} \simeq1.26$–$1.30$, $\alpha^{\mathrm{SSC}} \simeq1.38$–$1.43$, and $\Gamma\simeq2.18$–$2.20$ for both mechanisms, consistent with the measured LAT temporal and spectral indices.

A common electron index ($p \simeq2.35$–$2.4$) can therefore account approximately for the representative X-ray and LAT decay slopes through forward-shock synchrotron and SSC emission, respectively, above their corresponding characteristic and cooling frequencies. This interpretation does not fully reproduce the XRT spectral evolution. The closure relations alone also allow synchrotron emission in the LAT band; the highest-energy photons provide additional constraints on the emission mechanism (§ 6.1). The initial LAT excess is considered separately in § 5.4.

The X-ray rebrightening: connecting the early and late afterglow

The closure-relation comparisons in § 5.1 and § 5.2 suggest representative electron indices of $p \simeq2.73$ for the late afterglow and $p \simeq2.35$–$2.4$ for the early afterglow. Motivated by these comparisons and the intervening X-ray rebrightening, we model the connection between the two phases through energy injection and evolving shock microphysics.

Slower ejecta can catch up with the decelerating blast wave and transfer additional energy to the shocked flow through refreshed shocks ((Sari & Mészáros 2000)). Evolution of the effective electron and magnetic energy fractions could also modify the emission during this transition. Such variations are additional phenomenological assumptions, rather than a necessary consequence of energy injection ((Veres & Mészáros 2014)). Variable-microphysics models illustrate their possible effects on afterglow temporal and spectral behaviour (e.g., see Ioka et al. 2006; Fan & Piran 2006; Panaitescu et al. 2006; Lemoine 2013; Fraija et al. 2024b).

The rebrightening is shown in Figure 6 and described by the fourth segment of the Swift/XRT catalogue fit (Table 6). We adopt $t_a = 4700$ s and $t_b = 6200$ s after the GBM trigger as the boundaries of the modelled transition. We denote by $\epsilon_{e,f}$ and $\epsilon_{B,f}$ the fractions of the post-shock internal energy carried by accelerated electrons and magnetic fields, respectively, where the subscript $f$ refers to the forward shock. During $t_a < t < t_b$, we adopt

Here, $p$, $\epsilon_{e,f}$, and $\epsilon_{B,f}$ denote the pre-transition electron index and microphysical energy fractions of the early afterglow model. The indices $\alpha_e$ and $\alpha_B$ describe the temporal evolution of the electron and magnetic energy fractions, respectively. The exponent $s_p = \ln(p_{\mathrm{late}}/p)/\ln(t_b/t_a)$ connects $p$ to the adopted late-time value $p_{\mathrm{late}} = 2.75$, within the range inferred in § 5.1. For $t \leq t_a$, these quantities retain their pre-transition values; for $t \geq t_b$, they remain fixed at the values reached at $t_b$.

We denote by $E_{\mathrm{K,iso}}$ the isotropic-equivalent kinetic energy before the rebrightening and by $E_{\mathrm{K,iso}}(t)$ its subsequent evolution. Energy injection is restricted to $t_a < t < t_b$, such that

We adopt an energy increase by a factor $f_E \equiv E_{\mathrm{K,iso}}(t_b)/E_{\mathrm{K,iso}} = 2.52$, which gives $q_{\mathrm{inj}} \equiv\ln f_E/\ln(t_b/t_a) = 3.34$. This increase acts together with the prescribed microphysical evolution to connect the early and late afterglow flux levels.

A reverse-shock SSC interpretation of the early LAT excess

The interval-averaged 0.1–1000 GeV photon flux during 70.04–164.9 s exceeds the backward extrapolation of the power-law fit to the four subsequent LAT intervals by a factor of $2.25 \pm0.59$ (§ 3.2). To characterize this early excess and its steeper decline, we fitted a continuous broken power law with the break fixed at $t_{\mathrm{br}} = 164.9$ s, coinciding with the adopted end of the first bright GBM/BAT episode marked in Figure 1. The fit gives $\alpha_1 = 2.95^{+0.36}_{-0.35}$ and $\alpha_2 = 1.30^{+0.20}_{-0.19}$. Uncertainties correspond to 68.3 per cent profile intervals obtained with $\Delta\chi^2 = 1$, reoptimizing the normalization and the other temporal index while keeping $t_b = 164.9$ s fixed. The spectral models, exposure weights, and flux uncertainties were held fixed. The early index is inferred from a single interval-averaged flux under the imposed continuity and fixed break time.

This inferred early decline is steeper than the forward-shock synchrotron and SSC predictions considered in § 5.2, motivating an additional short-lived component. We describe the early LAT excess with reverse-shock SSC emission in the thick-shell regime, following the stratified-medium prescriptions of Fraija et al. (2024a). We adopt a fiducial reverse-shock crossing time of $t_\times= 150$ s after the GBM trigger. This model timescale is distinct from the fixed break time used to characterize the LAT light curve.

After the reverse shock crosses the ejecta, fresh electron injection ceases and the previously accelerated electrons continue to cool. As the SSC cutoff moves below the LAT band, the observed emission can become dominated by photons arriving from larger angles to the line of sight (Panaitescu & Kumar 2000; Kobayashi & Zhang 2003). For a power-law spectral segment, the curvature relation gives $\alpha= 2 + \beta$, with the temporal index evaluated using the appropriate time origin for the emitting episode.

The early LAT photon index, $\Gamma_{\mathrm{LAT}} = 2.21 \pm0.08$, corresponds to $\beta_{\mathrm{LAT}} = 1.21 \pm0.08$ and hence $2 + \beta_{\mathrm{LAT}} = 3.21 \pm0.08$, consistent numerically with the inferred early decay index. This numerical agreement does not constitute a test of the curvature relation, because the post-crossing decay is not separately resolved. The combined forward- and reverse-shock model is examined in § 5.5.

Constraints on the physical parameters

We use Markov chain Monte Carlo (MCMC) sampling to explore the physical parameters of the proposed external-shock interpretation. We first fit a forward-shock model to the selected LAT, X-ray, and optical/IR observations, excluding the first LAT interval to establish a reference solution for the longer-lived emission. The sampled parameters are $E_{\mathrm{K,iso}}$, $n_{\mathrm{dec}}$, $p$, $\log_{10}\Gamma_0$, $\varepsilon_{Bf}$, $\varepsilon_{ef}$, $\zeta_{ef}$, $\alpha_e$, and $\alpha_B$. Here, $\Gamma_0$ is the initial bulk Lorentz factor, and $\zeta_{ef}$ is the fraction of forward-shock electrons injected into the modelled non-thermal distribution. The parameters $\alpha_e$ and $\alpha_B$ describe the prescribed temporal evolution of the electron and magnetic energy fractions.

We then fit a combined forward- and reverse-shock model, now including the first LAT interval, with reverse-shock SSC accounting for the early excess discussed in § 5.4. This model introduces $\varepsilon_{Br}$, $\varepsilon_{er}$, and $\zeta_{er}$ as additional parameters, where the subscript $r$ denotes the reverse shock. The shared dynamical and forward-shock parameters are sampled jointly with the reverse-shock parameters.

We fix the redshift to $z = 3.51$ and calculate the luminosity distance for a spatially flat $\Lambda$CDM cosmology with $H_0 = 67.4\ \mathrm{km\ s^{-1}\ Mpc^{-1}}$ and $\Omega_M = 0.315$.

The circumburst number-density profile is parametrized as $n(R) = n_{\mathrm{dec}}(R/R_{\mathrm{dec}})^{-k}$, where $n_{\mathrm{dec}}$ is the ambient number density at the deceleration radius. We fix $k = 1.2$ and sample $n_{\mathrm{dec}}$ directly. Using the approximate deceleration criterion $\Gamma_0^2 M_{\mathrm{sw}}(R_{\mathrm{dec}})c^2 \simeq E_{\mathrm{K,iso}}$, we calculate

Here, $m_p$ is the proton mass and $c$ is the speed of light. This expression follows from $M_{\mathrm{sw}}(R_{\mathrm{dec}}) = 4\pi m_p n_{\mathrm{dec}}R_{\mathrm{dec}}^3/(3-k)$. The equivalent normalization in $n(R) = AR^{-k}$ is $A = n_{\mathrm{dec}}R_{\mathrm{dec}}^k$.

We explored the model parameters using emcee version 3.1.6 ((Foreman-Mackey et al. 2013)). The likelihood combined the Swift/XRT 0.3–10 keV, Fermi/LAT, and optical observations. Asymmetric measurement uncertainties were symmetrized by adopting the larger of the upper and lower errors for each data point. We adopted uniform priors over the ranges listed in Table 11, defined in $\log_{10}\Gamma_0$ for the initial bulk Lorentz factor and in linear space for the remaining free parameters. The walkers were initialized around the best-fitting parameter values obtained from a preliminary $\chi^2$ minimization.

Parameter	Priors	Values	Values
Shock:		Forward	Forward + Reverse
$E_{\mathrm{K,iso}}$ ($10^{55}$ erg)	$U(3.0, 9.16)$	$7.02^{+0.41}_{-0.40}$	$7.01^{+0.39}_{-0.37}$
$n_{\mathrm{dec}}$ (cm$^{-3}$)	$U(0.8, 7.2)$	$2.46^{+0.58}_{-0.57}$	$2.51^{+0.51}_{-0.46}$
$p$	$U(2.2, 2.5)$	$2.37^{+0.06}_{-0.06}$	$2.35^{+0.07}_{-0.07}$
$\log_{10}[\Gamma_0]$	$U(1, 3.1)$	$2.57^{+0.27}_{-0.28}$	$2.61^{+0.28}_{-0.29}$
$\varepsilon_{\mathrm{Bf}}$ ($10^{-5}$)	$U(1.0, 6.1)$	$1.50^{+0.34}_{-0.34}$	$1.42^{+0.39}_{-0.37}$
$\varepsilon_{\mathrm{ef}}$ ($10^{-2}$)	$U(0.5, 3.0)$	$2.06^{+0.20}_{-0.20}$	$2.11^{+0.19}_{-0.18}$
$\alpha_{\mathrm{B}}$	$U(0.92, 3.5)$	$3.00^{+0.34}_{-0.36}$	$3.00^{+0.34}_{-0.36}$
$\alpha_{\mathrm{e}}$	$U(1.14, 10.5)$	$3.89^{+0.54}_{-0.54}$	$3.89^{+0.54}_{-0.54}$
$\zeta_{\mathrm{ef}}$	$U(0.1, 1)$	$0.88^{+0.10}_{-0.10}$	$0.89^{+0.11}_{-0.11}$
$\varepsilon_{\mathrm{Br}}$ ($10^{-5}$)	$U(4.85, 12.1)$	—	$9.78^{+0.34}_{-0.34}$
$\varepsilon_{\mathrm{er}}$ ($10^{-1}$)	$U(1.5, 4.19)$	—	$2.04^{+0.21}_{-0.20}$
$\zeta_{\mathrm{er}}$ ($10^{-2}$)	$U(2.8, 4)$	—	$3.42^{+0.11}_{-0.11}$

For the forward-shock fit, we used 128 walkers and 10,000 steps per walker, discarding the first 1,000 steps as burn-in. The combined forward- and reverse-shock fit used the same sampling configuration, with the predicted flux given by the sum of the two emission components. Convergence was assessed by requiring the post-burn-in chain length to exceed 50 times the estimated integrated autocorrelation time for each sampled parameter. The mean acceptance fraction was approximately 0.3. The parameter traces were also inspected visually to check for stable sampling without persistent drifts.

Table 11 reports the posterior medians and 16th–84th percentile intervals for both fits. Figure 11 shows the posterior distributions for the forward-shock-only model. For brevity, we omit the corner plot of the combined fit, whose shared parameter estimates remain consistent with those of the forward-shock-only solution. Some posterior distributions approach the adopted prior boundaries; their quoted intervals are conditional on these bounds.

The kinetic energy and electron spectral index reported in Table 11 correspond to the pre-transition quantities $E_{\mathrm{K,iso}}$ and $p$. The tabulated forward-shock energy fractions correspond to $\varepsilon_{ef}$ and $\varepsilon_{Bf}$ in § 5.3. The sampled indices $\alpha_e$ and $\alpha_B$ govern their evolution between $t_a = 4700\ \mathrm{s}$ and $t_b = 6200\ \mathrm{s}$ after the GBM trigger. We fix $p_{\mathrm{late}} = 2.75$ and $f_E = 2.52$. After $t_b$, the kinetic energy and microphysical parameters remain constant.

Figure 12 compares the model light curves with the observations. The inset in the upper panel shows the unsmoothed XRT light curve, highlighting the rebrightening described through energy injection and evolving shock microphysics in § 5.3. The forward-shock-only solution underpredicts the first LAT interval, whereas the combined model accounts for this early excess through reverse-shock SSC emission. The inferred reverse-shock parameters are conditional on this emission mechanism and the adopted dynamical prescription.

Physical implications of the proposed interpretation

We discuss the implications of the multiwavelength observations and external-shock modelling for particle acceleration, ejecta magnetization, shock microphysics, and the circumburst environment. We also examine the prompt radiative efficiency and constraints on jet geometry, and compare the gamma-ray, X-ray, and optical emission with that of other energetic GRBs.

Maximum synchrotron energy and SSC emission

Figure 13 compares the energy–time distribution of the LAT events with the forward-shock maximum synchrotron energy and an acceleration-based reference for the reverse shock. The grey scale indicates the probability of association with GRB 260924A under the adopted LAT source and background model.

For the acceleration-based estimate, we adopt $t'_{\mathrm{acc}} = \eta_{\mathrm{acc}}2\pi r'_L/c$ with $\eta_{\mathrm{acc}} = 1$, corresponding to one gyration period. Here, $r'_L$ is the electron Larmor radius. Unless otherwise stated, photon energies and times are expressed in the observer frame; primed quantities refer to the comoving frame of the emitting fluid. Equating this timescale to the pitch-angle-averaged synchrotron cooling time gives $\gamma_{\mathrm{max},i} = [3e/(\eta_{\mathrm{acc}}\sigma_T B'_i)]^{1/2}$, where $i = f,r$ denotes the forward or reverse shock, $e$ is the elementary charge, $\sigma_T$ is the Thomson cross-section, and $B'_i$ is the comoving magnetic-field strength.

Using the photon-energy convention $E^{\mathrm{syn}}_i = heB'_i\gamma_i^2/(m_ec)$ and approximating the Doppler factor as $\delta_i \simeq\Gamma_i$, we obtain $E_{\mathrm{max},i}^{\mathrm{syn}}(t) \simeq5.57\ \mathrm{MeV}\,\Gamma_i(t)$ for $z = 3.51$ and $\eta_{\mathrm{acc}} = 1$. Here, $\Gamma_i$ is the bulk Lorentz factor of the corresponding emitting shocked fluid, $\hbar$ is the reduced Planck constant, and $m_e$ is the electron mass.

For an adiabatic, constant-energy forward shock decelerating in a medium with $n \propto R^{-k}$, this maximum energy evolves as $E_{\mathrm{max},f}^{\mathrm{syn}} \propto t^{-(1-k)/[2(4-k)]}$. For the adopted density index $k = 1.2$, this gives $E_{\mathrm{max},f}^{\mathrm{syn}} \propto t^{-9/28}$; the wind case, $k = 2$, recovers the $t^{-1/4}$ scaling. This relation applies during the decelerating, constant-energy phase. Before deceleration, an approximately constant bulk Lorentz factor gives an approximately constant acceleration-based maximum energy under the same assumptions. The LAT interval considered here precedes the prescribed energy-injection phase in § 5.3. The reverse-shock evolution must instead follow the dynamics and electron population of the shocked ejecta.

For the reverse shock, fresh electron injection ceases at the adopted crossing time $t_\times= 150\ \mathrm{s}$ (§ 5.4). Beyond this time, the blue curve serves only as an acceleration-based reference: it does not represent the evolving synchrotron cutoff of the cooling electron population or its delayed high-latitude emission.

Several high-probability LAT events exceed the synchrotron energies represented in Figure 13, including the 43.31, 24.95, and 14.24 GeV photons detected at 82.24, 122.92, and 304.15 s, respectively. Together with the light-curve modelling in § 5.5, this supports an SSC contribution from the forward and reverse shocks under the adopted prescriptions ((Sari & Esin 2001; Fraija et al. 2019b)).

The 14.24 GeV photon is particularly informative because it arrived between the two main GBM/BAT episodes. Forward-shock SSC could contribute at this epoch, while the reverse-shock contribution depends on its evolving electron distribution and SSC cutoff. The arrival time alone does not separate these possibilities or exclude a prompt contribution, since weaker activity persists between the bright episodes.

At $z = 3.51$, attenuation by the extragalactic background light (EBL) is relevant to the interpretation of the highest-energy LAT photons. To estimate this propagation effect, we use the EBL model of Saldaña-López et al. (2021) and its published optical-depth tables (Domínguez et al. 2024). The observed differential photon flux is related to the intrinsic flux by $\phi_{\mathrm{obs}}(E,t) = \phi_{\mathrm{int}}(E,t)\exp[-\tau_{\mathrm{EBL}}(E,z)]$, where $E$ is the observer-frame photon energy. Using the reference cosmology of these tables, the transmissions at 14.24, 24.95, and 43.31 GeV are approximately 0.89, 0.61, and 0.27, respectively. These transmissions characterize propagation losses and do not constitute an EBL-corrected refit of the broadband model. EBL absorption reduces the transmitted flux without changing the energies of the surviving photons, and therefore does not alter their comparison with the maximum synchrotron energy.

Klein–Nishina effects

The emission model includes Klein–Nishina effects for both the forward and reverse shocks. We evaluate the characteristic Klein–Nishina energies for both shock components near the arrival times of the highest-energy LAT photons.

For the forward shock, the observer-frame values are $E_{\mathrm{KN,c,f}}^{\mathrm{SSC}} \simeq86$ GeV and $E_{\mathrm{KN,m,f}}^{\mathrm{SSC}} \simeq252$ GeV at 80 s, and approximately 99 and 136 GeV at 300 s, respectively. The 43.31 GeV and 14.24 GeV photons therefore lie below both characteristic scales evaluated near their respective arrival times. These estimates do not indicate a strong Klein–Nishina suppression at those photon energies within the adopted forward-shock prescription.

For the reverse shock, the observer-frame Klein–Nishina characteristic energies at $t \simeq80$ s are $E_{\mathrm{KN,c,r}}^{\mathrm{SSC}} \simeq38$ GeV and $E_{\mathrm{KN,m,r}}^{\mathrm{SSC}} \simeq150$ GeV. At the adopted crossing time, $t_{\times} = 150$ s, these values are approximately 40 and 150 GeV, respectively. The 43.31 GeV photon detected at 82 s lies slightly above the lower of these scales. Klein–Nishina suppression may therefore affect its interpretation as reverse-shock SSC emission, although these characteristic energies do not represent sharp spectral cutoffs.

Ejecta magnetization

The reverse-shock SSC interpretation of the early LAT excess provides a way to examine the magnetic properties of the shocked ejecta. A magnetic field carried by the outflow could contribute to the synchrotron emission that supplies the seed photons for this SSC component.

For the combined forward- and reverse-shock fit (Table 11), the posterior median magnetic-energy fractions are $\varepsilon_{B,r} = 9.78 \times10^{-5}$ and $\varepsilon_{B,f} = 1.42 \times10^{-5}$. Their ratio is approximately 6.9, giving

This value is calculated from the reported posterior medians. Although the two components are constrained mainly at different epochs, the constant pre-transition forward-shock fraction adopted in Section 5.3 allows this comparison near shock crossing. In the hydrodynamic limit, where the post-shock internal-energy densities are comparable, $B'_r/B'_f \simeq R_B$.

The larger reverse-shock magnetic-energy fraction is compatible with a magnetic field retained by the ejecta (Zhang et al. 2003), although different efficiencies of magnetic-field amplification at the two shocks could also contribute to the contrast. This inference depends on the reverse-shock SSC interpretation of the early LAT excess.

For cold, unshocked ejecta, the magnetization is defined as $\sigma_{\mathrm{ej}} = B_{\mathrm{ej}}^{\prime2}/(4\pi\rho'_{\mathrm{ej}}c^2)$, where $B'_{\mathrm{ej}}$ and $\rho'_{\mathrm{ej}}$ are the comoving magnetic field and rest-mass density. This upstream quantity is distinct from $\varepsilon_{B,r}$ and $R_B$. Connecting them requires a shock model that accounts for compression, magnetic stresses, and energy dissipation (Zhang & Kobayashi 2005). The inferred contrast therefore supports magnetized ejecta within the adopted interpretation, but does not measure $\sigma_{\mathrm{ej}}$ or establish Poynting-flux dominance.

Electron energy partition and accelerated-electron fractions

For the combined forward- and reverse-shock fit in Table 11, the inferred electron energy fractions are $\varepsilon_{e,r} = 0.204^{+0.021}_{-0.020}$ and $\varepsilon_{e,f} = 0.0211^{+0.0019}_{-0.0018}$. The ratio of these fractions is $R_e = \varepsilon_{e,r}/\varepsilon_{e,f} \simeq9.7$. Within the adopted model, a larger fraction of the local post-shock internal energy is therefore assigned to non-thermal electrons in the shocked ejecta.

A useful comparison is provided by Fraija et al. (2024a), who modelled eight bright Fermi-LAT GRBs and obtained $2.8 \times10^{-3} \lesssim\varepsilon_{e,f} \lesssim0.12$ and $0.03 \lesssim\varepsilon_{e,r} \lesssim0.13$. Our forward-shock value lies within their reported range, while the reverse-shock value is approximately 1.6 times their largest reported value. These ranges provide modelling context rather than universal bounds. Earlier modelling of GRB 900123 also allowed substantially different microphysical parameters behind the two shocks (Panaitescu & Kumar 2004).

Using a sample of 49 radio afterglows, Duncan et al. (2023) combined radio-peak diagnostics with prompt-efficiency considerations to obtain approximate ranges $0.01 \lesssim\varepsilon_e \lesssim0.2$ and $0.1 \lesssim\zeta_e \lesssim1$. Our forward-shock values fall within these ranges. This comparison applies to the shocked ambient medium and does not impose the same bounds on the reverse-shock population.

The electron energy fractions should be interpreted together with the accelerated-electron fractions, $\zeta_r = 0.0342$ and $\zeta_f = 0.89$. Within the adopted prescription, approximately 3.4 per cent of the reverse-shock electrons carry 20.4 per cent of the local post-shock internal energy at injection. The corresponding ratios, $\varepsilon_{e,r}/\zeta_r \simeq5.96$ and $\varepsilon_{e,f}/\zeta_f \simeq2.37 \times10^{-2}$, describe how this energy is distributed among the accelerated electrons. At fixed dissipated energy per particle and electron spectral index, a smaller accelerated fraction increases the characteristic injection energy. These ratios do not directly determine $\gamma_{m,r}/\gamma_{m,f}$, since the two shocks have different strengths and injection conditions. Electrons outside the modelled power-law distribution need not be cold or radiatively inactive.

A relevant precedent is Moss et al. (2023), who adopted accelerated-electron fractions of 0.01 in the reverse shock and unity in the forward shock for GRB 030329, rather than fitting them independently. Allowing incomplete electron injection introduces degeneracies among the energy, density, and microphysical parameters in standard afterglow models (Eichler & Waxman 2005). Our inferred fractions are therefore conditional on the adopted priors and emissions sion prescriptions, including the reverse-shock SSC origin assigned to the early LAT excess.

Density profile of the circumburst medium

The temporal and spectral constraints discussed in § 5.1 favour an intermediate density gradient, $k \simeq1.0$–$1.3$. We use $k = 1.2$ within this range and hold it fixed during the parameter fitting. The combined forward- and reverse-shock fit gives $n_{\mathrm{dec}} = 2.51^{+0.51}_{-0.46}\ \mathrm{cm}^{-3}$ (Table 11), evaluated at the deceleration radius calculated for each sampled parameter set.

This stratification is intermediate between a homogeneous medium ($k = 0$) and a steady stellar wind ($k = 2$). Comparable profiles have been inferred for other GRBs. Yi et al. (2013) modelled the early optical afterglows of 19 GRBs and obtained $0.4 \leq k \leq1.4$, with a typical value near unity. Similarly, Liang et al. (2013) identified 38 onset components in a parent sample of 146 optical afterglows and inferred a characteristic profile $n \propto R^{-1}$. For GRB 130427A, Kouveliotou et al. (2013) obtained $k = 1.4 \pm0.2$ from the temporal and spectral behaviour within their synchrotron interpretation.

Physically, a profile shallower than $R^{-2}$ could reflect non-steady mass loss or modifications of the circumstellar wind structure. Interactions with the surrounding medium and the formation of a shocked-wind bubble provide possible ways to alter the density profile (van Marle et al. [94]). We therefore interpret $k = 1.2$ as an effective description of the external density over the radii probed by the observations, without identifying a unique mass-loss history or a wind-termination shock. This external profile governs the coupled forward- and reverse-shock dynamics; it does not describe the density distribution within the ejecta.

Prompt energetics and radiative efficiency

We estimate the prompt radiative efficiency using

The published prompt energy, $E_{\gamma,\mathrm{iso}} \simeq9.74 \times10^{54}\ \mathrm{erg}$ in the rest-frame 1–10,000 keV band (Gupta et al. 2026b), becomes $9.79 \times10^{54}\ \mathrm{erg}$ when rescaled to our adopted cosmology. An independent analysis of the GBM CSPEC data with RMfit, using non-overlapping time intervals, gives $E_{\gamma,\mathrm{iso}} \simeq9.31 \times10^{54}\ \mathrm{erg}$ for Band fits and $9.37 \times10^{54}\ \mathrm{erg}$ for Comptonized fits in the same rest-frame band. These estimates are approximately 4–5 per cent below the adopted value and provide an independent check of the burst’s extreme energetics. Combining the adopted value of $E_{\gamma,\mathrm{iso}} = 9.79 \times10^{54}\ \mathrm{erg}$ with $E_{\mathrm{K,iso}} = 7.01^{+0.39}_{-0.37} \times10^{55}\ \mathrm{erg}$ from the joint forward- and reverse-shock fit (Table 11), we obtain $\eta_{\gamma} = 0.123^{+0.006}_{-0.006}$. The quoted uncertainty propagates only the reported uncertainty in $E_{\mathrm{K,iso}}$, holding $E_{\gamma,\mathrm{iso}}$ fixed, and excludes model systematics.

This efficiency refers to the prompt emission and the pre-injection kinetic energy; it excludes the energy subsequently supplied through the prescription in § 5.3. The initial kinetic energy is approximately 7.2 times the prompt gamma-ray energy. Our inferred efficiencies lie within the range $0.01 \leq\eta_{\gamma} \leq0.36$ reported by Fraija et al. (2024a) for eight Fermi/LAT bursts. Within the adopted model, the extreme prompt energy is therefore accompanied by a large kinetic-energy reservoir rather than requiring an unusually high radiative efficiency.

An efficiency of this order is compatible with internal-shock calculations that allow a substantial dispersion in the Lorentz factors of the colliding shells (Guetta et al. 2001), although it does not uniquely identify the prompt dissipation mechanism. Its interpretation also depends on the adopted afterglow model and the treatment of radiative losses (Zhang et al. 2007). The forward and reverse shocks share the dynamical energy budget of the ejecta, so their kinetic energies are not added as independent reservoirs in this estimate.

The inferred prompt energy places GRB 260924A among the most energetic GRBs known, comparable in isotropic-equivalent energy to GRB 221009A, for which estimates of approximately $(1.0$–$1.2) \times10^{55}\ \mathrm{erg}$ have been reported (Lesage et al. 2023; Frederiks et al. 2023). Other exceptionally energetic events include GRB 080916C, GRB 090323, GRB 090902B, and GRB 220101A (Beniamini et al. 2015; Scotton et al. 2023). These comparisons place GRB 260924A at the extreme high-energy end of the observed population, although differences in spectral integration bands and temporal coverage preclude a precise ranking from the published values alone.

The afterglow modelling likewise implies a substantial isotropic-equivalent kinetic-energy reservoir. Estimates of several $10^{55}\ \mathrm{erg}$ have previously been obtained from the GeV afterglows of GRB 090323 and GRB 090902B under a forward-shock synchrotron interpretation (Beniamini et al. 2015). The quantitative comparison depends on the adopted radiation mechanism and microphysical parameters. Furthermore, the absence of a robust jet-break measurement in our analysis leaves the collimation correction uncertain. The extreme isotropic-equivalent energetics therefore do not, by themselves, establish an equally exceptional beaming-corrected energy.

Constraints on a possible jet break

The late-time XRT tests described in § 3.4 and Appendix § 13.1 do not establish a robust steepening. The improvement obtained by introducing a break depends on the data selection, and its epoch remains poorly constrained. The optical observations provide no evidence for a contemporaneous, sustained steepening that would support an achromatic jet break.

A jet break remains possible, but the available observations provide neither a measured break time nor a robust lower bound on it. Consequently, we do not infer a jet opening angle or a beaming-corrected energy from these data. The true energy budget therefore remains uncertain despite the extreme isotropic-equivalent energetics discussed in § 6.5.

Origin of the second bright episode

The second bright episode exhibits a pronounced GBM and BAT pulse accompanied by an XRT flare (Figure 9). Their temporal association suggests an additional emission component superposed on the underlying afterglow, although it does not establish a common emission mechanism across these bands. Internal dissipation through internal shocks or magnetic reconnection provides one possible explanation (Rees & Meszaros 1994; Zhang & Yan 2011).

The energy-dependent pulse profiles and hardness evolution provide complementary constraints. GBM and XRT show opposite average hardness trends over the covered episode, whereas the XRT hardness is consistent with remaining constant during the observed decline. These differences do not by themselves distinguish internal dissipation from an external contribution.

LAT emission is detected during 332.8–390 s, but the five-bin light curve does not resolve a distinct pulse comparable to those observed at lower energies. Within our afterglow model, the LAT emission can remain dominated by the external-shock component. The interval-averaged measurement nevertheless allows a shorter or weaker contribution associated with the second episode.

The GBM spectrum is described by an empirical Comptonized function (Table 1), which does not identify a specific radiation mechanism. Likewise, the XRT decay analysis in § 3.3 does not uniquely determine the emission origin: the fitted temporal index depends strongly on the assumed time origin, and coverage ends at 359.37 s, before the full flare decay is observed. These limitations prevent a decisive curvature-effect test ((Liang et al. 2006)).

Multiwavelength comparison with the GRB population

Figure 14–Figure 16 place GRB 260924A in the context of other GRBs across the X-ray, optical, and GeV bands. These comparisons use different observables: isotropic-equivalent rest-frame X-ray luminosities, optical magnitudes transformed to a common $z = 1$ system, and observed LAT photon fluxes.

In X-rays, GRB 260924A lies near the high-luminosity end of the Swift/XRT comparison sample (Figure 14). Table 10 quantifies its position at selected rest-frame epochs. At 3500, 7200, and 20400 s, its central luminosity estimate is the second highest among the 75, 62, and 35 events, respectively, satisfying the common-epoch coverage criteria. These nominal ranks do not establish statistically significant differences between neighbouring events, and the quoted luminosity uncertainties include statistical errors only. At a rest-frame time of 3500 s, its inferred luminosity is comparable to that of GRB 221009A: the statistical uncertainties overlap, while systematic differences between the spectral and extraction procedures have not been quantified. The nominal first-place position at 75 s occurs during the flare associated with the second bright episode and depends on the adopted mode-averaged spectral conversion; it does not characterize the underlying afterglow.

Table 10. X-ray luminosity of GRB 260924A and its nominal position within the comparison sample at selected rest-frame epochs. Luminosities refer to the absorption-corrected rest-frame 0.3–10 keV band.

$t_{\mathrm{rest}}$ (s)	Phase of GRB 260924A	$L_X$ (erg s$^{-1}$)	Rank/$N^{a}$	Brighter GRBs$^{b}$
60	Early clean interval	$5.26^{+0.60}_{-0.50} \times10^{50}$	8/191	060526, 181110A, 100302A, 140114A, 130514A, 090516A, 241030A
75	Flare$^{c}$	$6.71^{+1.45}_{-0.86} \times10^{50}$	1/206	None
820	Resumed WT coverage	$(3.52 \pm0.58) \times10^{49}$	3/70	150206A, 220101A
1000	Resumed PC coverage	$(2.27 \pm0.42) \times10^{49}$	2/69	050730
3500	Later decay	$(4.81 \pm0.89) \times10^{48}$	2/75	240529A
7200	Later decay	$(2.37 \pm0.38) \times10^{48}$	2/62	111209A
20400	Later decay	$(3.75 \pm0.89) \times10^{47}$	2/35	220101A

The optical afterglow is similarly remarkable. Figure 15 compares its evolution with the long-GRB sample of Kann et al. (2006, 2010, 2011, 2024) and with GRB 221009A using the observations of O’Connor et al. (2023). After transforming the light curves to $R$-band AB magnitudes in a common $z = 1$ system, GRB 260924A lies near the upper envelope of the sample out to approximately one day. In this reference frame, it is approximately 2 mag brighter than GRB 221009A at 0.3 d and 1 mag brighter near 1 d, corresponding to flux ratios of approximately 6.3 and 2.5. These differences depend on the extinction and spectral corrections applied to both events.

In the GeV band, GRB 260924A lies towards the bright end of the selected LAT comparison sample over approximately 70–1100 s, particularly during its first observed interval (Figure 16). The comparison uses observer-frame times and photon fluxes integrated above 100 MeV, with dataset-dependent upper energy limits of 100, 300, or 1000 GeV. For GRB 260924A, we use the 0.1–1000 GeV band. It therefore characterizes apparent brightness and temporal evolution without establishing a ranking in intrinsic GeV luminosity.

Together with the prompt energetics discussed in § 6.5, these comparisons show that GRB 260924A combines an extreme isotropic-equivalent gamma-ray energy with luminous X-ray and optical afterglows and prominent observed GeV emission. The broadband prominence need not arise from a single radiation mechanism: the afterglow luminosities also depend on the circumburst environment and shock microphysics.

Conclusions

Our multiwavelength study places GRB 260924A among the most energetic long GRBs observed. At $z = 3.51$, its isotropic-equivalent prompt gamma-ray energy is $E_{\gamma,\mathrm{iso}} \simeq9.79 \times10^{54}$ erg. The prompt emission comprises two bright episodes, the second beginning near 333 s, separated by weaker, energy-dependent activity. BAT and LAT detect emission between the episodes, including a 14.2 GeV photon at 304 s. The highest-energy photon associated with the burst, at 43.3 GeV, arrived at 82 s during the first episode.

LAT emission extends to approximately 1100 s. The four intervals after 164.9 s follow a decay with $\alpha_{\mathrm{LAT}} = 1.30^{+0.20}_{-0.18}$, whose backward extrapolation underpredicts the first-interval average by a factor of $2.25 \pm0.59$. The five interval-resolved spectra are compatible with a shared photon index, $\Gamma_{\mathrm{LAT}} = 2.24 \pm0.06$. The second bright GBM/BAT episode is accompanied by an X-ray flare, whereas the LAT light curve does not resolve a corresponding pulse. The incomplete XRT coverage of the flare decay limits constraints on its physical origin.

The X-ray and optical afterglows lie near the high-luminosity end of their respective comparison samples. In the common $z = 1$ system, the optical afterglow is approximately 2 mag brighter than GRB 221009A at 0.3 d and 1 mag brighter near 1 d under the adopted corrections. Its rest-frame 0.3–10 keV X-ray luminosity is among the highest at the selected post-flare epochs and is comparable to that of GRB 221009A at a rest-frame time of 3500 s. The LAT comparison likewise shows prominent observed GeV emission, without establishing an intrinsic luminosity ranking.

Under the single-synchrotron-segment working hypothesis, the late afterglow observations favour a stratified circumburst medium, which we model as $n(R) \propto R^{-1.2}$. We describe the broadband light curves through forward-shock synchrotron and synchrotron self-Compton emission, with reverse-shock SSC accounting for the early LAT excess. The highest-energy LAT photons exceed the predicted synchrotron energies under the adopted acceleration prescriptions, further supporting an SSC contribution. The later X-ray rebrightening is described through energy injection from trailing ejecta and evolution of the electron spectral index and microphysical energy fractions during 4700–6200 s. The kinetic energy and these parameters remain constant afterwards. This light-curve description does not resolve the early X-ray spectral discrepancy.

The inferred injection isotropic-equivalent kinetic energy, $E_{\mathrm{K,iso}} \simeq7.0 \times10^{55}$ erg, corresponds to a prompt radiative efficiency of $\eta_{\gamma} \simeq0.12$. Within the adopted model, the extreme prompt energy is accompanied by a large kinetic-energy reservoir without requiring an unusually high radiative efficiency. The magnetic contrast $R_B \simeq2.6$ inferred from the combined shock model is compatible with a magnetic field retained by the ejecta, but does not directly measure their upstream magnetization. The available optical and X-ray observations do not establish a robust jet-break epoch, leaving the jet opening angle and beaming-corrected energy undetermined.

GRB 260924A thus provides a valuable case for examining how a large relativistic energy reservoir, the circumburst environment, and shock microphysics jointly shape the broadband afterglow.

Acknowledgements

NF is grateful to UNAM-DGAPA-PAPIIT for the funding provided by grant IN112525. MGD acknowledges the support of the DoS and by JSPS Grant-in-Aid for Scientific Research (KAKENHI) (A), Grant Number JP25H00675.

DATA AVAILABILITY

The Fermi/GBM and LAT and Swift/BAT, XRT, and UVOT observations analysed in this work are available from the respective public mission archives. The derived light curves, interval definitions, fit summaries, and scripts used to produce the figures are available from the corresponding author upon reasonable request.

Times are relative to the GBM trigger. The integrated first-episode spectrum and its subintervals use overlapping data and are not independent measurements. Reported uncertainties correspond to nominal 68.27 per cent one-parameter profile intervals at $\Delta PG = 1$, conditional on the adopted model, background, and channel selection. $PG = 2S$, where $S = -\ln L$ up to data-dependent constants. $N_{\mathrm{ch}}$ and $N_{\mathrm{src}}$ give the numbers of fitted channels and free source parameters, respectively; background nuisance parameters are excluded from $N_{\mathrm{src}}$. These counts do not define calibrated degrees of freedom.

$A$ is the continuum prefactor at the fixed 100 keV pivot. For the BB prescription defined in Appendix A5, $A_{\mathrm{BB}} = 10^{-6}A_6\ \mathrm{ph\ cm^{-2}\ s^{-1}\ keV^{-3}}$. Normalizations refer to the indicated joint fit. Time and uncertainty conventions, including superscript $a$, follow Table 1. For Band+BB, $\Delta PG = PG_{\mathrm{Band}} - PG_{\mathrm{Band+BB}}$; for the second-episode Band alternative, $\Delta PG = PG_{\mathrm{Comp}} - PG_{\mathrm{Band}}$. These differences measure fit improvements, not detection significances.

Times are relative to the GBM trigger. $F_{\mathrm{ph}}$ and $F_E$ are interval-averaged photon and energy fluxes, respectively. Errors are nominal $1\sigma$ uncertainties from the local covariance matrix; energy-flux errors include normalization–index covariance. Both diffuse normalizations are free within $[0, 30]$ under the source and background-only hypotheses. Sensitivity tests are described in Appendix B.

Notes. The catalogue excludes 115–127 and 311–3647 s from the temporal fit; these windows are not measurements of flare duration. The fourth break lies within a coverage gap and is inferred from the fitted segments. The GBM reference is obtained using $t_{\mathrm{GBM}}=t_{\mathrm{BAT}}+4.71$ s.

$^{a}$ The catalogue displays the uncertainty on $\alpha_4$ as (+0.00, +2.15). We retain the best-fitting value without interpreting these entries as a conventional two-sided confidence interval.

Notes. The spectra were fitted jointly over 0.3–10 keV with TBabszTBabspowerlaw. The Galactic column was fixed at $N_{\mathrm{H,Gal}}=4.09\times10^{20}\,\mathrm{cm}^{-2}$. The shared intrinsic column at $z=3.51$ is $N_{\mathrm{H,int}}=(7.46^{+5.17}_{-4.88})\times10^{21}\,\mathrm{cm}^{-2}$; photon indices and normalizations vary independently. The joint fit gives $W=896.40$ for 1009 degrees of freedom. Errors correspond to one-parameter 90 per cent profile intervals, $\Delta W=2.71$. Flux-density uncertainties follow from the profiled power-law normalizations. Time intervals delimit the selected GTIs; exposures give the live integration times.

Additional details of the GBM spectral analysis

This appendix provides additional details of the response construction, background estimation, exposure calculation, and sensitivity tests used in the Fermi/GBM spectral analysis described in § 2.1.

Effective response construction

For each detector, we constructed an effective response by averaging the time-dependent RSP2 response over the corresponding spectral extraction interval. The response at the centre of each temporal bin was obtained by linear interpolation between the temporal centres of adjacent response matrices and weighted by the total observed counts in that bin, including both source and background counts.

For the first-episode spectra covering −6.1–50, −6.1–164.9, 50–115, and 115–164.9 s, the response weights were derived from the total observed counts in each detector’s fitted energy range, using 1.024-s bins. The weighting included the full bins intersecting each extraction interval, so their coverage could extend slightly beyond the spectral boundaries. For the two subdivisions of the −6.1–50 s interval, the weighting bins were instead clipped to the extraction boundaries.

For the second bright episode, the effective response was constructed separately for each detector. The weighting grid had a nominal spacing of 1.024 s and included the exact extraction boundaries and the response transition at $T_0 + 366.597954$ s. Counts were extracted directly from the TTE data using the full detector energy range adopted before applying the channel exclusions used in the principal spectral fit. This procedure combined RSP5 and RSP6 using detector-specific weights.

As a sensitivity test, we replaced the total-count weighting with weighting by dead-time-corrected exposure. We also compared the linearly interpolated responses with piecewise responses assigned according to their nominal validity intervals. The resulting spectral parameters were consistent with those obtained using the adopted response prescription.

Background estimation and uncertainty propagation

For the first-episode spectra, the background was fitted independently in each detector channel using first- and second-order polynomials on 1.024-s temporal bins over the nominal off-source intervals [−125, −35] and [250, 320] s relative to $T_0$. The fits used the two-pass weighted least-squares procedure implemented in the Gamma-ray Data Tools (GDT). Interval-averaged background rates and their uncertainties were computed from the fitted coefficients and their full covariance matrix.

The same two time-resolved subdivisions of the initial episode used the same off-source intervals. Their background estimates therefore share the same off-source data and have correlated background uncertainties. These correlations are relevant when comparing the two subintervals, although each spectrum was fitted independently.

For the second bright episode (332.8–390 s), the off-source intervals were 176–200, 224–240, and 410–476 s for the NaI detectors, and 190–240 and 410–476 s for b0, relative to $T_0$. The principal analysis used linear background models fitted independently in each channel by maximizing the Poisson likelihood on 2-s bins. Positivity was imposed on a grid spanning 176–480 s; this grid was not used as a continuous off-source interval. The coefficient covariance matrix was estimated from the inverse Fisher information matrix and propagated to the integrated background counts, accounting for the live exposure.

The principal spectral fit retained 391 channels: 119 from n3, 119 from n4, 90 from n5, and 63 from b0, spanning 8.73–884.51, 8.82–899.28, 52.19–858.48, and 318.35–4880.69 keV, respectively, with no internal channel exclusions within these ranges.

The resulting background estimates and variances supplied the Gaussian background constraints used in the spectral likelihood. The temporal polynomial coefficients were not refitted during spectral fitting. Correlations between detector channels and uncertainty associated with the choice of polynomial order were not included in these variances. We therefore additionally examined alternative polynomial orders as part of the sensitivity tests described in Section A6.

Exposure calculation and time-resolved extraction

For the −6.1–50, −6.1–164.9, 50–115, and 115–164.9 s spectra, we used the detector-specific live exposures returned by GDT in the extracted PHA products. These exposures were also used to convert the estimated background rates and their uncertainties into counts.

For the diagnostic subdivision of the initial bright episode, a single temporal boundary at $T_0 + 12.288$ s was selected from the 0.512-s GBM total-count light curve before fitting the two spectra. The boundary was chosen by visual inspection to separate the initial bright double-peaked complex and its rapid decline near 12–13 s from the subsequent weaker emission. No alternative cuts were tested. The quoted precision identifies the selected bin edge, rather than an uncertainty on the physical end of the pulse. Counts and live exposures were extracted independently for the two resulting intervals, while retaining the same spectral channel selection and off-source background intervals used for the integrated spectrum.

For each subinterval, the live exposure was calculated from the extraction duration after accounting for detector dead time from all channels. This convention differs by at most 0.163 per cent from the exposure recorded in the archived PHA products and has a negligible effect on the spectral results.

Poisson–Gaussian likelihood

The auxiliary background parameter was profiled analytically in each channel without a non-negativity constraint. The total Poisson mean was constrained to be non-negative and strictly positive in channels containing observed counts.

Blackbody implementation

We describe the BB component with a Planck photon spectrum (Ryde 2004; Ryde & Pe’er 2009)

where $E$ is the photon energy, $T$ is the blackbody temperature, and $k_{\rm B}$ is Boltzmann’s constant. Both $E$ and $k_{\rm B}T$ are expressed in keV, and $N_{\rm BB}(E)$ is the differential photon flux in ph cm$^{-2}$ s$^{-1}$ keV$^{-1}$. For the numerical implementation, we express its normalization as $A_{\rm BB} = 10^{-6} A_6$ ph cm$^{-2}$ s$^{-1}$ keV$^{-3}$, where $A_6$ is dimensionless. We imposed $A_6 \ge0$ and restricted the temperature to $5 \le k_{\rm B}T \le300$ keV. To test sensitivity to the allowed normalization range, we repeated the fits using upper bounds of $A_6 = 20$, 40, and 80.

Sensitivity tests

In addition to the response tests described in Appendix A1, we examined the sensitivity of the spectral results to linear and quadratic temporal background models and to the detector-channel selections used in the fits.

For the second bright episode, particular attention was given to the stability of the results under alternative background descriptions and channel selections because of the weaker signal and the more complex temporal environment surrounding this episode. The $T_0 + 250$ to $T_0 + 320$ s interval was not used to estimate the background for this episode.

As an additional background check, we repeated the Band and Band+BB fits over 62–68 s, excluding 250–320 s from the background selection. The alternative windows were [−125, −35] and [176, 200] s for NaI, and [−125, −35] and [190, 200] s for BGO, relative to the GBM trigger. The fit improvement persisted, with $\Delta PG = 27.13$ and 34.25 for linear and quadratic backgrounds, respectively, and $k_\mathrm{B}T \simeq41$–42 keV. This supports the persistence of the candidate under the tested background selections, although temporal background residuals remain.

Additional details of the LAT likelihood analysis

This appendix provides additional details of the likelihood implementation, source-model configuration, uncertainty propagation, and sensitivity tests used in the Fermi/LAT analysis described in § 2.2.

Likelihood implementation and source model

The analysis used the UnbinnedObs and UnbinnedAnalysis interfaces with the NewMinuit optimizer and the source configuration described in § 2.2. The GRB spectrum was implemented with PowerLaw2 over 0.1–1000 GeV. The archived model-building script specified the catalogue file gll_psc_v35.fit.⁶

Test statistic and diffuse-normalization bounds

We calculated $\mathrm{TS} = 2(\ln\hat{\mathcal{L}}_1 - \ln\hat{\mathcal{L}}_0)$ for models with and without the GRB, reoptimizing both diffuse normalizations within [0, 30] under each hypothesis. Sensitivity to these bounds is examined in Appendix B4.

The reported photon-flux and photon-index errors are nominal $1\sigma$ covariance uncertainties from the source-present fit; they exclude instrumental and background-model systematics.

Energy-flux calculation and uncertainty propagation

Energy fluxes were obtained by integrating the archived PowerLaw2 spectra over 0.1–1000 GeV. Flux uncertainties were propagated from the full-fit covariance matrix, retaining the normalization–index covariance with the diffuse normalizations free during fitting.

For Figure 4, photon and energy fluxes were also evaluated over 0.1–300 GeV using the same covariance treatment. This display conversion did not modify the original 0.1–1000 GeV likelihood fits.

Sensitivity to the background model

We examined the sensitivity of the LAT results to assumptions entering the background model and likelihood optimization. In particular, we repeated the relevant likelihood calculations using wider bounds on the Galactic and isotropic diffuse normalizations. This was especially important for the first principal interval, for which the inferred source significance showed the greatest sensitivity to the allowed diffuse-normalization range.

The catalogue-source configuration, event selection, GRB spectral model, and energy range were otherwise kept unchanged when evaluating the effect of the diffuse-normalization bounds. The resulting variations were used as a diagnostic of background-model sensitivity rather than incorporated into the formal statistical uncertainties quoted for the principal fits.

Accordingly, the reported TS values and covariance-based parameter uncertainties should be interpreted as conditional on the adopted source and diffuse-background model. No additional systematic uncertainty associated with alternative Galactic diffuse templates, instrumental response functions, or catalogue-source configurations was included.

For the five integrated intervals, the principal [0, 30] diffuse-bound TS values are 1058.36, 282.41, 63.55, 309.28, and 26.95. Under wider bounds, the corresponding stable values are 863.90, 281.62, 63.55, 309.27, and 26.95. The primary table retains the common [0, 30] convention. The larger differences in the first interval may reflect the improved source-removed hypothesis under wider bounds. Stability against increasing a parameter bound is not a test of the adequacy of the background templates.

Updated temporal fit with fixed normalized exposure weights

We fitted the final interval-averaged photon fluxes using independent Gaussian errors. For each time cell $j$ in interval $i$, the archived integrated response $w_{ij}$ defines a normalized weight. The prediction is

The response is treated as constant within each cell. The primary calculation retains the archived refined weights for the first four intervals and the 16-cell gtlexposure representation for the last interval. All cell boundaries match the final good-time intervals. The final spectral indices differ from those used to generate the weights by less than $9 \times10^{-5}$ in the first four intervals and by 0.00436 in the last. We retain the original normalized weights rather than claiming exposures regenerated at the final indices. Their sensitivity to this small index change has not been separately propagated.

Replacing the last-interval weights by the archived 16-cell $N_\mathrm{pred}/F$ weights changes the fitted slope from 1.2980 to 1.2969 and the first-bin excess ratio from 2.2487 to 2.2512. This comparison tests the effect of the previously audited difference in normalized temporal shape; it does not resolve the absolute response discrepancy between the two methods. The normalization $A$ is fitted freely. The first bin is excluded from fitting and used only to assess the extrapolation; the 250–320 s diagnostic and overlapping subdivisions are not additional independent measurements. The second-pulse interval remains included, so this fit does not isolate a pure afterglow contribution.

The covariance uses the absolute quoted flux errors without rescaling by reduced $\chi^2$. The asymmetric errors on $A$ and $\alpha$ are one-parameter $\Delta\chi^2 = 1$ profile intervals of this Gaussian flux fit, not LAT event-likelihood profiles. The extrapolation and excess-ratio errors use local covariance propagation. Instrumental, spectral, background-model, and exposure-weight uncertainties are not included. Analytical cell averages were checked against independent numerical quadrature.

Additional details of the BAT spectral analysis

As a reproducibility check, we subsequently applied the explicitly identified calibration file swbsyserr20030101v003.fits (version 3) to copies of the three fitted spectra using batphasyserr version 1.6 under HEASoft 6.37.1. The resulting SYS_ERR values matched those preserved in the fitted spectra exactly in all 80 channels of each spectrum.

This comparison provides a retrospective verification of the systematic uncertainties applied in the spectral analysis. The calibration file selected during the original processing was not recorded; consequently, the identification of swbsyserr20030101v003.fits refers to the calibration file used for this subsequent reproducibility check rather than to a contemporaneous record of the original processing.

For the 70–165 s and 334–390 s intervals, refining the response construction from segments of at most 5 s to at most 2.5 s changed the fitted photon indices and energy fluxes by substantially less than their statistical uncertainties. The 250–320 s interval retained a single response and was not subdivided in this test.

Additional details of the optical/NIR and XRT analysis

The UV, optical, and NIR photometry was placed on the AB magnitude system and corrected for Galactic foreground extinction adopting $E(B - V) = 0.041$ and $R_V = 3.1$. For each photometric band, the extinction correction was evaluated at the adopted effective or pivot wavelength of the corresponding filter using the extinction prescriptions described in § 2.5.

The integrated 0.3–10 keV XRT energy fluxes were converted to flux densities at 1 keV using power-law spectra with $\Gamma_{\mathrm{WT}} = 1.49$ and $\Gamma_{\mathrm{PC}} = 2.04$, normalized to the measured band-integrated flux in each mode.

For visualization, the densely sampled XRT light curve was rebinned in time. At $t < 0.01$ d, consecutive groups of measurements were combined, while at later times the observations were grouped into predefined temporal intervals. Within each bin, the flux was calculated using inverse-variance weighting, and the representative time was determined from the error-weighted geometric mean of the individual measurement times. The horizontal uncertainties span the temporal coverage of the measurements contributing to each bin. This rebinning was used only for visualization; all temporal fits were performed using the original XRT measurements.

For the SED reconstruction, the geometric mean of the boundaries of each temporal interval was adopted as the representative epoch. Optical/NIR measurements were shifted to the corresponding epoch using the adopted optical temporal decay, with its shared uncertainty propagated into the covariance. The reconstructed SEDs include NOT $r$ and $z$ measurements in the first epoch, GROWTH $r'$ and SYSU $J$ in the second, and COLIBRÍ $r$ and $z$ together with MDM $J$ in the third. Photometric upper limits were displayed but were not included in the spectral fits.

Because an XRT measurement was not necessarily available at each representative epoch, the 1 keV X-ray flux density was evaluated at that epoch using the best-fitting late-time XRT temporal power law. This procedure places the X-ray and optical/NIR measurements at a common characteristic epoch rather than averaging the X-ray emission over a different temporal interval.

Test for a late-time X-ray break

We used the public Swift/XRT PC-mode light curves provided by UKSSDC (Evans et al. (2019, 2010)). This independent check uses the native unabsorbed 0.3–10 keV light curves, whose count-to-flux conversion tracks the hardness ratio, rather than the fixed-index 1-keV reconstruction described above. The time origin is the adopted GBM reference, $t_{\mathrm{GBM}} = t_{\mathrm{Swift}} + 4.71$ s. We selected detection bins entirely later than $10^4$ s, with no upper time cut. The selections excluding and including exposures flagged for unreliable centroiding contain 69 and 72 detections, respectively, with bin edges spanning $1.49 \times10^4$–$2.94 \times10^5$ s and $1.49 \times10^4$–$8.80 \times10^5$ s. The products regroup surviving events, so their difference is not equivalent to removing three bins. The flags indicate uncertain centroiding, not necessarily erroneous measurements; we therefore report both selections.

We compared a single power law with a continuous broken power law having free normalization, pre- and post-break slopes, and break time. Both models used the normalization time $t_{\mathrm{ref}} = 10^5$ s. We examined both unrestricted slope changes and the steepening alternative $\alpha_2 \ge\alpha_1$; Table D1 and Figure D1 show the latter. We searched $t_{\mathrm{br}}$ between the first and last detection centres using profiles and multiple starting points. No break time or slope was imposed from another waveband or from the physical afterglow model.

Product	Flagged exposures	$N$	$\alpha_{\rm PL}$	$\chi^2_{\rm PL}$/dof	$\chi^2_{\rm BPL}$/dof	$\Delta\chi^2$
$F_{0.3-10\,\mathrm{keV}}$	Excluded	69	$1.55^{+0.04}_{-0.04}$	66.32/67	66.15/65	0.17
$F_{0.3-10\,\mathrm{keV}}$	Included	72	$1.61^{+0.03}_{-0.03}$	72.43/70	67.02/68	5.41
Count rate	Excluded	69	$1.56^{+0.04}_{-0.04}$	62.12/67	62.09/65	0.03
Count rate	Included	72	$1.61^{+0.03}_{-0.03}$	66.42/70	62.64/68	3.77

Predictions were averaged over the intersections of each bin with its actual good-time intervals (GTIs), weighting by GTI duration. Fits minimized squared residuals in linear flux, using the upper-reported uncertainty when the model exceeded the datum and the lower uncertainty otherwise. Quoted slope errors correspond to nominal profile intervals with $\Delta\chi^2 = 1$. Count-to-flux conversion uncertainties and covariances are not included. Controls in count rate and in the native observed 1-keV flux density, as well as alternative bin-averaging prescriptions, gave the same qualitative conclusion. These products are correlated and their statistics were not combined.

The unflagged band-flux selection provides essentially no improvement for a steepening, and its entire searched break-time range lies within $\Delta\chi^2 = 1$ of the minimum. Including flagged exposures yields a broad minimum near $t_{\mathrm{br}} = 2.26$ d, with $\alpha_1 = 1.54$ and $\alpha_2 = 1.99$; the nominal profile region $\Delta\chi^2 \le1$ extends approximately from 1.66 to 7.45 d. This is not a robust measurement of a break epoch. Allowing either steepening or flattening does not establish a late break in the unflagged data either. The unflagged count-rate upper limit near $8.59 \times10^5$ s is also consistent with the single-power-law extrapolation and was not treated as a Gaussian detection.

We do not assign a detection significance to $\Delta\chi^2$, because the break time is unidentified under the null hypothesis and zero steepening is a parameter boundary (Protassov et al. (2002)). The hardness ratios do not establish spectral constancy across the candidate break. These data neither require a jet break nor provide a robust break epoch for the physical model.