© The Authors 2025

1. Introduction

Cosmic explosions triggered by collapses of massive stars (Woosley 1993; Woosley & Bloom 2006) or binary neutron star mergers (Eichler et al. 1989; Narayan et al. 1992; Berger 2014; Abbott et al. 2017) produce gamma-ray bursts (GRBs), the most luminous transient phenomena in the Universe (Paczynski 1986). Thanks to several decades of observations, GRBs have been revealed to be associated with highly collimated ultra-relativistic jets launched in the aftermath of their progenitor’s explosion (see Salafia & Ghirlanda 2022, for a review).

The very fast prompt emission variability (e.g., Bhat et al. 2012) clearly suggests an emission site internal to the jet (Sari & Piran 1997). A natural explanation for the jet kinetic energy dissipation is through internal shocks (Rees & Meszaros 1994). However, sub-photospheric dissipation (Rees & Mészáros 2005; Pe’er 2008) or magnetic reconnection processes (Drenkhahn & Spruit 2002; Zhang & Yan 2011) are equally able to explain the GRB radiative output. The early emission observed soon after the explosion, during the so-called prompt emission phase, shows a broad and peaked spectral energy density (SED) that is possibly associated with nonthermal radiative processes. Historically, prompt emission spectra were modeled by the phenomenological Band function (Band et al. 1993); namely, two power laws smoothly connected around a peak at energy E_p. The systematic fit of this model to large samples of GRB spectra revealed that in some cases the low-energy photon index, α, was consistent with the one predicted by synchrotron emission in a fast-cooling regime, whereas some other GRBs exhibited very hard α, pointing toward quasi-thermal processes (Ghisellini & Celotti 1999; Lazzati et al. 2000; Mészáros & Rees 2000). The vast majority of the GRBs showed average behavior, with spectra that were too hard to be produced by synchrotron processes, but too soft to indicate thermal origins (Preece et al. 1998; Kaneko et al. 2006; Gruber et al. 2014).

This tension was partially mitigated when more complex models were employed and low-energy spectral breaks were identified in some GRB spectra at energies, E_b. It was shown that a single synchrotron spectrum in a marginally fast-cooling regime (Oganesyan et al. 2017, 2018; Ravasio et al. 2019), a slow-cooling regime (Zhang et al. 2009; Burgess et al. 2020), or even synchrotron self-absorption (Lloyd & Petrosian 2000; Lloyd-Ronning & Petrosian 2002) are able to account for the entire prompt emission spectrum, from the optical (Oganesyan et al. 2019) up to the GeV band (Macera et al., in prep.). Despite the synchrotron model seeming to be a viable explanation for most of the considered spectra, GRBs show diverse spectral and temporal behaviors, preventing a convincing unified interpretation.

A feature that appears to be shared among GRBs is that their prompt spectrum is harder when the GRB is brighter and more energetic. In fact, many sources show that the rest-frame spectral peak energy, E_p, z, is positively and tightly correlated with the isotropic-equivalent energy, E_iso (Amati et al. 2002), and luminosity, L_iso (Yonetoku et al. 2004). These two relations span several decades in both energies and luminosities, and despite being affected by selection and observational biases (Band & Preece 2005; Nakar & Piran 2005; Butler et al. 2007; Shahmoradi & Nemiroff 2011), they seem to hold for most of the GRB class (Ghirlanda et al. 2008; Nava et al. 2008), and also at different redshifts (Nava et al. 2012).

In particular, the E_p, z − E_iso relation (better known as the “Amati relation”) was initially found for long GRBs, with durations of T₉₀ > 2 s. Short GRBs (T₉₀ < 2 s) exhibit a similar trend but occupy a parallel track in the E_p, z − E_iso plane (Ghirlanda et al. 2009), showing a systematically lower E_iso for the same E_p, z. On the other hand, both long and short GRBs follow the same E_p, z − L_iso relation (better known as the “Yonetoku relation”, Yonetoku et al. 2010). After their first discoveries, many outliers to these relations were found; however, both the Amati and Yonetoku relations still hold when multiple time bins from the same bright GRB are considered through time-resolved spectral fits (Ghirlanda et al. 2010, 2011a,b).

At odds with the Amati relation, which is obtained comparing quantities averaged across the whole burst duration, the Yonetoku relation is found considering L_iso (but also E_p,z, e.g., Tsvetkova et al. 2017) estimates only during the brightest pulse; in other words, at the peak of the light curve. This makes the Yonetoku relation inherently connected with the physics taking place in single GRB pulses, and it appears particularly well suited for inquiring into prompt emission physics in depth.

For more than two decades, the observations carried out by gamma-ray instruments such as the Burst Alert Telescope (BAT: 15–150 keV) on board the Neil Gehrels Swift Observatory (Swift), the Gamma-ray Burst Monitor (GBM; 8 keV–40 MeV) on board the Fermi satellite, and Konus-Wind (KW, ∼20 keV–20 MeV) have further corroborated the validity of these correlations (e.g., Nava et al. 2012; Tsvetkova et al. 2017; Minaev & Pozanenko 2020; Wang et al. 2024). Despite all these findings pointing toward a possible common emission mechanism occurring in GRB central engines, a clear physical interpretation of these correlations is still missing (see Parsotan & Ito 2022, for a review).

In this work, we carry out spectral fits of GRB samples detected by Fermi/GBM using not only the phenomenological Band function, but also a physical synchrotron model. Our goals are (i) to test whether using a physical model allows one to derive the E_p,z−L_iso (Yonetoku) relation and (ii) to explore any possible connection between the Yonetoku relation and physical synchrotron parameters such as the characteristic minimum frequency, ν_m, and the cooling frequency, ν_c.

The paper is organized as follows. In Sect. 2, we describe the sample selection and the spectral modeling; in Sect. 3, we show the results of the spectral fit to the whole sample of GRBs (single-bin sample) and the time-resolved analysis of a subsample of bright GRBs (multiple-bin sample); in Sect. 4, we discuss the results; and finally, in Sect. 5, we summarize the conclusions. Throughout the paper, we assume standard cosmology with h = Ω_Λ = 0.7 and Ω_m = 0.3. If not stated differently, errors are reported at a 68% confidence level. As a convention, the subscript z refers to rest-frame parameters.

