© The Authors 2024

1. Introduction

Long-lasting gamma-ray bursts (GRBs) are thought to be produced by the collapse of the core of some hydrogen-stripped massive stars, in which a relativistic jet is launched by the newborn compact object, either a neutron star (NS) or a black hole (BH), and makes it to the stellar photosphere (Woosley 1993; Paczyński 1998; MacFadyen & Woosley 1999; Yoon & Langer 2005). Today, almost 30 years since the first discoveries of the long-lasting multiwavelength counterpart (from gamma rays, X-rays, and optical light to radio waves) following the initial gamma-ray emission (so-called afterglow), the wealth of knowledge of the GRB phenomenon, such as the energetics, the jet structure, its Lorentz factor, the density of surrounding environment, and the properties of the progenitor stars has grown remarkably (see e.g. van Eerten 2018; Salafia & Ghirlanda 2022; Levan et al. 2016 for reviews). However, the kind of mechanism and energy reservoir powering the GRB prompt emission, the nature of the dissipation process into gamma-rays, and the distance from the progenitor star at which it takes place, remain highly debated and hot topics (see Pe’er 2015; Zhang 2018 for reviews).

Even inexperienced eyes are caught by the variety of prompt emission light curves (LCs), which range from simple and smooth single pulses all the way up to very complex spiky profiles including many peaks and occasionally quiescent times in between, during which the signal can temporarily drop below instrumental sensitivity. Ever since, attempts to interpret the variability and apparent lack of systematic temporal evolution of the most variable LCs have resulted in mechanisms based on internal dissipation of some kind of energy (either kinetic or magnetic) into gamma-rays through shocks within a wind of relativistic shells (see e.g. Kumar & Zhang 2015 for a review). However, the rich wealth of information hidden in the observed complexity of GRB LCs remains mostly unintelligible. In particular, an open question is whether this variety could be the result of a common stochastic process. If positive, a precise characterisation of this process would reveal how GRB inner engines work, shedding light on the nature of the powering mechanism and on the nature of the engines themselves. Under this assumption, a way to gain clues is by studying the distributions of a number of properties, such as the duration, energy, and luminosity of individual peaks that make up GRB LCs, as well as their waiting times. Possible evidence for self-organised criticality was found and discussed in analogy with other astrophysical transient sources, such as solar flares or magnetar bursts, which might also be ruled by similar kinds of instabilities (Wang & Dai 2013; Lyu et al. 2020; Li et al. 2023; Li & Yang 2023; Maccary et al. 2024).

The distribution of the number of peaks within long GRB LCs has been studied and modelled in detail for the first time only recently (Guidorzi et al. 2024, hereafter G24). This distribution was studied in four different experiments: the Burst And Transient Source Experiment (BATSE; Paciesas et al. 1999), on board the Compton Gamma-Ray Observatory; the Burst Alert Telescope (BAT; Barthelmy et al. 2005) on board the Neil Gehrels Swift Observatory; the BeppoSAX Gamma-Ray Burst Monitor (GRBM; Frontera et al. 1997); and the Insight-HXMT High Energy instrument (HE; Liu et al. 2020). The result, which emerged from the analysis of these four independent catalogues, is that the distribution of the number of peaks is modelled by a mixture of two exponentials (M2E). The two components are characterised by two different average numbers of peaks per GRB: 2.1 ± 0.1 and 8.3 ± 1.0, nicknamed ‘peak-poor’ and ‘peak-rich’, respectively. The result was unexpected, given also the diversity of the four experiments in terms of energy passbands and effective areas. The existence of two distinct components was interpreted as evidence of two correspondingly different dynamical regimes, through which long GRB inner engines work.

In their work, G24 could not extend the same analysis to the entire GRB catalogue of the Fermi Gamma-ray Burst Monitor (GBM; Meegan et al. 2009), but reported an analogous preliminary result obtained over a sample of nearly 400 GRBs. In the present work we report the same kind of analysis applied to the full GRB catalogue of GBM currently available. The importance of this additional investigation cannot be overstated: not only does it represent a further test on independent data, but it is also statistically very sensitive, given that the GBM data set outnumbers the data sets previously analysed by G24. In addition, since GBM and BAT share the largest number of GRBs, it is possible to further explore to what extent the number of peaks of any given GRB depends on the detector.

