Spectral index-flux relation for investigating the origins of steep decay in γ-ray bursts

γ-ray bursts (GRBs) are short-lived transients releasing a large amount of energy (1051 − 1053 erg) in the keV-MeV energy range. GRBs are thought to originate from internal dissipation of the energy carried by ultra-relativistic jets launched by the remnant of a massive star’s death or a compact binary coalescence. While thousands of GRBs have been observed over the last thirty years, we still have an incomplete understanding of where and how the radiation is generated in the jet. Here we show a relation between the spectral index and the flux found by investigating the X-ray tails of bright GRB pulses via time-resolved spectral analysis. This relation is incompatible with the long standing scenario which invokes the delayed arrival of photons from high-latitude parts of the jet. While the alternative scenarios cannot be firmly excluded, the adiabatic cooling of the emitting particles is the most plausible explanation for the discovered relation, suggesting a proton-synchrotron origin of the GRB emission.

T he prompt emission of GRBs is characterized by an erratic superposition of several pulses, whose spectrum typically peaks in the keV-MeV energies. Its physical origin is still matter of discussion and the main open questions concern the composition of the jet (matter-1 or magnetic 2 -dominated), the energy dissipation mechanisms (sub-photospheric emission 3 , internal shocks 4 or magnetic reconnection 5 ), and the nature of particle radiation. Once the prompt emission ceases, the light curve usually presents a steep decay (SD) phase [6][7][8][9] (tail), which can be well monitored in the X-ray band. The duration of the SD is around 10 2 -10 3 s and it is characterized by a typical decay power-law slope of 3-5. After the prompt emission, the jet interacts with the interstellar medium, producing the so called afterglow emission [10][11][12] . The afterglow models cannot account for such steep slopes and the origin of the SD is attributed to the fade-off of the emission mechanism that is responsible of the prompt phase.
Considering that the emitting surface of the jet is curved, an on-axis observer first receives photons from the line of sight and later photons from higher latitudes [13][14][15] , which are less Doppler boosted. This gives rise to the so called high-latitude emission (HLE). Under the assumption of a single power-law spectrum (F ν ∝ ν −β ), the HLE predicts that the flux decays as F ν ðt obs Þ / ν Àβ t Àðβþ2Þ obs . On the other hand, if the spectrum is curved, the HLE can also lead to the transition of the spectral peak across the observing band 16 , causing a spectral evolution, as often observed in the soft X-rays 17,18 .
In this work, we find a unique relation between the spectral index and the flux. Here, we systematically analyze the X-ray spectral evolution during the SD phase as motivated by fact that temporal and spectral evolution during the tail of prompt pulses can provide clues about emission and cooling processes in GRB jets. Given the same trend followed by all the GRBs of our sample, we search for a common process at the basis of the spectral relation. We find that the standard HLE model cannot account for the observed relation, implying that efficient cooling of particles is disfavored. We test several assumptions about the dominating cooling mechanisms and we find that the combined action of adiabatic cooling of particles and magnetic field decay robustly reproduces our data. We conclude discussing the implications for the physics of GRB jets, their composition, and radiation mechanisms.