2. Methods

2.1. Sample selection

We analyzed GRBs observed by the Fermi satellite before September 2023. We selected the GRBs with a firm redshift estimate, provided by the MPE GRB online catalog¹. We considered only the GRBs that are bright enough during their main pulse to perform a robust spectral analysis in that temporal bin. Specifically, we included in our sample only the GRBs whose 1s peak photon flux in the energy range 50–300 keV is P  ≥  3.5 photons s⁻¹ cm⁻², which corresponds to nearly seven times the flux sensitivity of Fermi/GBM. We excluded GRB221009A, GRB230307A, and GRB130427A, the three brightest GRBs, since their emission during the peak are affected by pileup effects and their GBM data cannot be used for spectral analysis.

This leads to a total sample of 74 GRBs, with a spectroscopic redshift estimate for 70 of them. Three GRBs (GRB100816A, GRB200826A, and GRB211211A) have a redshift estimate through host galaxy association, while one GRB (GRB200829A) has a photometric redshift estimate.

2.2. Data reduction

We downloaded the Fermi/GBM data (8 keV–40 MeV) of all the 74 GRBs in our sample from the Fermi GBM Burst catalog², and performed a standard data reduction using the Fermi science tool GTBURST³. In particular, we considered data from the two sodium iodide (NaI, 8–900 keV) and one bismuth germanate (BGO, 0.3–40 MeV) detectors with best observational conditions (i.e., lowest viewing angles). In three GRBs (GRB121128A, GRB131231A, and GRB140508A), the second best NaI was exhibiting artefacts in the low-energy data (< 30 keV), leading to spectral deviations with respect to the best NaI. For this reason, for these three GRBs we considered the best and third-best NaI, together with the best BGO detector.

The background analysis was performed through GTBURST by using, whenever the polynomial fit was converging, the background intervals reported by Fermi collaboration (von Kienlin et al. 2020). In case this selection was leading to a non-convergent background polynomial fit, larger custom background intervals were selected.

We selected the source interval to coincide with the peak of the GRB light curve. The time bin is the one reported in the bcat file, a file containing basic burst parameters provided by Fermi collaboration for every GRB. The width of the time bin is 1.064 s for long GRBs (T₉₀ > 2 s) and 256 ms for short GRBs (T₉₀ < 2 s), in order to pinpoint their particularly short peak phase. The outputs of the data reduction are the source, background, and weighted response spectral files that are used for the spectral analysis through the Heasarc package XSPEC⁴ (Arnaud 1996).

The addition of high-energy (> 30 MeV) data to the prompt emission proved to be crucial in determining spectral parameters that would otherwise be inaccessible or difficult to constrain. For this reason, we decided to include in our dataset Large Area Telescope (LAT) Low Energy (LLE, 30–100 MeV) data whenever available. We retrieved the data from the Fermi LAT Low-Energy Events Catalog⁵. The LLE data were reduced using the same tools employed in the Fermi/GBM analysis, selecting the source in the same time bin in order to perform a simultaneous joint time-resolved analysis during the brightest peak of the emission. In our current sample, only 22 GRB have available LLE data during their peak. For those GRBs, we produced spectral files similarly to the GBM ones.

2.3. Spectral fit routine

After the data reduction, we modeled the spectra in our sample. For the fitting process, we used PYXSPEC⁶, the Python interface to the XSPEC spectral-fitting program. We implemented the Python-based software Bayesian X-ray analysis (BXA, Buchner 2016), which allows for Bayesian parameter estimation and model comparison through nested-sampling algorithms in PYXSPEC.

We ignored the energy channels outside 8–900 keV for the NaI detectors, as well as the 30–40 keV band in order to avoid the iodine K-edge line at 33.17 keV (Meegan et al. 2009). We selected the energy range 300 keV–40 MeV and 30–100 MeV for BGO and LLE data, respectively. For each dataset, consisting of GBM and LLE data (when available), we tested two models: the phenomenological Band function (grbm in XSPEC notation, Band et al. 1993) and an idealized physical synchrotron model. We included the presence of a cross-calibration constant, constant in XSPEC notation, allowing for a 30% variation for each dataset. We applied Poisson-Gaussian statistics (pgstat) to GBM data and Cash statistics (cstat) to LLE data. For each model, we measured the logarithm of the 10 keV–10 MeV bolometric flux in the source rest frame log F using the convolutional model cflux in XSPEC.

Moreover, in two GRBs (GRB090902B and GRB190114C) the prompt spectrum exhibits in NaI data a power law component (Abdo et al. 2009; Ajello et al. 2020; Ursi et al. 2020) in addition to the main emission one, also at low energies. Therefore, we added a power law in both the Band and synchrotron models, powerlaw in XSPEC notation. In these two cases, we computed the flux of the main component only (Band or synchrotron), since the physical nature of this additional power law was already discussed in the literature and it is outside the scope of this work.

For the synchrotron case, we used the same table model presented in Oganesyan et al. (2019). Synchrotron spectra were produced by a population of nonthermal electrons injected as a power law, dN_e/dγ ∝ γ^−p, with a characteristic minimum Lorentz factor, γ_m. The final spectrum was obtained by convolving the single electron spectrum with the electron population distribution after the particles cooled through synchrotron losses down to the cooling Lorentz factor, γ_c. Table spectra were produced by varying log(γ_m/γ_c) between –1 and 2, allowing for both fast (γ_m > γ_c) and slow (γ_m < γ_c) cooling regimes, and by varying p between 2 and 5 and the normalization, N_Sync, between 0 and 10²⁰ (in cgs units). However, this table model assumes a constant cooling frequency, ν_c = 1 keV. To let this parameter vary, we multiplied the table model with the XSPEX convolutional model zmshift, which shifts the spectrum along the energy axes. The resulting photon spectrum depends on four parameters: the ratio between the Lorentz factors, γ_m/γ_c, the power law slope, p, of the injected electron distribution, the spectral shift parameter, z_shift (univocally linked with the observed cooling frequency, ν_c), and the normalization, N_Sync.