The paper is organised as follows. Section 2 describes the data sample and reduction. Section 3 presents the results of the statistical analysis, whose discussion and conclusions are reported in Sect. 4.

2. Data analysis

2.1. Data set

The GBM is composed of 12 NaI scintillators sensitive in the 8–1000 keV energy range and of two BGO scintillators, sensitive to the high energy tail of a typical GRB spectrum (150 keV–30 MeV), overlapping with NaI detectors at low energies. We started from 3091 long GRBs detected by GBM from 14 July 2008 to 4 February 2024. We excluded short GRBs by discarding those with T₉₀ < 2 s, where T₉₀ is the time interval encompassing from 5% to 95% of the total fluence, as well as those with clear evidence of a compact merger origin, such as 211121A and 230307A¹ (Levan et al. 2023; Troja et al. 2022; Gompertz et al. 2023; Levan et al. 2024). Very bright GRBs that saturated the GBM detectors, such as 130427A and 221009A (Preece et al. 2014; Burns et al. 2023), were also left out. From the remaining GRBs, 154 have a known (spectroscopic) redshift, and 410 are in common with BAT.

Gamma-ray bursts affected by the simultaneous occurrence of a solar flare were discarded. Solar flares were identified by observing the relative intensity between the flare and the GRB on different units and on different energy ranges, taking into account the direction of the Sun. Being spectrally softer than GRBs, most solar flares are barely visible above 40 keV.

We finally rejected all the GRBs whose LC was lacking data points within the T₉₀ interval. Eventually, we ended up with 2971 GRBs, which is hereafter referred to as the GRB sample used for the analysis.

2.2. Background subtraction

Gamma-ray burst LCs were extracted in the full NaI energy range (8 − 1000 keV), using a 64 ms bin time. We subtracted the background using the GBM data tools² (Goldstein et al. 2022) following standard prescriptions: we selected one time window preceding the burst and another one following, both having comparable or longer duration than that of the interval containing the burst³, and interpolated the background with a polynomial function of order up to 3. For each GRB we looked into the ‘scat detector mask’ entry on the HEASARC catalogue⁴ to identify the detectors used by the GBM team. We used the TTE (Time Tagged Event) data whenever they covered the whole GRB, from the start of its T₉₀ interval to the end.

We removed the spikes caused by charged particles hitting the detectors as follows: all the bins whose counts exceeded the neighbouring counts by ≥9σ were tagged as spike candidates. Once visual inspection of different GBM units confirmed the particle, and thus the spurious nature of a given spike candidate, its counts were replaced with the mean of the neighbouring bins.

The quality of the background subtraction was assessed by computing mean and standard deviation of the normalised residuals⁵ of the time windows used for background interpolation. In particular, we made sure that the background-subtracted counts satisfied the uncorrelated Gaussian noise assumption required by MEPSA (Guidorzi 2015; see Sect. 2.3) first by carrying out a normality test on the normalised residuals (Kolmogorov-Smirnov), to ensure that the noise was Gaussian, and then by examining the autocorrelation function (ACF) of the normalised residuals and by carrying out a runs test, to ensure they were uncorrelated.

2.3. Identification of peaks

In line with G24, all the 64 ms background-subtracted LCs were sifted with MEPSA, a well-calibrated code tailored to identify statistically significant GRB peaks in uniformly sampled time series affected by Gaussian uncorrelated noise. We discarded peaks with signal-to-noise ratio S/N < 5. We ended up with 9625 peaks from 2954 GRBs having at least one significant peak. In parallel, we repeated the same selection by imposing two further thresholds on S/N: 8822 peaks from 2790 GRBs for S/N > 7; 7371 from 2388 GRBs for S/N > 9.

3. Results

3.1. Fermi-GBM sample

We computed the distribution of the number of peaks per GRB, following the procedure by G24, and tested the various models considered therein (simple exponential, simple and broken power law, and stretched exponential). As did G24, we found that the data can only be satisfactorily fitted by a M2E model given by

$$\begin{array}{c}\hfill f(n)=k(e^{-n/n_1}+\xi e^{-n/n_2}),\end{array}$$(1)

where k is a normalisation constant, n_i is the characteristic number of peaks per GRB of the ith component (i = 1, 2), and ξ is a relative normalisation parameter. The expected number of peaks per GRB of the ith component is given by