Results
In order to investigate the spectral evolution during the SD phase, we select a sample of GRBs from the archive of the X-ray Telescope (XRT, 0.3-10 keV) on-board the Neil Gehrels Swift Observatory (Swift) 19 . We restrict our study to a sample of GRBs (eight in total) whose brightest pulse in the Burst Alert Telescope (BAT, 15-350 keV) corresponds to the XRT peak preceding the X-ray tail (see as example the Fig. 1a). We perform a timeresolved spectral analysis of the tail in the E = (0.5-10) keV band assuming a simple power-law model for the photon spectrum N γ ∝ E −α (see "Methods:" "Sample selection," "Time-resolved spectral analysis," and "Spectral modeling"). We represent the spectral evolution plotting the photon index α as a function of the flux F integrated in the E = (0.5-10) keV band, hereafter referred to as the α − F relation. The flux is normalized to the peak value of the X-ray tail. This normalization makes the result independent of the intrinsic brightness of the pulse and of the distance of the GRB.
We find a unique α − F relation for the analyzed GRBs as shown in Fig. 1b. This is consistent with a systematic softening of the spectrum; the photon index evolves from a value of α~0.5-1 at the peak of the XRT pulse to α~2-2.5 at the end of the tail emission, while the flux drops by two orders of magnitude. The initial and final photon indices are consistent with the typical low-and high-energy values found from the analysis of the prompt emission spectrum of GRBs, namely~1 and~2.3 [20][21][22] , respectively. The α − F relation can be interpreted as being due to a spectral evolution in which the spectral shape does not vary in time, but the whole spectrum is gradually shifted toward lower energies while becoming progressively dimmer (see Fig. 2). The consistent spectral evolution discovered in our analysis is a clear indication of a common physical mechanism responsible for the tail emission of GRBs and the corresponding spectral softening.
Testing HLE. We first compare our results with the expectations from the HLE, which is the widely adopted model for interpreting the X-ray tails of GRBs. When the emission from a curved surface is switched off, an observer receives photons from increasing latitudes with respect to the line of sight. The higher the latitude, the lower the Doppler factor, resulting in a shift toward lower energies of the spectrum in the observer frame. Through an b a Fig. 1 The steep decay phase and the correspondent spectral evolution. In a we show an example of a light curve of an X-ray tail selected from our sample, taken from the GRB 161117A. We show on the same plot the XRT (orange) and the BAT (blue) flux density at 1 and 50 keV, respectively. The XRT light curve decays less steeply than BAT because of the evolution of the peak energy. The error bars represent 1σ uncertainties and they are derived from the Swift archive. In b we report the spectral evolution of the X-ray tail for all the GRBs in the first sample (shown with different colors). The photon index α is represented as a function of the reciprocal of the normalized flux F max =F. Time flows from left to right. The error bars represent 1σ uncertainties, calculated via spectral fitting in XSPEC. In the legend, we report the name of each GRB.
accurate modeling of HLE (as described in "Methods:" "HLE from infinitesimal duration pulse") we derive the predicted α − F relation under the assumption of an abrupt shutdown of the emission, consistent with particles cooling on timescales much smaller than the dynamical timescale. We first consider a smoothly broken power-law (SBPL) comoving spectrum. Regardless of the choice of the peak energy, the bulk Lorentz factor or the radius of the emitting surface, the HLE predicts an α − F relation whose rise is shallower than the observed one (Fig. 3). We, additionally, test the Band function, commonly adopted for GRB spectra 23 , and the physically motivated synchrotron spectrum 4 , obtaining similar results (Supplementary Figs. 11 and 12a): the HLE softening is too slow to account for the observed α − F relation. We further relax the assumption of an infinitesimal duration pulse, i.e., considering a shell that is continuously emitting during its expansion and suddenly switches off at radius R 0 24 (see Supplementary Note 1). The contributions from regions R < R 0 are subdominant with respect to the emission coming from the last emitting surface at R = R 0 , resulting in a spectral evolution still incompatible with the observations (Supplementary Fig. 1). An interesting alternative is the HLE emission from an accelerating region 25 taking place in some Poynting flux dissipation scenarios 5 . Even though it can explain the temporal slopes observed in the X-ray tails, also this scenario fails in reproducing the α − F relation (see Supplementary Fig. 2). Our results on HLE are based on the assumption of a common comoving spectrum along the entire jet core. Even changing the curvature (or sharpness) of the spectrum or assuming a latitude dependence of the spectral shape, the disagreement with the data remains, unless we adopt a very fine-tuned structure of the spectrum along the jet core, which is not physically motivated (see "Methods:" "HLE from infinitesimal duration pulse"). Alternative models, such as anisotropic jet core [26][27][28] or subphotospheric dissipation 29 , can hardly reproduce our results (see Supplementary Note 3).
Adiabatic cooling. Since the standard HLE from efficiently cooled particles and its modified versions, as well as alternative scenarios, are not able to robustly interpret the observed α − F relation, we consider a mechanism based on an intrinsic evolution of the comoving spectrum. The most natural process is the adiabatic cooling of the emitting particles 30 . Here we assume conservation of the entropy of the emitting system hγi 3 V 0 throughout its dynamical evolution, where 〈γ〉 is the average random Lorentz factor of the emitting particles and V 0 / R 2 ΔR 0 the comoving volume 31 . We consider both thick and thin emitting regions, i.e., a comoving thickness of the emitting shell ΔR 0 ¼ const or ΔR 0 / R, respectively. We assume a power-law radial decay of the magnetic field B ¼ B 0 ðR=R 0 Þ Àλ , with λ > 0, and synchrotron radiation as the dominant emission mechanism. Here R 0 is the radius at which adiabatic cooling starts to dominate the evolution of the emitting particles. We compute the observed emission taking also into account the effect of HLE by integrating the comoving intensity along the equal arrival time surfaces (EATS) (see "Methods:" "Adiabatic cooling"). In this scenario, contrary to HLE alone, the emission from the jet is not switched off suddenly, but the drop in flux and the spectral evolution is produced by a gradual fading and softening of the source, driven by adiabatic cooling of particles. The resulting spectral evolution and light curves are shown in Fig. 4.
Adiabatic cooling produces a much faster softening of α as a function of the flux decay, with respect to HLE alone, in agreement with the data. Assuming a different evolution of the shell thickness, the behavior of the curves changes only marginally (see Supplementary Fig. 3). For large values of λ, the evolution of α flattens in the late part of the decay (see Fig. 4a), Fig. 2 Illustration of the spectral evolution caused by a shift of the spectrum towards lower energies. The transition of the spectral peak through the XRT band explains the observed spectral softening. The spectra in (b) colored in blue, green, and red correspond to the three temporal bins shown in (a) with the same colors. The inset in (b) shows how the local spectral slope evolves as observed in the XRT band. Since in the b we plot the flux density, the local slope in the XRT band is given by 1 − α, where α is the photon index. Both the x and y axes in (a) and (b) have arbitrary units. The comoving spectrum is assumed to be a SBPL. The several colors indicate the observed peak frequency at the beginning of the decay. The error bars represent 1σ uncertainties, calculated via spectral fitting in XSPEC. In the legend, we report the name of each GRB. NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-24246-x ARTICLE indicating that the spectral evolution becomes dominated by the emission at larger angles, rather than by adiabatic cooling in the jet core. Adiabatic cooling can also well reproduce the light curve of X-ray tails (Fig. 4b). For comparison, in the same plot, we show the light curve given by pure HLE, adopting the same value of R 0 and Γ.
In order to fully explore the parameter space of the adiabatic cooling model, we used a Monte Carlo Markov Chain (MCMC) algorithm for the parameter estimation. We consider the joint temporal evolution of flux and photon index and we find agreement of the model with data (see "Methods:" "Parameter estimation via MCMC" and Table 1). In Supplementary Figs. 9 and 10, we show for each burst the observed temporal evolution of photon index and normalized flux in comparison with the curves produced with 500 random draws from the posterior sample set of the MCMC. We obtain a value of λ in the range 0.4-0.7 (except for 090621A which prefers λ~2). On average, these values of λ are smaller than those expected in an emitting region with a transverse magnetic field (λ = 1 or λ = 2 for a thick or a thin shell, respectively) or magnetic field in pressure equilibrium with the emitting particles (λ = 4/3 or λ = 2 for a thick or a thin shell, respectively 31 ).
The typical timescale of adiabatic cooling τ ad = R 0 /2cΓ 2 , i.e., the observed time interval during which the radius doubles, is equal to the HLE timescale 13,32 and radically affects the slope of X-ray tails. Therefore, the comparison between the model and the observed light curves allows us to constrain the size R 0 of the emitting region as in HLE 33,34 . We find values in the range 0.3 s ≲ τ ad ≲ 24 s. These values are quite larger than the typical duration of GRB pulses (<1 s 35 ), which can be due to the following reason. For the spectral analysis to be feasible, we had to choose only tails that are long enough (to be divided down into a sufficient number of temporal bins). Moreover, the prompt emission is usually interpreted as a superposition of several emission episodes: the SD observed in XRT is likely dominated by the tails with the slowest decay timescales. Since, in our model, the decay timescale is τ ad = R/2cΓ 2 , this could indicate that the emission radius of the pulses that dominate the tail is systematically larger than that of pulses that dominate the prompt emission. If this is the case, a lower magnetic field is also expected, which goes well along with the long radiative timescale and slow (or marginally fast) cooling regime, in agreement with our results. For the range of τ ad obtained from the analysis, the corresponding range for the emission radius is 1.8 × 10 14 (Γ/100) 2 cm ≲ R 0 ≲ 1.4 × 10 16 (Γ/100) 2 cm. A different prescription for adiabatic cooling has been suggested in the literature 30 , in which the particle's momentum gets dynamically oriented transverse to the direction of the local magnetic field. In this case, HLE is the dominant contributor to the X-ray tail emission, which is again incompatible with the observed α − F relation.
Extending the sample. In order to further test the solidity of the α − F relation, we extend our analysis to a second sample of GRBs (composed by eight elements), which present directly a SD at the beginning of the XRT light curve, instead of an X-ray pulse (see Fig. 5a), often observed in early X-ray afterglows 8,9 . We require that the XRT SD is preceded by a pulse in the BAT light curve (the brightest since its trigger time). We add the data of this second sample to the α − F plot, estimating the peak flux by the extrapolation of the XRT light curve backwards to the peak time a b Fig. 4 Spectral and temporal evolution in case of adiabatic cooling. In a we show the α − F relation expected in the case of adiabatic cooling (solid lines). The theoretical curves are computed taking also into account the effect of HLE. The value of λ specifies the evolution of the magnetic field. We adopt a SBPL as spectral shape with α s = − 1/3 and β s = 1.5, an initial observed peak frequency of 100 keV and a thickness of the expanding shell that is constant in time. The dot-dashed line is the evolution expected in case of HLE without adiabatic cooling, assuming the same spectral shape and initial observed peak frequency. The error bars represent 1σ uncertainties, calculated via spectral fitting in XSPEC. In b we show the temporal evolution of normalized flux expected in case of adiabatic cooling. δt obs + 100 s is the time measured from the peak of the decay shifted at 100 s, the typical starting time of the tail emission detected by XRT. We adopt the same parameters as in a, assuming R 0 = 2 × 10 15 cm and Γ = 100. The dot-dashed line is the corresponding HLE model without accounting for adiabatic cooling. τ ad = R 0 /2cΓ 2 indicates the timescale of adiabatic cooling, which is the same of HLE. The vertical error bars represent 1σ uncertainties and they are calculated via spectral fitting in XSPEC, while horizontal error bars represent the width of the time bin. In the legend we report the name of each GRB. The confidence intervals and the lower limits represent the 16th, 50th, and 84th percentiles of the samples in the marginalized distributions (i.e., 1σ level of confidence).
of the BAT pulse, under the assumption that BAT and XRT peaks were simultaneous (see "Methods:" "Extrapolation of F max "). We find that these GRBs follow the overall α − F relation (Fig. 5b), confirming that a common physical process is governing the spectral evolution of X-ray tails. Adiabatic cooling is still capable of reproducing the data of this second sample (Fig. 6a), provided that we assume a slightly softer high-energy intrinsic spectrum (α~3 instead of α~2.5). Alternatively, the introduction of an exponential cutoff in the spectral shape at ν = ν c~νm can also reproduce the data (see Fig. 6b), where ν c and ν m are the synchrotron characteristic frequencies. The cutoff is formed by a combined action of adiabatic cooling and mild synchrotron cooling (see Supplementary Note 4). We specify that the limited size of our samples is related to the selection requirements, which are necessary for an appropriate time-resolved spectral analysis. Thus our results are proved for X-ray tails firmly connected to prompt emission pulses.