2.4. Parameter estimation, model comparison, and goodness-of-fit

For each free parameter of the two models, we defined broad and uninformative priors (Table 1). The analysis returns posteriors of the parameters fit to our dataset together with the Bayesian evidence.

Being interested in physical parameters too, we derived them from the inferred ones. In the case of the Band model, we evaluated the rest frame peak energy defined as $E_{p,z}^{\mathrm{Band}}=(2+\alpha )·E_{\mathit{ch}}·(1+z)$, where α is the low-energy photon index and E_ch is the Band characteristic energy.

In the case of the synchrotron model, we defined the rest-frame synchrotron frequency associated with γ_c as ν_c, z = (1 + z)/(1 + z_shift) and with γ_m as ν_m, z = ν_c, z ⋅ (γ_m/γ_c)² and the rest frame peak energy as $E_{p,z}^{\mathrm{Sync}}=max({\nu }_{m,z},{\nu }_{c,z})$. Here, z_shift is the spectral shift parameter.

In both cases, we inferred the logarithm of the 10 keV – 10 MeV rest frame flux, log F, from which we computed the flux, F = 10^{log F}, and the total isotropic-equivalent luminosity, L_iso = 4πD_L²(z)⋅F, where D_L(z) is the luminosity distance from the burst. For each relevant parameter, we defined the best-fit value as the median on the posterior distribution, and lower and upper errors were derived from the 16th and 84th percentile of the posterior distribution, respectively. To ensure the good convergence of the fit, we modified the prior of the cooling frequency for GRB170607A between ν_c ∈ (8, 500) keV; that is, z_shift ∈ ( − 0.998, −0.8).

From the Bayesian fit of both models to each GRB dataset, we get the Bayesian evidence, log Z. This provides a tool with which to compare which model better describes the data. In fact, the model that has substantially larger evidence is the best-fit model. To assess this, we defined the Bayes factor for these two models as log B = log(Z_Sync/Z_Band), and we set the threshold to 0.5 according to Jeffrey’s scale (Kass & Raftery 1995). If |log B|> 0.5, we can identify a best-fit model, which is the synchrotron one if log B > 0 or the Band function if log B < 0. If |log B|< 0.5, the Bayesian evidence of one model is not enough to statistically exclude the other.

In addition, we aim to establish whether or not the synchrotron model provides a good fit to the data. To do so, we performed posterior predictive checks (PPCs) for each synchrotron fit to our sample. By fitting a given model to a dataset, the best-fit model is the realization that maximizes the relative likelihood, and therefore returns the smallest statistics. To assess whether that model realization is a good fit to the data, we simulated 1000 fake spectra out of the best-fit model, using the XSPEC command fakeit. We estimated, for each fake spectrum, the likelihood of the best-fit model out of the fake spectrum dataset.

By repeating this computation for all the 1000 fake spectra, we can define a (Gaussian-like) probability density function relative to the statistics distribution. We measured the goodness-of-fit through a comparison of the simulated statistics distribution and the measured likelihood value. Specifically, we computed the p value of the real statistics on the simulated distribution, and we considered it a good fit if p > 0.05. We defined the p value as the integral of the simulated statistics probability density function between the best-fit statistics value and +∞ or −∞, depending on whether the measured statistics is greater or lower than the maximum of the distribution, respectively.

2.5. Linear fit routine

To investigate the Yonetoku relation, we fit the following power law relation to the data of our samples:

$$\begin{array}{c}\hfill \frac{E_{p,z}}{100\mathrm{keV}}=K·{\left(\frac{L_{\mathrm{iso}}}{{10}^{52}\mathrm{erg}/\mathrm{s}}\right)}^m,\end{array}$$(1)

where m and K are the power law slope and the normalization, respectively. More practically, we used the logarithmic form of Eq. (1) in order to perform a linear fit:

$$\begin{array}{c}\hfill log{E}_{p,z}=m·(log{L}_{\mathrm{iso}}-52)+logK-2.\end{array}$$(2)

We normalized both E_p,z and L_iso in order to perform the linear fit closer to the data barycenter, providing better m and K estimates. Values and 1σ errors for log E_p, z and log L_iso were derived directly from E_p,z and L_iso parameter distributions, respectively.

We performed the fit using a Bayesian approach. We defined a likelihood that assumes Gaussian parameter distributions for both log E_p,z and log L_iso and takes into account errors on both parameters with the addition of an intrinsic Gaussian noise term, σ_sc (D’Agostini 2005):

$$\begin{array}{cc}& -2ln\mathcal{L}(m,K,{\sigma }_{\mathit{sc}}|\{x_i,{\sigma }_{x_i},y_i,{\sigma }_{y_i}\})\hfill \\ \hfill & ={\displaystyle \sum _{i=0}^N}[ln({\sigma }_{\mathit{sc}}^2+{\sigma }_{y_i}^2+m^2{\sigma }_{x_i}^2)+\frac{{(y_i-m(x_i-52)-logK+2)}^2}{({\sigma }_{\mathit{sc}}^2+{\sigma }_{y_i}^2+m^2{\sigma }_{x_i}^2)}],\hfill \end{array}$$(3)

where N is the number of GRBs in a given sample, (x_i, σ_xi) are log L_iso mean and standard deviation relative to the i-th GRB, and (y_i, σ_yi) are log E_p,z mean and standard deviation relative to the i-th GRB, respectively. With this approach, together with m and K, we estimated the intrinsic scatter, σ_sc, which constitutes a new free parameter.