$$\begin{array}{c}\hfill ⟨n^{(i)}⟩=\frac{{\sum }_{n=1}^{+\infty }n{e}^{-n/n_i}}{{\sum }_{n=1}^{+\infty }{e}^{-n/n_i}}=\frac{1}{1-e^{-1/n_i}},\end{array}$$(2)

and the fraction of GRBs contributing to the second component is given by

$$\begin{array}{c}\hfill {\overline{w}}_2=k\xi {\displaystyle \sum _{n=1}^{+\infty }}{e}^{-n/n_2}.\end{array}$$(3)

The observed distribution along with the best fit model are shown in Fig. 1, while Table 1 reports the best fit values of the model parameters and the corresponding p-value of the χ² test.

Fig. 1. Distribution of the number of peaks per GRB for the Fermi/GBM catalogue. The observed and expected counts are displayed in orange and blue, respectively. The histogram bins were grouped to ensure a minimum number of 15 expected counts per bin. The dashed black line shows the best-fitting model of mixture of two exponentials. The blue and grey regions respectively show the 50% and 90% model confidence interval obtained from sampling the posterior distribution of the model parameters computed through Markov chain Monte Carlo simulations.

The M2E model fits the data very well (e.g. p-value of the two-tail χ² test of 0.80), while other models were unable to describe it (e.g. a fit with a stretched exponential yielded a p-value of 0.0026), thus confirming the previous results. The mean number of peaks for the peak-poor and for the peak-rich component is $⟨n^{(1)}⟩=2.{10}_{-0.10}^{+0.10}$ and $⟨n^{(2)}⟩=7.{61}_{-0.84}^{+0.97}$, respectively. The fraction of peak-rich GRBs is ${\overline{w}}_2=0.{21}_{-0.04}^{+0.04}$.

3.2. Swift-BAT/Fermi-GBM common sample

Fermi/GBM shares 410 GRBs with Swift/BAT, 361 of which were covered by BAT in burst mode. As G24 did for the common sample between CGRO/BATSE and BeppoSAX/GRBM, here we used the BAT-GBM common sample to study how the number of peaks of any given GRB depends on the detector; the ultimate goal was to test the robustness of our result.

Following G24, we considered the difference Δn = n_GBM − n_BAT between the number of peaks detected by GBM and that of BAT for each common GRB. We replicated the same analysis assuming the three different thresholds on S/N: S/N > 5; S/N > 7; and S/N > 9. Figure 2 shows the distributions of Δn for the three different thresholds on S/N. The median value of Δn is zero for all S/N thresholds. The interval |Δn|≤2 collects from 85 to 88% of the complete common samples, whereas from 92 to 94% common GRBs have |Δn|≤3. These values prove that both experiments are equivalently sensitive to the number of peaks, in spite of their different energy passbands. As noted by G24, Δn for most GRBs is significantly smaller than the difference between the mean values of peak-rich and peak-poor GRBs, $⟨n^{(2)}⟩-⟨n^{(1)}⟩=5.5_{-0.7}^{+0.9}$ (90% confidence), which makes the identification of the two families of long GRBs robust and essentially detector-independent.

Fig. 2. Distribution of Δn for three different thresholds on S/N: 5, 7, and 9.

3.3. Combined sample

We merged the results of the five catalogues and treated them as if they were a single result. The GBM result reported in the present work, is clearly dominant as its 2954 GRBs outnumber those from the other catalogues (which included 1457, 1277, 820, and 202 GRBs). We obtained an acceptable fit (p-value of 0.06) with the same but more accurate value for the fraction of peak-rich GRBs: ${\overline{w}}_2(\mathrm{all})=0.{20}_{-0.02}^{+0.02}$, and similar values for the mean number of peaks for the peak-poor and peak-rich component $⟨n^{(1)}⟩=2.{25}_{-0.07}^{+0.08}$ and $⟨n^{(2)}⟩=9.{08}_{-0.66}^{+0.72}$, respectively.

3.4. Impact of selection effects and the Malmquist bias

The number of peaks per GRB clearly depends on the S/N of each peak. As such, the selection of GRBs with suitable peaks as well as their classification of peak-poor versus peak-rich depends on S/N, too. Any given GRB that would have been observed at higher-than-actual redshift, could have been either tagged with fewer peaks or missed altogether. Consequently, our results are potentially affected by a significant Malmquist bias. As we also did in G24, repeating the analysis for a range of different S/N thresholds, and thus mimicking the effect of moving given GRBs further away, our aim was to explore the robustness of our results against any distance-related bias. Two interesting facts emerge from this analysis:

1. The fraction ${\overline{w}}_2$ of peak-rich GRBs essentially remains constant for all catalogues and for the different S/N thresholds, and therefore appears to be a robust property that is least affected by anything that depends on peak brightness. In the following, we delve further into this aspect, which seems to clash with the Malmquist bias;

2. The average number of peaks per GRB changes only slightly for different thresholds on S/N: from 2.1 (S/N > 5) to 1.8 (S/N > 9) for the peak-poor class, and, correspondingly, from 7.6 to 7.2 for the peak-rich class (Table 1). This is understood by the following argument: not only does the number of peaks decrease by increasing the threshold on S/N of any given GRB (mimicking the consequence of moving it to higher redshifts), but the number of GRBs featuring at least one significant peak also decreases. As a result, the average number of peaks per surviving GRB decreases less abruptly than the overall number of peaks.

In principle, a direct way to understand the impact of the Malmquist bias is by focusing on a sample of GRBs with measured redshift. Unfortunately, this sample is currently too small to make a statistically sound analysis feasible. Nonetheless, we explored in more detail the reasons for the first consideration above as follows. We addressed two questions: By moving a peak-rich GRB to progressively higher redshifts, we wanted to know (a) at what relative distance it becomes peak-poor, thus contributing to altering the true distribution, and (b) at what relative distance the GRB disappears, an event coinciding with the brightest peak sinking below threshold?

The answer to (b) depends on the intensity of the brightest peak, while we assume that the answer to (a) depends on the intensities of the third or fourth brightest peaks, once we take three or four as a reasonable proxy of the number of peaks for discriminating between peak-poor and peak-rich GRBs.

For each of the peak-rich GRBs with at least 6 S/N > 5 peaks, we took the four brightest ones. We ended up with 399 GRBs. For each GRB we calculated the three following ratios: r₂₁ ≥ r₃₁ ≥ r₄₁, respectively corresponding to the ratio of the intensities of the second, of the third, and of the fourth brightest peaks to the intensity of the brightest peak.

As a result, on average it is r₄₁ ∼ 2/3, with 0.56–0.77 as interquartile, so the range of intensities of the four brightest peaks is relatively narrow. Specifically, from the disappearing of the fourth peak to the disappearing of the brightest peak, the threshold on S/N increases only by ∼50% on average. Translating the relative change of the S/N threshold in terms of the relative change in the maximum volume within which a given peak can be detected, the volume increases by a factor of (2/3)^−3/2 ≈ 1.8, so 80% (we used the scaling V ∝ p^−3/2 of volume V on peak flux p). This is a relatively short range compared with the range of redshifts spanned by the GRB populations that make up the observed catalogues and the GBM catalogue in particular. Replacing the fourth with the third brightest peak, the relative volume change shrinks to 60%, so the previous conclusion holds true a fortiori. This property of most peak-rich GRBs may explain why the Malmquist bias does not strongly affect the result of a constant peak-rich fraction throughout the different catalogues and GRB populations.

4. Discussion and conclusions

Guidorzi et al. (2024) modelled for the first time the distribution of the number of peaks in long GRB LCs through the analysis of four independent catalogues (CGRO/BATSE, Swift/BAT, BeppoSAX/GRBM, and Insight-HXMT/HE), which differ from each other in terms of energy passbands and effective areas. Only a mixture of two exponentials was found to model all four distributions, supporting the existence of two families, called peak-poor and peak-rich, characterised by an average number of peaks per GRB of around 2.1 ± 0.1 and 8.3 ± 1.0, respectively. It was also found that peak-poor GRBs make up about 80% of the GRB populations observed by past and present experiments.

In the present work we carried out the analogous investigation over the Fermi/GBM GRB catalogue, which currently includes about 3000 GRBs and therefore offers a further opportunity to test the previous results with a high degree of statistical accuracy. Remarkably, not only did we confirm that the distribution can only be modelled by a mixture of two exponentials, but we also found that the best-fitting parameter values are fully consistent with what was found by G24 from different and independent data sets. Specifically, the expected numbers of peaks for the two families of peak-poor and of peak-rich GRBs are remarkably similar to the G24 results. Analogously, the fraction of peak-rich GRBs is found to be 0.21 ± 0.04, fully consistent with what was found by G24. Although they were ignored in the present analysis, long-merger candidates would belong to the hard tail of the peak-rich family: 211211A has 48 peaks and 230307A has 53 peaks. Should peak richness be confirmed in future analogous events, it would provide a clue to the way the inner engines of these mergers work.