Discussion
The α − F relation, found in our analysis, requires a mechanism that produces the X-ray tails of GRBs with a unique law of flux decay and spectral softening. Although other scenarios cannot be ruled out, we find that adiabatic cooling of the emitting particles, together with a slowly decaying magnetic field, is the most plausible scenario able to robustly reproduce this relation. Our results suggest an efficient coupling between a slowly decaying magnetic field and the emitting particles. Our findings are generally in agreement with moderately fast and slow cooling regimes of the synchrotron radiation, which is able to reproduce the overall GRB spectral features 36 . In the adiabatic cooling scenario, most of the internal energy is not radiated away before the system substantially expands. If electrons are responsible for the emission, an extremely small magnetic field would be required [37][38][39] , which is unrealistic for this kind of outflows. Protons radiating Fig. 5 The steep decay phase and the correspondent spectral evolution for the extended sample. In a we show an example of a light curve of an Xray tail selected for our extended sample, taken from GRB 150323A. We report on the same plot the XRT (orange) and the BAT (blue) flux density at 1 and 50 keV, respectively. The peak flux F max is estimated extrapolating the X-ray tail back to the BAT peak. The error bars represent 1σ uncertainties and they are derived from the Swift archive. In b we show the spectral evolution of our extended sample of GRBs, which present a steep decay at the beginning of the XRT light curve, preceded by the brightest BAT pulse since the trigger time. The evolution of α lies on the same region of the plane occupied by the original sample, indicated in gray. The error bars represent 1σ uncertainties, calculated via spectral fitting in XSPEC. In the legend, we report the name of each GRB. a b Fig. 6 Spectral evolution expected in case of adiabatic cooling (solid lines) superimposed to the extended sample. In a, the theoretical curves are computed considering adiabatic cooling and inefficient synchrotron cooling, taking also into account the effect of HLE. The value of λ specifies the evolution of the magnetic field. We adopt a SBPL as spectral shape with α s = −1/3 and β s = 2.0, an initial observed peak frequency of 100 keV and a thickness of the expanding shell that is constant in time. In b we show the spectral evolution expected in case of combined adiabatic cooling and mild synchrotron cooling. The adopted spectral shape is a SBPL plus an exponential cutoff. The initial peak frequency is 100 keV. The theoretical curves are computed taking also into account the effect of HLE. In both panels, the error bars represent 1σ uncertainties, calculated via spectral fitting in XSPEC, and the dot-dashed line is the evolution expected considering only HLE, assuming the same spectral shape and initial observed peak frequency. In the legend, we report the name of each GRB. through synchrotron emission can solve this problem 40 . Due to their larger mass, they radiate less efficiently than electrons, explaining why adiabatic cooling dominates the spectral evolution.
In conclusion, our results indicate that adiabatic cooling can play a crucial role for the collective evolution of the radiating particles in GRB outflows and consequently for the determination of spectral and temporal properties of prompt emission episodes. The coupling between particles and magnetic field ensures the intrinsic nature and hence the universality of this process, whose effects are independent of the global properties of the system, such as the luminosity of the GRB or the geometry of the jet.