Given the shape of the log L_iso parameter distribution derived from the spectral fit, we note that the hypothesis of a Gaussian distribution constitutes a good approximation. For this reason, we defined σ_xi as the quadrature sum of its upper and lower errors. On the other hand, the shapes of log E_p,z parameter distributions are asymmetric and deviate from a classic normal distribution. To account for this inconsistency, we still approximated log E_p,z as a Gaussian, but we considered σ_yi to be equal to the log E_p,z upper or lower error if the data point of the i-th GRB was above or below the best-fit line, respectively.

We adopted uniform priors for the three parameter, spanning the ranges m ∈ (0, 5), K ∈ (0.01, 100), and σ_sc ∈ (0, 100). We sampled the posterior probability density with a Markov chain Monte Carlo approach using the emcee Python package (Foreman-Mackey et al. 2013), employing N_walk = 32 walkers. We initialized the walkers in a small three-dimensional ball around a point in our parameter space, (m, K, σ_sc) = (0.5, 1.8, 10), and we performed N_iter = 5000 iterations, for a total of N_iter × N_walk = 160 000 samples.

3. Results

The main results of the spectral analysis fit to our single-bin sample are shown in Tables A.1 and A.2. From the Bayes factor comparison, we see that 21 GRBs are best fit by the synchrotron model and 40 GRBs are best fit by the Band model, while the two GRBs with a second component during the prompt emission peak prefer a Band model together with a power law. For 11 GRBs, the Bayes factor modulus, |log B|, is not large enough to prefer one model over the other.

However, the comparison between Band and the synchrotron models is not straightforward. Band is a phenomenological model, and it is able to accommodate a variety of spectra by providing measures of photon indices and peak energy. On the other hand, the synchrotron model hinders many physical assumptions, providing precious insights into the nature of the emission, but still having less freedom in describing an observed spectrum.

Despite the Band function being the best-fit model for most of the GRBs in our sample, we notice that the synchrotron model provides a good fit to the data in the majority of the cases. For this reason, we split the sample in two groups: the synchrotron sample (with 58 GRBs out of 74, Table A.1) and the complementary Band sample (with 16 GRBs out of 74, Table A.2). The latter includes all the GRBs that meet at least one of the following conditions: (1.) the synchrotron model provides a good fit of the data according to the PPC criterium; (2.) the best-fit model is the synchrotron one; and (3.) the Bayes factor is not able to prefer one model to the other one. For these GRBs, we report parameter estimates obtained from the synchrotron fit.

Interestingly, all of the GRBs in this sample are in a fast-cooling regime (ν_m/ν_c >  1); therefore, in Table A.1 we report L_iso, the peak energy, E_p,z = ν_m, z and the cooling frequency, ν_c, z, both in the rest-frame. On the other hand, the former sample includes all the GRBs where the synchrotron model does not provide a good fit to the data, and therefore we report parameter estimates obtained from the Band model fit.

In the synchrotron subsample, all the spectra are well fit by the synchrotron model according to the PPC criterium, except for three GRBs: GRB131231A, GRB150517A, and GRB211211A. For the first two GRBs, p − value < 0.5 because the best-fit statistics is particularly smaller than the average simulated statistics; that is, it has a significantly larger likelihood. In the case of GRB211211A, the Bayes factor prefers the synchrotron model despite the very large statistics compared to the simulated ones. In our synchrotron analysis, we get a comparable statistics with respect to the one reported by Gompertz et al. (2023) fitting a double smoothly broken power law function to the peak spectrum of this burst. Because of the model comparison, the results driven by analyzing the single-bin synchrotron subsample are valid under the assumption of synchrotron emission producing the prompt emission spectra.

3.1. E_p,z−L_iso relation

We first tested the E_p,z−L_iso relation obtained from the analysis of the whole sample. We computed the Spearman’s rank correlation coefficient, ρ, and its associated chance probability, P_chance. We report the results in Table 2, together with the ones of the linear fit. In Fig. 1 (left panel), we show that the data are scattered around the best-fit line, suggesting the absence of any correlation between L_iso and E_p,z (ρ ∼ 0.5). In addition, the slope we obtain deviates from the one expected in the Yonetoku relation.

Fig. 1. Yonetoku relation obtained from the single-bin sample. Black stars represent GRBs from the Band subsample, while colored dots represent GRBs from the synchrotron subsample. The color scale is associated with the frequency ratio, νm/νc, derived from the synchrotron fit. We fit a power law to the whole single-bin sample (left-hand panel) and to the GRBs in the synchrotron subsample in an intermediate cooling regime (right-hand panel). The straight black line represents the best-fit line from the linear fit, and the green-shaded area the relative 3σ scatter region. For comparison, the dashed black line represents the best-fit line from Yonetoku et al. (2004).

Nonetheless, we observe that GRBs with different frequency ratios, ν_m/ν_c, populate the E_p,z−L_iso plane differently. In fact, GRBs with a high ratio (i.e., in a fast-cooling regime, ν_m/ν_c ≳ 10) tend to deviate from the Yonetoku relation, whereas GRBs with a lower ratio (i.e., in an intermediate-cooling regime, ν_m/ν_c ≲ 3) tend to cluster along the Yonetoku relation.

By restricting the synchrotron sample to the GRBs in the intermediate-cooling regime (1< ν_m/ν_c< 3) and fitting the relation, we find that the data points are tightly correlated, and the linear fit returns values similar to the ones obtained in Yonetoku et al. (2004) (Fig. 1, right panel). The results of this linear fit are also shown in Table 2.

3.2. A physical relation: ν_c, z − L_iso

The E_p,z−L_iso analysis previously discussed shows that, when the physical synchrotron model is used, the Yonetoku relation does not hold for the whole sample, but for only a subset of GRBs that show a synchrotron spectrum where 1< ν_m/ν_c< 3 (Fig. 1, left panel). This constitutes an important connection between the Yonetoku relation and a physical parameter; namely, the ratio between the rest-frame cooling frequency, ν_c, z, and the minimum frequency, ν_m, z. In particular, the GRBs exhibiting a relation between their E_p,z and L_iso are the ones whose spectral break energy is particularly close to their peak energy.