Figure 3 shows a visual comparison between the GBM results reported here and the results obtained by G24. Interestingly, while the fraction of peak-poor GRBs is basically the same for all the five catalogues, the average numbers of peaks per each GRB for both families of the only BATSE is somewhat higher. As discussed in G24, this is ascribed to the larger effective area of BATSE compared with the other instruments, including GBM.

Fig. 3. This figure describes our present results compared with those of G24. Top panels, left to right: violin plots of the posterior distributions of the three parameters ⟨n(1)⟩, ⟨n(2)⟩, and ${\overline{w}}_2$, respectively representing the expected number of peaks per GRB of the peak-poor and the peak-rich classes, and the fraction of peak-rich GRBs obtained for the GBM sample. The horizontal bars span 5–95% quantiles. Bottom panels, left to right: comparison between the results obtained for the same parameters in this work, with those obtained by G24 from other catalogues.

In conclusion, we confirmed that the existence of these two families is a robust result that is not sensitive to the experiments and their energy passbands. As discussed by G24, these two different dynamics ruling long GRB prompt emission might point to either two different families of engines left over by the collapse of the core, either a NS or a BH, or two different regimes in which GRB engines release their energy. As noted in G24, the distinctive property of peak-rich GRBs is the presence of sub-second variability that adds to the slow-varying component, which is instead observed in both kinds of GRBs.

The existence of two distinct dynamics in the observed population of GRBs can be interpreted within different models of how dissipation into gamma-rays takes place (see G24 for a detailed discussion). Consequently, it does not directly help to discriminate between different interpretations. Nevertheless, it provides a solid clue to how GRB engines work, with the possibility that they could operate close to a critical regime. In this way, a completely different dynamical behaviour can result from a small difference in some key properties, such as the degree of magnetisation of the relativistic ejecta, as discussed in G24. Should the peak richness of GRBs be found preferentially associated with other key properties in future investigations, the implications would help to narrow down the identity of GRB engines and the nature of the dissipation process that governs the prompt emission.

Acknowledgments

We thank the anonymous reviewer for providing insightful comments that improved the manuscript. R.M. and M.M. acknowledge the University of Ferrara for the financial support of their PhD scholarships. R.M. acknowledges the Fermi team for covering flight costs for attending the Fermi summer school 2023. A.T. acknowledges financial support from ASI-INAF Accordo Attuativo HERMES Pathfinder operazioni n. 2022-25-HH.0 and the basic funding program of the Ioffe Institute FFUG-2024-0002. L.A. acknowledges support from INAF Mini-grant programme 2022. A.C. is partially supported by the 2023/24 “Research and Education” grant from Fondazione CRT. The OAVdA is managed by the Fondazione Clément Fillietroz-ONLUS, which is supported by the Regional Governement of the Aosta Valley, the Town Municipality of Nus and the “Unité des Communes valdotaines Mont-Emilius”. L.F. acknowledges support from the AHEAD-2020 Project grant agreement 871158 of the European Union’s Horizon 2020 Program.

References

All Tables

All Figures

Fig. 1. Distribution of the number of peaks per GRB for the Fermi/GBM catalogue. The observed and expected counts are displayed in orange and blue, respectively. The histogram bins were grouped to ensure a minimum number of 15 expected counts per bin. The dashed black line shows the best-fitting model of mixture of two exponentials. The blue and grey regions respectively show the 50% and 90% model confidence interval obtained from sampling the posterior distribution of the model parameters computed through Markov chain Monte Carlo simulations.

In the text

	Fig. 2. Distribution of Δn for three different thresholds on S/N: 5, 7, and 9.
In the text

Fig. 3. This figure describes our present results compared with those of G24. Top panels, left to right: violin plots of the posterior distributions of the three parameters ⟨n(1)⟩, ⟨n(2)⟩, and ${\overline{w}}_2$, respectively representing the expected number of peaks per GRB of the peak-poor and the peak-rich classes, and the fraction of peak-rich GRBs obtained for the GBM sample. The horizontal bars span 5–95% quantiles. Bottom panels, left to right: comparison between the results obtained for the same parameters in this work, with those obtained by G24 from other catalogues.

In the text