Methods
Sample selection. We define the SD segment 6-9 as the portion of the light curve that is well approximated by a power law, F ∝ t −α with α > 3. Such criterion allows us to exclude a decay coming from a forward shock [10][11][12] . In order to determine the presence of a SD, we analyze the light curve of the integrated flux in the XRT (0.3-10) keV band.
From the Swift catalog 41 as of the end of 2019, we selected all GRBs with an XRT peak flux F XRT p > 10 À8 erg cm À2 s À1 . We selected the brightest pulses in order to have a good enough spectral quality as to perform a time-resolved spectral analysis. The peak flux is computed taking the maximum of F(t i ), where F(t i ) are the points of the light curve at each time t i . Among these GRBs, we selected our first sample according to the following criteria: (1) The XRT light curve shows at least one SD segment that is clean, i.e., without secondary peaks or relevant fluctuations.
(2) If we call F 1 and F 2 , the fluxes at the beginning and at the end of the SD, respectively, we require that F 1 F 2 > 10. This requirement is necessary to have a sufficient number of temporal bins inside the SD segment and therefore a well sampled spectral evolution.
(3) The beginning of the SD phase corresponds to a peak in the XRT light curve, such that we have a reliable reference for the initial time. We stress that the identification of the SD starting time in XRT is limited by the observational window of the instrument. This means that, if the XRT light curve starts directly with a SD phase, with no evidence of a peak, the initial reference time is possibly located before and its value cannot be directly derived. (4) The XRT peak before the SD has a counterpart in BAT, whose peak is the brightest since the trigger time. This requirement is necessary to ensure that XRT is looking at a prompt emission episode, whose typical peak energy is above 100 keV. In a quantitative way, we define two times, t p and t stop 90 , where the first indicates the beginning of the peak that generates the SD, while the second is the end time of T 90 42 , with respect to the trigger time. We require t stop 90 > t p in order to have an overlap between the last prompt pulses (monitored by BAT) and the XRT peak that precedes the SD phase. Namely, such requirement ensures that a considerable fraction of the energy released by the burst goes into the pulse that generates the X-ray tail. It is possible that more than one peak is present in the XRT light curve, each with a following SD. In this case we consider only the SD after the brightest peak. If two peaks have a similar flux, we consider the SD with the larger value of F 1 F 2 . We then define a second sample of GRBs that satisfy the first two points listed before, but have a SD at the beginning of the XRT light curve, namely no initial peak preceding the SD is present. In addition, we require that a BAT pulse precedes the XRT SD and is the brightest since the trigger time. The BAT pulse enables us to constrain the starting time of the SD.
The selection criteria limit the size of our sample, but they are unavoidable to perform a well-targeted analysis of X-ray tails and to achieve robust conclusions about their origin.
Time-resolved spectral analysis. For each GRB, we divided the XRT light curve in several time bins, according to the following criteria: (1) Each bin contains only data in windowed timing (WT) mode or in photon counting (PC) mode, since mixed WT + PC data cannot be analyzed as a single spectrum. (2) Each bin contains a total number of counts N bin in the E = (0.3-10) keV band larger than a certain threshold N 0 , which is chosen case by case according to the brightness of the source (see below). The definition of the time bins is obtained by an iterative process, i.e., starting from the first point of the light curve we keep including subsequent points until where N(t n ) are the counts associated to each point of the light curve, while t i and t f define the starting and ending time of the bin. Then the process is repeated for the next bins, until t f is equal to the XRT ending time. Due to the large range of count rates covered during a typical XRT light curve, the choice of only one value for N 0 would create an assembly of short bins at the beginning and too long bins toward the end. Therefore, we use one value of N 0 for bins in WT mode (N WT 0 ) and a smaller value of N 0 for bins in PC mode (N PC 0 ). In our sample, the SD is usually observed in WT mode, therefore we adjust N WT 0 in order to have at least 4-5 bins inside the SD. A typical value of N WT 0 is around 1500-3000, while N PC 0 is around 500-1000. Using these values, we verified that the relative errors of photon index and normalization resulting from spectral analysis are below~30%.
(3) For each couple (N i , N j ) of points inside the bin, the following relation must hold: where σ i and σ j are the associated errors. Such requirement avoids large flux variations within the bin itself. (4) The duration of the bin is larger than 5 s, in order to avoid pileup in the automatically produced XRT spectra. It is possible that condition 3 is satisfied only for a duration of the bin T bin < T 0 , while condition 2 is satisfied for T bin > T Ã 0 , but T Ã 0 > T 0 , meaning that they cannot be satisfied at the same time. In this case, we give priority to condition 3, provided that N bin is not much smaller than N 0 .
Due to the iterative process that defines the duration of the bins, it is possible that the last points in WT and PC mode are grouped in a single bin with a too small N bin , giving a too noisy spectrum. Therefore, they are excluded from the spectral analysis.
Spectral modeling. The spectrum of each bin is obtained using the automatic online tool provided by Swift for spectral analysis (see "Data availability"). Each spectrum is analyzed using XSPEC 43 , version 12.10.1, and the Python interface PyXspec. We discard all photons with energy E < 0.5 keV and E > 10 keV. The spectra are modeled with an absorbed power law, and for the absorption, we adopted the Tuebingen-Boulder model 44 . If the GRB redshift is known, we use two distinct absorbers, one Galactic 45 and one relative to the host galaxy (the XSPEC syntax is tbabs*ztbabs*po). The column density N H of the second absorber is estimated through the spectral analysis, as explained below. On the other hand, if the GRB redshift is unknown, we model the absorption as a single component located at redshift z = 0 (the XSPEC syntax is tbabs*po) and also in this case the value of N H is derived from spectral analysis.
For the estimation of the host N H , we consider only the late part of the XRT light curve following the SD phase. At late time with respect to the trigger, we do not expect strong spectral evolution, as verified in several works in the literature 46,47 . Therefore, for each GRB, the spectrum of each bin after the SD is fitted adopting the same N H , which is left free during the fit. Normalization and photon index are also left free, but they have different values for each spectrum. We call N late H the value of N H obtained with this procedure. In principle, the burst can affect the ionization state of the surrounding medium, but we assume that such effects are negligible and N H does not change dramatically across the duration of the burst 48 . Hence, we analyzed separately all the spectra of the SD using a unique value of N H ¼ N late H , which is fixed during the fit. Normalization and photon index, instead, are left free.
An alternative method for the derivation of N H is the fitting of all the spectra simultaneously imposing a unique value of N H that is left free. On the other hand, since N H and photon index are correlated, an intrinsic spectral evolution can induce an incorrect estimation of N H . For the same reason we do not fit the spectra adopting a free N H , since we would obtain an evolution of photon index strongly affected by the degeneracy with N H .
In this regard, we tested how our results about spectral evolution depend on the choice of N H . On average, we found that the fits of the SD spectra remain good (stat/dof ≲1) for a variation of N H of about 50%. As a consequence, the photon index derived by the fit would change at most of 30%. Therefore the error bars reported in all the plots α − F are possibly underestimated, but even considering a systematic error that corresponds to~30% of the value itself would not undermine the solidity of the results.
Extrapolation of F max . We explain here how we extrapolated the F max for the GRBs of the second sample, for which the XRT light curve starts directly with a SD. We consider the peak time T BAT p of the BAT pulse that precedes the SD. In the assumption that the SD starts at T BAT p , we can derive F max using the following procedure. We consider the (0.5-10) keV flux F(t i ) for each bin time t i in the SD, derived from spectral analysis. Then we fit these points with a power law with s > 0 and imposing that t 0 ¼ T BAT p . Finally we derive the best fit value of F max with the associated 1σ error. The error of F max has a contribution coming from the error associated to s and another associated to t 0 , as well as from the assumption of a power law as fitting function. The value of T BAT p is obtained fitting the BAT pulse with a Gaussian profile. Since usually the BAT pulse can have multiple sub-peaks and taking also into account possible lags between XRT and BAT peaks, we adopt a conservative error associated to T BAT p equal to 5 s. HLE from infinitesimal duration pulse. We assume that an infinitesimal duration pulse of radiation is emitted on the surface of a spherical shell, at radius R 0 from the center of the burst. Such treatment implicitly assumes particles that cool on timescales much smaller than the dynamical timescales. Therefore, all the X-ray tail emission is dominated by photons departed simultaneously from the last emitting surface. The jet has an aperture angle ϑ j and it expands with a bulk Lorentz factor Γ. We assume also that the comoving spectrum is the same on the whole jet surface. The temporal evolution of the observed flux density is given by ref. 49 F ν ðt obs Þ / S ν 0 ðν=DðϑÞÞD 2 ðϑÞ cosðϑÞ ð 4Þ with S ν 0 ðν=DðϑÞÞ the comoving spectral shape, DðϑÞ the Doppler factor, and ϑ the angle measured from the line of sight, which is assumed to coincide with the jet symmetry axis. The observer time t obs is related to the angle ϑ through the following formula: where t em is the emission time. Equation (4) is valid for ϑ < ϑ j , while for ϑ > ϑ j the emission drops to zero. This implies that for t obs > t em ð1 À β cos ϑ j Þ, the flux drops to zero. At each time t obs (ϑ), the observer receives a spectrum that is Doppler shifted by a factor DðϑÞ with respect to the comoving spectrum. If the comoving spectrum is curved, i.e., if d 2 dν 0 2 S ν 0 ≠ 0, then also the photon index is a function of time 16 . The shape of the resulting curve α − F is determined only by the spectral shape and the comoving peak frequency ν 0 p , while it is independent on the emission radius R 0 and the bulk Lorentz factor Γ.
We notice that the observed photon index goes from 0.5-1.0 up to 2.0-2.5, consistent with the slopes of a synchrotron spectrum before and after the peak frequency. Indeed for a population of particles with an injected energy distribution N(γ) ∝ γ −p that has not completely cooled, the expected shape of the spectrum is F ν~ν 1/3 (α = 2/3) for ν < ν c and F ν~ν −p/2 (α = p/2 + 1) for ν > ν m ≳ ν c . Hereafter, if not otherwise specified, we assume a spectral shape given by a smoothly broken lower law, which well approximates the synchrotron spectrum below and above the peak frequency. The form of the adopted spectral shape is with α s = −1/3 and β s = 1.5. The peak frequency ν p of the energy spectrum νS ν is related to ν 0 through the following relation: At each arrival time, we compute the flux and the photon index in the XRT band using Eq. (4). In particular, the XRT flux is given by where h is Planck's constant, while the photon index is computed as 16,24 αðt obs Þ ¼ 1 À log½F ν¼10keV=h ðt obs Þ=F ν¼0:5keV=h ðt obs Þ logð10 keV=0:5 keVÞ ð9Þ This method for the evaluation of photon index is valid in the limit of a spectrum that can be always approximated with a power law as it passes through the XRT band, which is the case for typical prompt emission spectra. In addition to the SBPL, we test HLE also using other spectral shapes. We first adopt a Band function 23 with the following form: where ϵ = ν/ν 0 . In this case, the peak of the energy spectrum is at ν p = (2 + α s )ν 0 . The resulting spectral evolution is very similar to the case of SBPL, as visible in Supplementary Fig. 12a. As a final attempt, we use synchrotron spectrum emitted by a population of particles with an initial energy distribution N(γ) ∝ γ −p . Synchrotron is considered the dominant radiative process in prompt emission of GRBs 4,36 . In the fast cooling regime, the particle distribution becomes The only three parameters that define the shape of the synchrotron spectrum are ν m / γ 2 m , ν c / γ 2 c , and p. For the computation of the spectrum we use 50 with where B is the magnetic field and K 5/3 (x) is the modified Bessel function of order 5/3. The resulting spectral evolution for values of ν m /ν c = 1 and ν m /ν c = 10 is reported in Supplementary Fig. 11. A value of ν m /ν c~1 is expected in the marginally fast cooling regime [37][38][39] , which is favored by broad-band observations of GRB prompt spectra [51][52][53][54][55][56] . Finally, we test how the sharpness of the spectral peak can affect our results. In particular we consider again a SBPL and we generalize the formula adding a sharpness parameter n S ðnÞ where larger values of n correspond to sharper spectral peaks. As visible in Supplementary Fig. 12b, where we have adopted n = 4, the shape of the curves becomes flatter at the beginning and at the end of the decay, but with no substantial steepening of the intermediate part. This is attributable to HLE that imposes an evolution of the observed peak frequency like t À1 obs . Thus, while the initial and final values of photon index are dictated by the spectral shape, the steepness of the transition from the initial to the final value is governed by HLE and is independent on the spectral shape. In conclusion, no one of the alternative spectral shapes that we tested is able to reconcile HLE with the observed spectral evolution.
We finally test how the α − F relation changes if we assume a structured jet with an angle-dependent comoving spectrum. In particular, we consider a spectral peak energy that is nearly constant inside an angle ϑ c (measured with respect to the line of sight) and starts to decrease outside it. Regardless of the choice of the specific law for the angular dependence (e.g., Gaussian or power law), the HLE can reproduce the α − F relation only if all the analyzed GRBs have a fine-tuned value of ϑ c < 1 ∘ . Such a small value of ϑ c , on the other hand, would imply a very short SD, in contradiction with observations. Adiabatic cooling. In this section, we derive the effect of adiabatic cooling of the emitting particles 57 on the light curve and the spectral evolution of X-ray tails. We assume that the emission is dominated by a single species of particles that can be treated as a relativistic gas in adiabatic expansion. We assume also that there is no interaction with other species of particles. If the particles are embedded in a region of comoving volume V 0 , an adiabatic expansion satisfies the equation where 〈γ〉 is the average Lorentz factor of the emitting particles in the comoving frame. The last equation is valid in the limit in which the adiabatic cooling timescale is smaller than the cooling time of other radiative processes, such as synchrotron or inverse Compton. Namely, particles radiate only a negligible fraction of their internal energy during the expansion of the system. Regarding the radial dependence of the volume V 0 , we distinguish two cases: (1) thick shell, with a comoving width ΔR 0 that does not evolve with time, hence V 0 / R 2 ΔR 0 / R 2 (2) thin shell, with a comoving width ΔR 0 that evolves linearly with R, hence We assume that the dominant radiative process is synchrotron. The evolution of the spectrum in the observer frame is therefore fully determined once we know how the spectrum normalization F ν p and the peak frequency ν p evolve in time. These two quantities, under the assumption of constant total number of emitting particles and constant bulk Lorentz factor Γ, take the following form: where B is the magnetic field (assumed tangled) as measured by a comoving observer. As described in the main text, we adopt the following parametrization for the magnetic field where λ ≥ 0, under the reasonable assumption that magnetic field has to decrease or at most remain constant during the expansion.
Parameter estimation via MCMC. In order to fully explore the parameter space of the adiabatic cooling model, we performed a MCMC, using the emcee algorithm 58 . The setup of our analysis is described in the following points: (1) The model contains as free parameters E p , λ and τ ad = R/2cΓ 2 , which are the peak energy at the beginning of the SD, the decay index of the magnetic field, and the adiabatic timescale. The inclusion of Γ as free parameter returns a flat posterior distribution, indicating that the model is insensitive to it. Therefore we performed the analysis fixing Γ = 100.
where ϕ ¼ F=F max , α is the photon index, and with ϕ, α we indicate the value predicted by the model at each time t n . Moreover, s 2 ϕ;n ¼ σ 2 ϕ;n þ f ϕ Á ϕðt n Þ; s 2 α;n ¼ σ 2 α;n þ f α Á αðt n Þ ð20Þ where σ ϕ and σ α are the errors, while f ϕ and f α are introduced to take into account possible underestimation of the errors. Since keeping f α ≠ 0 leads to a posterior distribution of f α peaked around~10 −7 , the parameter estimation was performed fixing f α = 0. Instead, we keep f ϕ ≠ 0 taking into account that the error on the flux resulting from the fitting of time-averaged spectrum may not represent the true flux error over the time bin. (3) The MCMC runs until the number of steps exceeds 100 times the autocorrelation time (its maximum) and the averaged autocorrelation time (over 100 steps) becomes constant within 1% accuracy. The burn-in is chosen as twice of autocorrelation time. As an example, we show in Supplementary Fig. 7, the evolution of the autocorrelation time as a function of the steps. The resulting parameter estimation is summarized in Table 1. The uncertainties are reported based on the 16th, 50th, and 84th percentiles of the samples in the marginalized distributions (1σ level of confidence). An example of corner plot obtained via MCMC is shown in Supplementary Fig. 8. In Supplementary Figs. 9 and 10, we show for each burst the observed temporal evolution of photon index and normalized flux in comparison with the curves produced with 500 random draws from the posterior sample set of the MCMC.
We performed an analogous MCMC analysis adopting the model of HLE from an instantaneous emission. However, the algorithm is unable to converge, demonstrating that the model cannot successfully match with the observations. The only way to obtain converged chains by this model is to admit extreme and unrealistic values of f ϕ , of the order 10 4 -10 8 . The only exception is GRB 090621, which can be fitted by HLE alone. This is the only case where it is meaningful to compute the Bayes factor between HLE and AC, which results to be~200. Thus we prove that the adiabatic cooling model is strongly preferred for all the analyzed cases.
The model comparison (adiabatic cooling + HLE against HLE only) is also done assuming different spectral shapes: SBPL, Band, and synchrotron. The spectral parameters are the same as those adopted before. For synchrotron, we use ν c = ν m (the case ν c ≠ ν m does not improve the goodness of fit). In order to compare the goodness of fit of the two models, we used the Akaike information criterion (AIC), which is defined as AIC ¼ 2k À 2ln ðLÞ, where k in the number of parameters of the model, and L is the best fit likelihood, that is, 2ln ðLÞ ¼ Àχ 2 . In Table 2, we show the value of Δ AIC = AIC HLE − AIC AC for each spectral shape. For all cases, the adiabatic cooling is significantly favored with respect to HLE.

Data availability
Raw data are public and available in the UK Swift Science Data Centre at the University of Leicester. The light curve data are taken at this link: https://www.swift.ac.uk/ xrt_curves/GRB_ID/flux.qdp where GRB_ID is the GRB observation ID. The spectra are obtained at this link: https://www.swift.ac.uk/xrt_spectra/addspec.php?targ=GRB_ID where GRB_ID is the ID number of the GRB. The details of the automatic spectral analysis can be found here: https://www.swift.ac.uk/xrt_spectra/docs.php. Derived data are available from the corresponding author on request.