Therefore, we tested the possibility of having a correlation between the break energy, namely the cooling frequency, ν_c, z, and the isotropic-equivalent luminosity, L_iso. We employed the same methodology used to investigate the E_p,z−L_iso relation, but applied to the synchrotron subsample of the single-bin sample (Table A.1).

We find that all of the GRBs in the synchrotron subsample, regardless of their ν_m/ν_c, show a tight ν_c, z − L_iso correlation (Fig. 2). A Spearman’s rank correlation test returns a coefficient of ρ = 0.81 with P_chance = 1.91 × 10⁻¹⁴. The slope of the power law is consistent with both the one expected from the Yonetoku relation and the one from the E_p,z−L_iso relation of GRBs in the synchrotron subsample. The results are reported in Table 2.

Fig. 2. νc, z − Liso relation obtained from the synchrotron single-bin sample. The color scale is associated with the frequency ratio, νm/νc, derived from the synchrotron fit. The straight black line represents the best-fit line from the linear fit, and the green-shaded area the relative 3σ scatter region.

3.3. Multiple-bins analysis

Of the 74 GRBs in our total sample, we selected the 15 GRBs with the highest fluence (reported by Fermi/GBM catalog) in order to perform time-resolved spectral analyses across their whole duration. Among these 15 GRBs, we found GRB090902B and GRB190114C; namely, the two GRBs exhibiting a second power law spectral component in addition to the main one. We removed these two GRBs from the sample, and the remaining 13 GRBs constitute the new multiple-bin sample. The time-resolved spectral analysis relative to this sample includes both Band and synchrotron models, similar to the procedure described in Sect. 2. We visually selected the time bins from the background-subtracted light curve of the best NaI detector in GBM. The time-bin start and stop were chosen in order to match as closely as possible pulses or single emission episodes shown in the GRB light curve. The time-bin width was chosen in order to collect enough signal for the spectral analysis. In case of particularly bright pulses, they were divided into multiple time bins given the high signal. In this sample, every time bin has a signal-to-noise ratio of SNR > 30. We included LLE spectral data in any time bin where they were available.

The multiple-bin sample covers a total of 125 spectra. The GRB list, their time bins, and the spectral fit results are reported in Table B.1. According to the Bayes factor, 65 spectra are best fit by synchrotron model, 48 spectra are best fit by Band, and for 12 spectra the Bayes factor cannot determine which is the best-fit model. For the 48 spectra best fit by Band, we report the E_p,z, L_iso, and α estimates obtained from the Band model, whereas for the remaining 77 spectra we report E_p,z, L_iso, and ν_c, z obtained from the synchrotron model.

Similarly to the analysis of the single-bin sample, we observe that the fit with the synchrotron model does not return any E_p,z−L_iso relation. However, restricting the multi-bin sample to its synchrotron subsample (i.e., the 77 GRBs for which we report synchrotron results), we observe again a very tight correlation between ν_c, z and L_iso (Fig. 3). The Spearman’s rank correlation coefficient is ρ = 0.83 with chance probability of P_chance = 1.83 × 10⁻²⁰ (Table 2).

Fig. 3. νc, z − Liso relation obtained from the synchrotron GRBs in the multi-bin sample. We show multiple time bins from different GRBs with circles of the same colors. The number inside the circles represents the bin number that data point is associated with. The straight black line represents the best-fit line from the linear fit, while the dashed-line the best-fit line obtained from the νc, z − Liso relation in the single-bin sample. The green-shaded area shows the 3σ scatter region of the relation.

We notice that the ν_c, z − L_iso relation derived in this way is described by a slightly flatter best-fit line, with a slope of m = 0.41 ± 0.03 and normalization of K = 1.12 ± 0.09, whereas the intrinsic scatter is comparable (σ_sc = 0.31). It is worth mentioning that, if the GRB211211A time-resolved dataset is considered alone, it is consistent with a slope of m ∼ 0.5. If excluded from the multiple-bin sample, the relative ν_c, z − L_iso relation has the same slope as the one obtained from the single-bin one.

4. Discussion

We found a novel ν_c, z − L_iso relation by employing a physical synchrotron model to describe the prompt emission of GRBs observed by Fermi/GBM. The fact that this new relation is found with high significance in the single-bin and multiple-bin samples, where both the emissions during the peak and across the whole burst duration are explored, strongly points toward its possible physical nature. In the following, we provide an interpretation of these results.

In the internal shocks framework, synchrotron spectra have a peak energy, E_p,z, which depends on many other parameters, most relevantly the bulk Lorentz factor, Γ, and the dynamical timescale, t_var (Rees & Mészáros 2005). The presence of a E_p,z−L_iso correlation would imply that Γ and t_var are the same for all the GRBs. This condition is not supported by any observation in large GRB samples. Therefore, this scenario is not in accordance with the Yonetoku relation, and our findings point toward this conclusion.

As a matter of fact, when time-resolved spectral analysis is performed, E_p,z exhibits a strong spectral evolution, which can track the intensity of the pulse (e.g., Golenetskii et al. 1983) or show a hard-to-soft transition (e.g., Lu et al. 2012). The latter case is not consistent with the Yonetoku relation when multiple pulses are present. On the other hand, Ravasio et al. (2019) performed time resolved analysis of a sample of ten bright long GRBs, fitting their spectra with a double smoothly broken power law (2SBPL) function, namely a Band function with two power laws before the peak. They find that the break energy, consistent with the synchrotron break in a marginally fast-cooling regime, does not vary much across the GRB duration. This finding suggests that the break energy may be more “standard” than the peak energy, or at least have minor dependencies on other parameters. However, the values of peak and break energies obtained through 2SBPL fits do not necessarily coincide with synchrotron breaks due to the shape of synchrotron spectra.

The question arises of whether and how the ν_c, z − L_iso relation can agree with almost two decades of observations pointing toward correlations with the peak energy, E_p,z. A synchrotron SED in fast-cooling regime (ν_m/ν_c >  1) exhibits two typical energies: a peak around ν_m and a break around ν_c. Before the peak, we expect a hard photon index of α₁ = −2/3 for E < hν_c, while for hν_c < E < hν_m the emission softens with a photon index of α₂ = −3/2. The smaller the frequency ratio, ν_m/ν_c, is, the closer the two breaks are.

Historically, most of the spectral fits that lead to the Yonetoku relation were carried out using phenomenological models, which allowed for only one power law segment at low energies, with a single spectral break coinciding with the peak. In the hypothesis that GRBs prompt spectra are produced by synchrotron processes, the fit of a single-break model (such as the Band function) to a spectrum that intrinsically presents two breaks (such as a synchrotron one) can lead to two main observational biases:

• (1.)

  It appears that, when a single-break model like Band is fit to an intrinsic double-break model, it returns $E_b<E_p^{\mathrm{Band}}<E_p$, where $E_p^{\mathrm{Band}}$ is the Band peak energy, whereas E_b and E_p are break and peak energies. Therefore, the fit of a two-breaks model to a suitable spectrum leads to systematically larger E_p,z.

• (2.)

  The low-energy spectral index, α, which one obtains fitting a single-break model, is a weighted average of the two spectral indices, α₁ and α₂, before and after the break, respectively (Toffano et al. 2021). Therefore, the hardness of the spectrum depends on the proximity of E_b to E_p.

By comparing results from the Band and synchrotron models related to the single-bin synchrotron sample, we observe a similar pattern. In fact, when the synchrotron fit returns a fast-cooling spectrum, its peak energy, $E_{p,z}^{\mathrm{Sync}}$, is systematically larger than the value obtained by fitting the Band model, $E_{p,z}^{\mathrm{Band}}$. Conversely, when 1< ν_m/ν_c <  3, the two frequencies are very close and break and peak energies are hardly distinguishable. The Band and synchrotron fits return similar peak energies; that is, $E_{p,z}^{\mathrm{Sync}}/E_{p,z}^{\mathrm{Band}}\sim 1$ (Fig. 4, left panel).

Fig. 4. Ratio between the rest-frame peak energy from synchrotron model fits, $E_{p,z}^{\mathrm{Sync}}$, and Band model fits, $E_{p,z}^{\mathrm{Band}}$, as a function of $E_{p,z}^{\mathrm{Band}}$ (left-hand panel) and low-energy break, α, histograms (right-hand panel), both relative to GRBs in the synchrotron single-bin sample. Parameters from GRBs associated with fast-cooling spectra (νm/νc >  3) are shown in red, while the ones associated with intermediate-cooling spectra (1< νm/νc <  3) are shown in blue.

In Fig. 4 (right panel), we show the distributions of the photon index, α, associated with the fast-cooling and intermediate-cooling spectra. While fast-cooling spectra span a large range of values pointing to softer spectra, the intermediate-cooling spectra exhibit harder photon indices close to α = −0.67, which is the value expected for E < E_b, z = hν_c, z.

We can conclude that the effect of fitting a single-break model such as the Band function to a synchrotron spectrum leads to systematically lower values of E_p,z. These are the premises by which the Yonetoku relation was vastly studied in the literature. According to these findings, it is not surprising that when a physical synchrotron model is used in testing the Yonetoku relation, the latter holds only for GRBs where 1< ν_m/ν_c <  3, namely when peak and break energies nearly coincide, and therefore Band fits return correct E_p,z estimates. This suggests that, at least for the GRBs that can be correctly described by an idealized synchrotron model, the underlying physical relation that holds is the ν_c, z − L_iso one. When spectra are analyzed with the Band function, $E_{p,z}^{\mathrm{Band}}$ is not far from hν_c, z, returning the notorious Yonetoku relation. However, when a second break is accounted for, the Yonetoku relation does not hold anymore, and the new ν_c, z − L_iso appears to be the underlying relation.

In Fig. 5, we show an example of this phenomenon for GRB160625B, one of the brightest sources in the multiple-bin sample. Its time-resolved analysis shows that, during the time bins in which the spectrum is in the fast-cooling regime, the corresponding E_p,z and L_iso are outliers of the Yonetoku relation. As soon as its spectrum shows an intermediate-cooling state, its E_p,z approaches the relative ν_c, z and the corresponding E_p,z and L_iso follow the Yonetoku relation.

Fig. 5. Fermi/GBM (8–900 keV) light curve and time bins of the associated time-resolved analysis of GRB160625B (upper panel). In the lower panel, we show the relative Ep,z−Liso plane. The color scale is associated with the time bin in which the spectral analysis was performed. We show, for this GRB, the results of the synchrotron fit, and we represent the bins with fast-cooling spectra as triangles, and the bins with intermediate-cooling spectra as squares. The straight black line represents the Yonetonu relation obtained from Yonetoku et al. (2004).

Nonetheless, Guiriec et al. (2013, 2015a) discuss an opposite scenario. They fit a two-components model to time-resolved GRB prompt spectra. The presence of a second peaked component at lower energies resembles (in the photon flux representation) a total spectrum with two breaks. Despite finding systematically larger E_p,z with respect to single-break model fits, they claim a tighter Yonetoku relation when a second component is accounted for in the fit. This behavior was proven in a sample of five GRBs (e.g., Guiriec et al. 2015b). However, the extension of this to larger samples would contrast with the findings of E_p,z−L_iso relations in numerous samples when single-break models are used. In fact, if both the E_p,z increase and the Yonetoku relation from single-break models constitute valid observational evidence, the fit of a second thermal component (or a double-break model) would systematically bring GRBs out of the relation. In this scenario, the only way to have a tighter relation would be for these GRBs exhibiting breaks to be all systematically beneath the best-fit line; in other words, to be scattered away from the best-fit relation with systematically lower E_p,z for the same L_iso. Nonetheless, this hypothesis is hard to test in a small sample. In this work, testing a large sample of GRBs, we instead observed that GRBs in synchrotron subsamples do not preferentially lie underneath the best-fit line, and this leads to the absence of the Yonetoku relation in both our samples.

It is important to stress that the results derived from the single-bin sample are obtained under the assumption of synchrotron being the main radiative process producing prompt emission. This hypothesis is verified through model comparison only for 21/74 GRBs. On the other hand, as is also discussed in Toffano et al. (2021), the comparison between physical and phenomenological models is quite tricky, the peculiar synchrotron features being lost in low signal-to-noise spectra, often favoring a Band function. However, the conclusions that we draw in this analysis are valid as long as the synchrotron model provides a good fit to the data, hence returning reliable estimates of both ν_m and ν_c. Nevertheless, the results drawn from the multiple-bin sample go beyond this assumption, selecting only synchrotron spectra through model comparison.

In this regard, one important result of this work is the fact that the synchrotron model we used to fit the data is not able to describe all the spectra well. The interpretation can be twofold: it appears that prompt dissipation sites and/or radiative processes are not the same for all the GRBs (therefore synchrotron emission cannot explain all GRB spectra), or that the synchrotron model we use is too simplistic to account for the GRB spectra.

From the single-bin Band subsample (Table A.2), we notice that the spectra that reject the synchrotron model show a relatively hard spectral index, α. The presence of these hard prompt emission spectra are in tension with a synchrotron interpretation, pointing to other emission mechanisms possibly produced in the optically thick region. However, we performed an additional fit of GRB spectra in the single-bin Band subsample, restricting the dataset between 8–30 keV and fitting a power law model. We excluded from this analysis GRB190114C and GRB090902B, since their low-energy spectrum at a few keV is dominated by a second power law component.

We observe that the spectral indices obtained through the power law fit, α_PL, is much smaller (therefore the spectra are softer) than the ones obtained through Band fits, α_Band (Fig. 6). In particular, the α_PL values obtained are close to the value of α = −0.67, expected from synchrotron spectra at energies of E < E_b. The reason for this mismatch could be the fact that the photon index, α, of the Band function is degenerate with the smoothness and position of the spectral peak. This finding suggests that α may be a biased indicator of GRB hardness (Burgess et al. 2015).

Fig. 6. Low-energy spectral index, α, obtained by fitting a Band model (in red) or a power law model considering data between 8–30 keV (in blue). The fit was performed to GRBs in the single-bin Band subsample. The dashed vertical line corresponds to α = −0.67, the value expected in synchrotron spectra for E < Eb. Errors are reported with a 68% confidence level.

Still, the synchrotron model used in this work is not able to provide a good representation of the data for 13 GRBs. We propose that this fact is due not to the violation of the single electron spectrum, as is often stated in the literature, but rather the shape and smoothness of the spectral profile around the peak. In these regards, there are many processes not accounted for in our modeling, such as inverse Compton scattering in the Klein-Nishina regime (Derishev 2009; Bošnjak et al. 2009; Daigne et al. 2011) provided by decaying magnetic fields (Pe’er & Zhang 2006; Derishev 2007; Zhao et al. 2014; Daigne & Bošnjak 2024), anisotropic pitch angle distribution (Medvedev 2000; Sobacchi et al. 2021), and so on. Accounting for all these processes may lead to spectra more in accordance with the observations, without the need to introduce different radiative processes.

When instead a simple synchrotron model also describes the data well, we see a particular trend. All the GRBs in the synchrotron subsample are in (marginally) fast-cooling regime, as was expected from other works (Ravasio et al. 2019).

Interestingly, the majority of them are in the intermediate-cooling regime, which is a specific stage of a marginally fast-cooling regime in which the two main frequencies, ν_m and ν_c, are coincident. Nearly 60% of GRBs in the single-bin synchrotron sample are in this state, at odds with the ∼25% of GRBs in the multiple-bins synchrotron sample. This suggests that this state is reached more frequently in the brightest pulses.

5. Conclusions

In this work, we have performed spectral analysis of GRBs in two samples: the single-bin sample (with 74 GRBs, one spectrum with 1s exposure per GRB during their peak) and the multiple-bin sample (with 13 bright GRBs, multiple spectra with exposures of a few seconds per GRB for a total of 125 spectra).

In both these samples, we fit two main models. The first one is the Band function, which phenomenologically describes GRB spectra without any assumption about the underlying physics. In addition, we used a simple physical model that accounts for power law-distributed relativistic electrons cooling down through synchrotron radiation.

We aim to test the validity of the E_p,z−L_iso (Yonetoku) relation using both Band and synchrotron models. While the former was the one historically used to derive GRB empirical relations, the introduction of a synchrotron model, which is a viable explanation for GRB prompt spectra, allows us to better investigate the nature of this emission and the relative empirical relations. We summarize the main results of this work in the following:

1. We show that most of the GRBs in the single-bin sample are well fit by the synchrotron model (58/74). However, there are GRBs whose spectrum cannot be described by a synchrotron model, preferring instead a Band function (16/74). Nonetheless, GRBs that prefer a Band model, when fit with a power law at E < 30 keV, exhibit a photon index of α_PL ≲ −0.7; in other words, they do not violate the single electron spectrum.

2. It appears that all the GRBs for which we report results from the synchrotron fit are in a fast-cooling regime (ν_m/ν_c >  1). The majority of these ones show particular ν_m/ν_c values, where ν_m (i.e., the spectral peak) is only slightly larger than ν_c (i.e., the synchrotron break). We refer to these cases as GRBs in an intermediate-cooling state, a particular case of marginally fast-cooling spectra where the break and the peak energies are basically coincident (1< ν_m/ν_c <  3).

3. We show that the fit of the Band model (single-break spectrum) to a possibly intrinsic synchrotron spectrum (double-break spectrum) leads to systematically lower E_p,z and the proximity of E_b, z to E_p,z influences the hardness of the spectrum; namely, the low-energy spectral index, α, fit with the Band model.

4. We observe that GRBs in the single-bin sample do not follow the Yonetoku relation unless they are in an intermediate-cooling state; that is, the peak and break energies are extremely close to each other.

5. We find a novel tight relation between the rest-frame cooling frequency, ν_c, z, and the isotropic-equivalent luminosity:

   $$\frac{{\nu }_{c,z}}{100\mathrm{keV}}=(0.97\pm 0.14)\times {\left(\frac{L_{\mathrm{iso}}}{{10}^{52}\mathrm{erg}/\mathrm{s}}\right)}^{(0.53\pm 0.06)}.$$

   This relation is found in both single-bin and multiple-bin samples, with slightly different slopes.

The underlying assumption of this work is that prompt emission spectra are produced through synchrotron cooling of power-law distributed particles. Nonetheless, in the multiple-bin sample analysis we report synchrotron fit results whenever the model comparison does not discard the synchrotron model. This points toward the physical nature of the ν_c, z − L_iso relation.

We show that, if GRB prompt emission is produced through synchrotron radiation, the fit using a Band function does not always represent a good spectral description. For example, a Band fit would return a peak energy estimate, $E_{p,z}^{\mathrm{Band}}$, in between hν_c, z and $h{\nu }_{m,z}=E_{p,z}^{\mathrm{Sync}}$. This suggests that, for synchrotron spectra, the fundamental physical relation is between L_iso and ν_c, z. The Yonetoku relation, which was found through Band fits, can be the effect of a single-break model fit to a double-break spectrum. Given the occurrence of GRB spectra in intermediate-cooling regimes (or in general in marginally fast-cooling regimes), $E_{p,z}^{\mathrm{Band}}$ provides a proxy of ν_c, z. The higher the intrinsic ratio, ν_m/ν_c, the higher the scatter of that GRB relative to the Yonetoku relation.

This novel interpretation of the Yonetoku relation constitutes a step forward in understanding the physics behind empirical relations and prompt emission. For the prompt synchrotron spectra, we found a new tight relation that connects the cooling regime with the total luminosity of the burst: the less it cools, the brighter it is.

This behavior is counterintuitive. One would expect, in the optically thin synchrotron framework, more radiative output when particles cool down, interacting with stronger magnetic fields. This new relation can constitute a benchmark in testing GRB prompt emission models, and push forward the possibility of using GRBs to constrain cosmological properties. The recent launch of the SVOM and Einstein Probe missions can provide, in synergies with past observatories, a variety of GRB observations with measured redshift and extend the empirical correlations study with a particularly soft burst emitting in the X-rays.

Data availability

Fermi/GBM raw data are publicly available in the Fermi GBM Burst Catalog: https://heasarc.gsfc.nasa.gov/W3Browse/fermi/fermigbrst.html. Fermi/LAT Low Energy (LLE) data are publicly available in the Fermi LAT Low-Energy Events Catalog: https://heasarc.gsfc.nasa.gov/W3Browse/fermi/fermille.html. The result of the single-bin and multiple-bin samples, together with the spectral fit plots, are available at the following link: https://zenodo.org/records/14205680

Acknowledgments

The authors thank B. Banerjee, D. Bjørn Malesani, M. Branchesi, A. Celotti, G. Ghirlanda, G. Ghisellini, L. Nava and O.S. Salafia for the fruitful discussions about this work. AM thanks the Cosmic Dawn Center in Copenhagen and the Observatory of Brera in Merate for the hospitality during the development of this project.

References

Appendix A: Single-bin sample tables

Appendix B: Multiple-bin sample tables

All Tables

All Figures

Fig. 1. Yonetoku relation obtained from the single-bin sample. Black stars represent GRBs from the Band subsample, while colored dots represent GRBs from the synchrotron subsample. The color scale is associated with the frequency ratio, νm/νc, derived from the synchrotron fit. We fit a power law to the whole single-bin sample (left-hand panel) and to the GRBs in the synchrotron subsample in an intermediate cooling regime (right-hand panel). The straight black line represents the best-fit line from the linear fit, and the green-shaded area the relative 3σ scatter region. For comparison, the dashed black line represents the best-fit line from Yonetoku et al. (2004).

In the text

	Fig. 2. νc, z − Liso relation obtained from the synchrotron single-bin sample. The color scale is associated with the frequency ratio, νm/νc, derived from the synchrotron fit. The straight black line represents the best-fit line from the linear fit, and the green-shaded area the relative 3σ scatter region.
In the text

Fig. 3. νc, z − Liso relation obtained from the synchrotron GRBs in the multi-bin sample. We show multiple time bins from different GRBs with circles of the same colors. The number inside the circles represents the bin number that data point is associated with. The straight black line represents the best-fit line from the linear fit, while the dashed-line the best-fit line obtained from the νc, z − Liso relation in the single-bin sample. The green-shaded area shows the 3σ scatter region of the relation.

In the text

Fig. 4. Ratio between the rest-frame peak energy from synchrotron model fits, $E_{p,z}^{\mathrm{Sync}}$, and Band model fits, $E_{p,z}^{\mathrm{Band}}$, as a function of $E_{p,z}^{\mathrm{Band}}$ (left-hand panel) and low-energy break, α, histograms (right-hand panel), both relative to GRBs in the synchrotron single-bin sample. Parameters from GRBs associated with fast-cooling spectra (νm/νc >  3) are shown in red, while the ones associated with intermediate-cooling spectra (1< νm/νc <  3) are shown in blue.

In the text

Fig. 5. Fermi/GBM (8–900 keV) light curve and time bins of the associated time-resolved analysis of GRB160625B (upper panel). In the lower panel, we show the relative Ep,z−Liso plane. The color scale is associated with the time bin in which the spectral analysis was performed. We show, for this GRB, the results of the synchrotron fit, and we represent the bins with fast-cooling spectra as triangles, and the bins with intermediate-cooling spectra as squares. The straight black line represents the Yonetonu relation obtained from Yonetoku et al. (2004).

In the text

	Fig. 6. Low-energy spectral index, α, obtained by fitting a Band model (in red) or a power law model considering data between 8–30 keV (in blue). The fit was performed to GRBs in the single-bin Band subsample. The dashed vertical line corresponds to α = −0.67, the value expected in synchrotron spectra for E < Eb. Errors are reported with a 68% confidence level.
In the text