Distributions of Energy, Luminosity, Duration, and Waiting Times of Gamma-Ray Burst Pulses with Known Redshift Detected by Fermi/GBM

Discovered more than 50 years ago, gamma-ray burst (GRB) prompt emission remains the most puzzling aspect of GRB physics. Its complex and irregular nature should reveal how newborn GRB engines release their energy. In this respect, the possibility that GRB engines could operate as self-organized critical (SOC) systems has been put forward. Here, we present the energy, luminosity, waiting time, and duration distributions of individual pulses of GRBs with known redshift detected by the Fermi Gamma-ray Burst Monitor. This is the first study of this kind in which selection effects are accounted for. The compatibility of our results with the framework of SOC theory is discussed. We found evidence for an intrinsic break in the power-law models that describe the energy and the luminosity distributions.


INTRODUCTION
Gamma-ray prompt emission is the first electromagnetic radiation observed from a gamma-ray burst (GRB) source.It is caused by at least two different types of catastrophic events: (i) the merger of a binary system of compact objects (Eichler et al. 1989;Paczynski 1991;Narayan et al. 1992;Abbott et al. 2017); (ii) the core collapse of fast rotating, hydrogen-stripped massive stars (dubbed as "collapsars", Woosley 1993;Paczyński 1998;MacFadyen & Woosley 1999;Yoon & Langer 2005).Most events of the first (second) class are associated with short (long) GRBs (SGRBs/LGRBs), with some notable exceptions: short with extended emission GRBs (Norris & Bonnell 2006;Gehrels et al. 2006), short GRBs coming from collapsars (Rossi et al. 2022;Ahumada et al. 2021), and long GRBs coming from mergers (Rastinejad et al. 2022;Troja et al. 2022;Gompertz et al. 2023;Yang et al. 2023).In the light of the emerging complexity, nowadays families (i) and (ii) are often referred to as merger and collapsar GRBs, or, equivalently, as Type-I and Type-II GRBs, respectively (Zhang 2006).
The nature of the central engine that powers a GRB is still debated: either a hyper-accreting stellar mass black hole (BH; Lei et al. 2013;Janiuk et al. 2017;Sado 2019) surrounded by a thick accretion disk, or a fastrotating, highly magnetized neutron star (NS; known as a "millisecond magnetar", Usov 1992;Wheeler et al. 2000;Thompson et al. 2004;Metzger et al. 2011).In the case of a hyper-accreting BH, the jet can be produced either by a neutrino-dominated accretion flow (NDAF) in which neutrinos tap the thermal energy of the accretion disk via neutrino-antineutrino annihilation, launching a thermally-dominated fireball (Popham et al. 1999;Di Matteo et al. 2002), or through the Blandford-Znajek effect (BZ, Blandford & Znajek 1977), which converts the BH spinning energy into a Poyntingflux dominated outflow.Also in the case of a msmagnetar, the rotational energy could be the source of energy (Duncan & Thompson 1992;Metzger et al. 2011).Magneto-hydrodynamic (MHD) instabilities can also be considered to launch the jet, possibly for both engines (Bromberg & Tchekhovskoy 2016).
GRB temporal profiles exhibit a remarkable variety in intensity, duration, number of pulses, and variability in general.Most profiles can be seen as a superposition of fast-rise exponential decay pulses (FRED, Fishman & Meegan 1995).Prompt emission is often characterized by phases of intense activity separated by periods of prolonged inactivity (quiescent times).
The study of the waiting times (WTs) indicates that GRB central engine could emit pulses according to a non-stationary Poisson process (Guidorzi et al. 2015).The diversity of all observed bursts could be understood as different realizations of a common stochastic process.Should this be the case, an open question is whether this process can be characterized in a more detailed way, so as to reveal the dynamics that rules the energy release.If positive, this would strongly constrain the way GRB engines work and, ultimately, help identify their nature and the powering mechanism(s).
Complex systems with many interacting elements, presenting erratic on-off intermittency (such as earthquakes, neuronal activity, stock market, evolution of species) are interpreted in the self-organized criticality (SOC) framework (Sornette & Sornette 1989;de Arcangelis 2012;Bartolozzi et al. 2005;Bak & Sneppen 1993).Invoked to explain the ubiquitousness of 1/f noise spectra, SOC was originally applied to describe sand pile avalanches (Bak et al. 1987(Bak et al. , 1988)).A SOC system can be seen as an out-of-equilibrium nonlinear dynamical system in which energy steadily brought by an external source drives the system towards a critical point, at which the energy is liberated through scale-free avalanches (Aschwanden 2014).
One of the signatures of a SOC lies in its scale invariance: the released energy, luminosity (or peak flux, depending on the context), duration, and waiting time of individual avalanches are power-law (PL) distributed, with a set of relations between the different PL indices (see Section 4).
SOC was pointed out to explain a wide range of astrophysical phenomena, such as solar flares, lunar craters, Saturn rings, auroral emission from Earth's magnetosphere, pulsar glitches, blazars, stellar flares (see Aschwanden et al. 2018 for a review; Aschwanden & Güdel 2021).One of the main accomplishments of SOC theory is the successful and exhaustive description of solar flares (Lu & Hamilton 1991;Charbonneau et al. 2001;Aschwanden & Aschwanden 2008a,b;Aschwanden & McTiernan 2010;Wang & Dai 2013).A body of evidence for SOC behavior was also found in the case of soft gamma-ray repeaters or magnetars (SGRs; Nakazato 2014), which are also GRB engine candidates.It is worth noting that SOC behavior in the case of a BH surrounded by an accretion disk was also investigated as well as for the Galactic Centre Sgr A* X-ray flaring activity (Li et al. 2015;Yuan et al. 2018).
Do GRB engines work as SOC systems?Cellular Automata (CA) simulations of 1D and 2D magnetized outflow demonstrated that such systems could reach a SOC state (Harko et al. 2015;Dȃnilȃ et al. 2015).This is consistent with the picture of prompt emission being the result of fast magnetic reconnection events occurring in a Poynting outflow powered by a highly magnetized NS.Such CA simulations were also performed in the case of a hyper-accreting BH (Mineshige et al. 1994), another leading candidate as GRB progenitor.Thus, the SOC interpretation could actually be consistent with radically different theoretical pictures of prompt emission.In the case of a hyper-accreting BH, SOC avalanches could be triggered either by thermal instability of the accretion disk (Janiuk et al. 2007;Janiuk & Yuan 2010;Taylor et al. 2011) or by magnetic instabilities (as the kink mode instability Bromberg et al. 2019;Giannios & Spruit 2006, or the tearing mode instability Del Zanna et al. 2016; Yang 2019), while only magnetic instabilities are expected in the ms-magnetar scenario.
Evidence for SOC in GRB X-ray flares was suggested by Wang & Dai (2013), whose results could be understood as a 1D SOC system with normal diffusion.Multipulsed LGRBs (Lyu et al. 2020), SGRBs (Li et al. 2023), and LGRB precursors (Li & Yang 2023) seem roughly compatible with a 3D one.In these works, GRB properties are PL distributed, even if the agreement with SOC theory is not always compelling (see Table 5).Moreover, observables such as fluence or peak flux of pulses are often considered, which might not necessarily reflect an intrinsic property of GRB engines, but strongly depend, among other things, on redshift.
In this paper we used a sample of Type-II GRBs detected by the Fermi Gamma-ray Burst Monitor (GBM; Meegan et al. 2009) with known (spectroscopic) redshift to compute the distributions of the isotropic-equivalent released energy, of the luminosity, of the duration, and of the waiting time of the individual pulses that compose GRB light curves.We performed a thorough estimation of the selection effects, splitting our sample into redshift bins and estimating the detection efficiency separately for each redshift bin.The redshifts of our data sample were determined through spectroscopic observations of either the optical afterglow or the host galaxy.
Section 2 describes the GRB samples and the data analysis, results are reported in Section 3 and discussed in Section 4, with conclusions in Section 5.In Appendix A the impact of performing a time-averaged analysis instead of a time-resolved one is discussed.

Sample selection
The GBM consists of 12 NaI scintillators sensitive in the (8 − 1000 keV) energy range and 2 BGO scintillators, sensitive to higher energies (150 keV − 30 MeV), overlapping at low energies with NaI detectors.
We first selected all the GRBs detected by the GBM from August 2008 to July 2022, with redshift measured spectroscopically, except for a few constraining photometric estimates. 1We excluded 130427A and 221009A from our analysis since the main burst saturates all NaI detectors of the GBM (Ackermann et al. 2014;Lesage et al. 2023).
To ensure that the population of GRB progenitors is as homogeneous as possible, we selected only Type-II GRBs.To this aim, we first considered events with T 90 ≥ 2 s (Kouveliotou et al. 1993).We used T 90 , after excluding the LGRBs that have been reported as credible Type-I candidates.The boundary value of T 90 = 2 s was originally derived from the catalog of the Burst And Transient Source Experiment (BATSE; Paciesas et al. 1999): this threshold depends on the energy passband of the instrument and, as such, is consequently different for the BATSE (∼ 2.4 s), the Neil Gehrels Swift Observatory BAT (Gehrels et al. 2004; 0.8 s, see Bromberg et al. 2013) and the GBM (∼ 4.2 s, see von Kienlin et al. 2020).We therefore assumed the more conservative threshold of T 90 > 4.2 s, as it is more suitable for the GBM.T 90 values were taken from the GBM catalog (von Kienlin et al. 2020).
Despite its T 90 = 2.6 ± 0.6 s, we included 141004A because two independent telescopes (the Gran Telescopio Canarias and the Reionization and Transients Infrared Camera) detected a possible brightening in the optical counterpart, that was interpreted as an emerging supernova (Littlejohns et al. 2014;Schulze et al. 2014b).Same for 200826A with T 90 = 1.14±0.13s, for which robust evidence for a collapsar origin was found (Ahumada et al. 2021;Zhang et al. 2021;Rossi et al. 2022).Among the long-lasting GRBs we excluded 211211A, for which compelling evidence for a Type-I burst was reported despite its T 90 of 34 s and a multi-peaked light curve that deceptively looks like that of a typical Type-II GRB 1 They are 120922A, 151111A, and 200829A with respectively z = 3.1 ± 0.2, z = 3.5 ± 0.3, and z = 1.25 ± 0.2.(Rastinejad et al. 2022;Troja et al. 2022;Gompertz et al. 2023;Yang et al. 2022).We also excluded two GRBs (090510 and 161129A) among the so-called short with extended-emission GRBs (hereafter SEE GRBs; Norris & Bonnell 2006), that are likely Type-I GRBs in spite of their long duration due to a spectrally soft, long-lasting tail.As a result, we were finally left with a sample of 142 Type-II GRBs.The sample is described in Table 1.

Light curve extraction and background interpolation
The GBM has 12 NaI detectors oriented so as to provide a nearly uniform coverage of the sky.Only a fraction of the 12 detectors have a good view of any given GRB, typically the ones with a small angle between the normal of the detector and the direction of the GRB.To increase the signal-to-noise ratio (S/N), one conveniently sums up the most illuminated detectors' light curves.Summing too many detectors would mainly add noise and end up with a worse S/N.To find the optimal combination, we applied the following strategy: 1) we preliminarily selected the detectors with a viewing angle θ < 60 • , as recommended by the GBM team (Bhat et al. 2016), unless no detector fulfilled the previous condition or when these detectors had been used by the GBM team to compute the T 90 (von Kienlin et al. 2020); 2) we generated all combinations of detectors and their corresponding light curves (e.g. if we selected detectors n1, n2, n3 at step 2, then we considered the different combinations n1+n2, n2+n3, n1+n3 and n1+n2+n3); 3) we took the combination with most pulse candidates (with S/N > 6).In case of equality, we opted for the combination of detectors that maximizes the total S/N of all identified pulse candidates.For each GRB, we extracted the Time-Tagged Event (TTE) light curves for the selected NaI detectors in three different energy channels: 8-1000, 30-1000, and 40-1000 keV.TTE data typically cover from −30 to 300 s with respect to the trigger time.For some extended triggers, TTE data coverage extend from −100 to 450 s.For very long GRBs (T 90 ≳ 600 s), TTE data do not cover the full extent of the events, as in the case of the ultralong 091024 (Virgili et al. 2009;Gruber et al. 2011a), so we used instead the CSPEC data2 .For the ultra-long 160625B (e.g., Zhang et al. 2018) we used TTE data for the first part of the event (t < 400 s) and CSPEC data for the second part (t > 400 s).The light curves were extracted using the GBM tools publicly available3 (Goldstein et al. 2022).We binned these light curves with a time resolution from 16 to 1024 ms.Some GRBs had particularly short timescales that required the higher resolution of 4 ms4 (see Section 2.3 for more details on the identification of the optimal bin time for each GRB).
In order to interpolate the background in the "onsource" time interval (where the burst is present), we chose two "off-source" windows and modeled the background with a polynomial function with order up to 3.
To check the quality of the interpolation, we computed the normalized residuals ϵ i in the off-source windows defined as: where r i is the count rate in the i-th bin, σ ri is the count rate uncertainty in the i-th bin, given by: obtained assuming a Gaussian regime in counts per bin time ∆t, and b i is the interpolated background count rate.An optimal background modeling is characterized by normally distributed residuals with null mean value and unity standard deviation N (0, 1) (standardized normal distribution).
The background-subtracted light curve was considered for the next steps in the analysis, provided that the following conditions on the mean value µ ϵ and standard deviation σ ϵ were fulfilled: |µ ϵ | < 0.2 and |σ ϵ − 1| < 0.2.Otherwise, the corresponding detector was ignored.
We obtained a background-subtracted light curve independently for each detector and energy passband and then we added them through all the combinations mentioned above to obtain as many background-subtracted light curves.We selected for the following study the light curve associated with the best combination of detectors.

Pulse identification
Pulses were identified using mepsa5 (Guidorzi 2015), that is a peak search algorithm designed and calibrated to find statistically significant local maxima in GRB light curves.Each pulse candidate is characterized in terms of the pulse peak time t 0 , the pulse peak rate A, the S/N, and the detection timescale at which the pulse was detected with the highest S/N, hereafter denoted ∆t det .For the goal of the present analysis we rejected all the pulse candidates with S/N < 6 to ensure a high purity.
Since mepsa applies to uniformly binned light curves, we preliminarily had to determine the optimal bin time for any given GRB.On the one hand, using an unnecessarily short bin may lead to relatively numerous statistical fluxes.On the other hand, a too coarse resolution would wash out genuine temporal structures with a consequent loss of information.We opted for the longest bin time ∆t for which all the pulses detected by mepsa are resolved.This condition was implemented by asking ∆t ≤ ∆t det,i /2 for any generic i-th pulse candidate.As for the choice of the combination of detectors, we took the one with most pulse candidates, as explained in Section 2.2.
We applied mepsa to the light curves in the three energy passbands mentioned in Section 2.1 and with a bin time in the range 16-1024 ms, with a maximum rebin factor of 100 (see Guidorzi 2015 for details), except for the few GRBs which required a 4-ms resolution (see Section 2.2).In harder channels, pulses tend to be narrower (Fenimore et al. 1995): we exploited this property to improve our ability in identifying and separating partially overlapped pulses.
Since most photons are relatively soft due to GRB spectral shapes, neglecting the softer energy channels turns into a statistically poorer S/N.We therefore chose to merge the results obtained with mepsa on the three different energy ranges, whenever visual inspection suggested the presence of blended pulses that were missed by mepsa in the analysis of individual energy channels.We made sure to count just once the pulses that were identified in more than one energy band.To this aim, we adopted the following strategy: let t 0i and ∆t det,i be respectively the pulse peak time and the detection timescale of the i-th pulse candidate detected by mepsa in the 8-1000 keV band.Let t ′ 0j and ∆t ′ det,j be the analogous quantities of the j-th pulse candidate from one of the harder channels.The two pulse candidates are tagged as distinct ones if their intervals do not overlap: (3) The condition (3) eliminates most of multiple detections of the same pulse, except for a few cases, which were corrected after visual inspection.

Light curve modeling
We modeled each identified pulse candidate with a FRED profile, which was found to describe most GRB pulses (Norris et al. 1996): where A is the pulse peak rate, t 0 the pulse peak time, τ r the pulse rise time, ξ the decay-to-rise ratio (note that the pulse decay time is τ d = ξτ r ), and ν the so-called peakedness, which determines the pulse sharpness.As first noted by Norris et al. (1996) in the analysis of BATSE bursts, ξ mostly ranges between 2 and 3, whereas ν lies between 1 (pure exponential) and 2 (Gaussian).Finally, we investigated the impact of the choice of the pulse model by alternatively adopting the one proposed in Norris et al. (2005).In this model two parameters (two timescales) instead of three are needed to define the shape of a pulse.
We used a non-linear least square algorithm 6 to fit the light curves with a superposition of Norris pulses, as in Equation ( 5): The number of pulses N p was given by the results of mepsa.We used the pulse peak times and the pulse peak rates given by mepsa as starting values for the fit.We set some boundaries on the 5N p parameters (5 for each pulse).The pulse peak time of a generic i-th pulse was constrained to be in the interval [t 0i − ∆t det,i /2 ; t 0i + ∆t det,i /2].The pulse peak rate was allowed to vary within the interval where σ Ai is the pulse peak rate error given by mepsa, whereas n 2 can be adjusted by the user.For well separated pulses we used n 2 = 1.Instead, for partially overlapped pulses we allowed larger values for n 2 .For the GRB light curves consisting of a forest of overlapped pulses, we left the pulse peak rate unconstrained.To avoid unrealistic modeling of pulses, we constrained ξ and ν in the following ranges: 0.1 < ξ < 10 and 0.4 < ν < 4, as suggested by the results by Norris et al. (1996) in their systematic analysis of BATSE GRBs, which shared the energy passband with GBM NaI detectors.
Similarly to the validation of the background modeling (Section 2.2), we here verified the quality of the fit by computing the corresponding normalized residuals defined as: where r i is the i-th count rate, σ ri its associated uncertainty, f (t i ) the count rate predicted by the model of Equation ( 5) at time t i .
We computed the mean, median, and standard deviation of {ϵ i } over the same time interval that was used for the background interpolation plus the one that includes the GRB profile.Ideally, ϵ i is distributed as a standardized Gaussian N (0, 1).A non-zero mean would reveal unaccounted trends in the modeling, whereas a value of σ > 1 (σ < 1) would indicate under(over)-fitting.A poor modeling often turned out to be due to the presence of weak and/or blended unaccounted pulses: in such cases we applied the strategy described in Section 2.3 to add pulses to the signal until we reached an acceptable solution.In some cases we had to add pulses by visual inspection, until the quality of the fit improved sufficiently.Figure 1 shows an example.For a few challenging GRBs we could only come up with a relatively poor fit.More details are reported in Section 3.
For each GRB we obtained a set of 5N p best-fit parameters µ fit with an associated covariance matrix Σ fit , which was used to estimate the uncertainties on the best-fit parameters.The uncertainties on the best fit parameters are given by the square root of the diagonal of the covariance matrix.For each GRB we randomly generated a sample of 10 3 sets of parameters around the best-fit solution according to a multivariate Gaussian distribution N (µ fit , Σ fit ) truncated along the positivedefinite parameters.Once the best-fit profile of Equation (4) is known for every pulse, we computed the integral counts of the i-th pulse: The uncertainty on C i was estimated as the standard deviation of the corresponding distribution of values generated by random realizations of the set of parameters using N (µ fit , Σ fit ).

Isotropic-equivalent energy and luminosity
In the following, energy flux (erg cm −2 s −1 ) and fluence (erg cm −2 ) are denoted with F and Φ, respectively.
We estimated the fluence of the i-th pulse as Φ i = C i f tot , where f tot = Φ tot /C tot is the conversion factor obtained from the corresponding time-integrated quantities.Φ tot is the fluence of the whole GRB and C tot the corresponding counts calculated on the same time interval.
Φ tot was estimated by taking the best fit model and its associated spectral parameters provided by the GBM GRB catalog 7 .For 091024 we used instead fluence values on the whole event provided by Gruber et al. (2011b).K-corrections, needed to estimate the fluence in the common rest-frame energy passband 1-10 4 keV, were computed in the following way: This choice complies with what is usually adopted in the GRB literature (e.g., Amati et al. 2002).Therefore, we computed Φ tot as Φ tot = k c Φ 8−1000 .This procedure assumes a negligible impact of the spectral evolution throughout the GRB: in Appendix A we show that the impact of our approximation on the resulting pulse energy distribution does not affect our conclusions significantly.
The corresponding released isotropic-equivalent energy is calculated as where d L is the luminosity distance 8 .The corresponding isotropic-equivalent luminosity is computed as where F i = A i f tot .

RESULTS
mepsa detected 974 pulses in the 8-1000 keV band, while 13 additional pulses were detected in the harder channels.We had to add at least one pulse by visual inspection for 45 GRBs.Overall, we added 159 undetected pulses to the 987 detected by mepsa, so that the final sample consists of 1146 pulses, whose 14% had to be identified through visual inspection.
We split the GRB sample in two groups, depending on the quality of the modeling and on the presence of pulses with extreme values of the parameters which define the pulse shape (ξ and ν).The Silver Sample 9 (hereafter,

SS) includes both (i)
GRBs with a standard deviation of the normalized residuals deviating from 1 by at least 0.1: |σ ϵ −1| > 0.1, and (ii) GRBs for which the best-fit model contains at least one pulse with an extreme shape in terms of parameters ξ and ν (either ξ < 0.11 or ξ > 8.99 and ν < 0.41 or ν > 3.99 respectively).The Golden Sample (hereafter, GS) includes all the remaining GRBs.
The GS (SS) consists of 119 ( 23) GRBs.In terms of pulses, the GS (SS) includes 696 (450) pulses.On average, SS GRBs have more pulses per GRB, with a median of 15 against 3 of the GS.They are also brighter and more energetic with a median fluence ∼ 4.5 times higher (4.4×10 −5 vs. 9×10 −6 erg cm −2 ).The isotropicequivalent energy of SS GRBs is on average ∼ 7 times bigger (4.3 × 10 53 vs. 6.2 × 10 52 erg).The fact that, on average, the GRBs having the best signal are also the most problematic ones to model, reveals the degree of complexity of GRB time profiles as well as the limits of our approach.
Our mean values of ⟨ξ⟩ and ⟨ν⟩ for the whole sample (GS+SS) are respectively 2.6 and 1.6, which is broadly compatible with the corresponding distribution obtained by Norris et al. (1996) for BATSE GRBs (Figure 2).

Detection efficiency
To account for the completeness limit of the energy and of the luminosity distributions, we estimated the detection efficiency as a function of the pulse peak rate A. In the case of energy, the detection efficiency is a function of the pulse counts I.
In either case we estimated mepsa efficiency to detect features on samples of simulated pulses with different shapes for given pulse peak rates/pulse counts.We carried out two complementary kinds of simulations: (i) single-and (ii) multi-pulse.
For (i), for every GRB we took its interpolated background profile and added a randomly generated FRED pulse.We finally added Poisson noise.
For (ii), for every GRB we randomly generated and added a FRED pulse to the best-fit model of the real GRB light curve plus the modeled background curve.Hence, we ended up with as many replicates of the original GRB light curves, each of which included one additional synthetic pulse.This was done in order to mimic the more realistic situation in which a given pulse might have occurred in conjunction with the other ones that were already identified.We finally added Poisson noise.
Each simulated pulse was generated assuming for t r , ξ, and ν the corresponding log-normal distributions that best fit the observed distributions of the same parameters in the real samples.These best-fit log-normal distributions were validated through a two-population Kolmogorov-Smirnov (KS) test. 10A random variable X is log-normally distributed if ln(X) ∼ N (µ, σ 2 ).The best fit log-normal distribution parameters (µ, σ) are (−0.69,1.5), (0.66, 0.97), and (0.35, 0.62) for t r expressed in s, ξ, and ν distributions, respectively.
In the case of energy, the simulated pulse peak rate A was computed as a function of a given I, which is the total counts assigned to a given pulse.The rationale behind this choice is to describe the detection efficiency for pulses with a given I, which is the most important parameter for a pulse to be detected (against a given background rate).Since I scales linearly with A (see Equation ( 7)), this is easily done: where we set A 1 = 1.A is then calculated as follows: In the case of luminosity, as pulses of different shapes are grouped by count rates, we simply assigned a given pulse peak rate to a given group of simulated pulses.For (i), the pulse peak time was set to t 0 = 0 s by convention, as its position is not relevant in the singlepulse case.For (ii), the pulse peak time was drawn from a uniform distribution with the constraint to obtain pulses with separability (defined as the distance between the simulated pulse and its nearest neighbor divided by its FWHM) between 1 and 10.
For each GRB we generated a sample of 100 simulated FRED-like pulses with parameters [t 0 , t r , ξ, ν] randomly chosen from the corresponding best-fit log-normal distributions.
We repeated the process for evenly spaced logarithmic values of I (A), in the range 10 2 -10 4 counts (10 2 -10 4 counts s −1 ).We used the same sample of FRED 10 The KS test was done using scipy.stats.ks2samp.
parameters for each value of I or A. This procedure allowed us to model the detection efficiency as a function of the pulse counts and of the pulse peak rate of a generic pulse.We finally obtained the detection efficiency as a function of either E iso or L iso , using Equations ( 9) and (10) separately for 9 bins in redshift11 : these were obtained through an evenly spaced logarithmic sampling of the luminosity distance range.
The detection efficiency η z (E iso ) (η ′ z (L iso )) was estimated as the fraction of simulated pulses with E iso (L iso ) (taken from the grid of simulated values) that had been identified by our procedure at the given redshift bin.Finally, η z (E iso ) (η ′ z (L iso )) was obtained for any value of E iso (L iso ) through interpolation.12

Modeling of energy and luminosity distributions
Having modeled the redshift-dependent selection effects that are present in the observed E iso and L iso distributions, we now aim to test if/which given conjectural intrinsic differential distribution f (E iso ) = dN/dE iso (f ′ (L iso ) = dN/dL iso ) may explain the results, once the selection effects are properly simulated.To this aim, we started with the simplest distribution, that is a PL: , with s and κ being the energy and luminosity PL indices.As an alternative, we considered a broken power law (BPL): iso for E iso ≥ E b , where s 1 and s 2 are the low-and high-energy indices, respectively, and E b the break energy.We firstly tried to fix the low-energy index to zero, s 1 = 0, such that the resulting BPL is similar to a thresholded PL (Aschwanden 2015), although not mathematically equivalent to it.Similarly, we considered a BPL distribution for the luminosity distribution: f ′ (L iso ) ∝ L −κ1 iso for L iso < L b and f ′ (L iso ) ∝ L −κ2 iso for L iso ≥ L b , where κ 1 and κ 2 are the low-and high-luminosity indices, respectively, and L b the break luminosity.
For each redshift bin, we assigned a detection probability p = η z (E iso ) (p ′ = η ′ z (L iso )) to any given value of E iso (L iso ), that was randomly generated from f (E iso ) (f ′ (L iso )).Each simulated value of E iso (L iso ) was kept, provided that a Bernoulli trial 13 with the probability of success given by η z (E iso ) (η ′ z (L iso )) turned out to be 1.This process was iterated until we ended up with as many pulses as in the observed distribution for that red-shift bin.Consequently, the resulting simulated distribution has the same number of events and same redshift distribution as the observed distribution.To appreciate the impact of the selection effects, Figure 3 displays the simulated energy distribution for each of the 9 different redshift bins, along with the one predicted assuming dN/dE iso ∝ E −2 iso and the same number of GRBs per each redshift bin.

Validation of a putative intrinsic differential distribution
To assess the probability that the observed E iso (L iso ) distribution is the result of an assumed intrinsic differential distribution f (E iso ) ( f ′ (L iso )) net of selection effects, we carried out a few tests: likelihood ratio test (LRT), χ 2 , two-population KS, and Anderson-Darling 14 (AD) tests.
We divided the range of E iso (L iso ) values into M = 6 (6) logarithmic evenly spaced bins, ensuring that each bin contained ≥ 20 pulses.Let C o,i (C s,i ) be the number of pulses in bin i of the observed (simulated) distribution, either in the case of E iso or L iso .Let the total number of pulses in the observed and simulated distributions, respectively.For any simulated distribution s the χ 2 s and the LRT s are obtained by computing the following quantities: Both metrics were calculated for a sample of 100 simulated distributions.We then took the median values of both.In Equation ( 14), the numerator represents the maximum likelihood under the assumption that the two sets are independently distributed, whereas the denominator is the maximum likelihood under the alternative assumption that both sets share the parent distribution.Under the assumption that the two sets have a common distribution, LRT is χ 2 M −1 distributed: 15 we therefore carried out a one-tail test, calculating the p-value as P (χ 2 M −1 ≥ LRT).
14 The two-sample AD test was done using scipy.stats.andersonksamp 15 This is the reason why we incorporated a factor of 2 in Equation ( 14).
In the case of the E iso distribution, we were not able to find a satisfactory solution for the PL model according to all of the four tests: in particular, the best PL model yielded a value of ∼ 60 for both χ 2 and LRT (Equations ( 13) and ( 14)).Consequently, a simple PL model is confidently rejected.
Moving to the BPL model, given that two additional parameters come into play, we initially considered two alternative approaches: (1) we fixed s 2 = 1.5 and let both s 1 and E b free to vary; (2) we fixed s 1 = 0 and let both s 2 and E b free.This was done to make a preliminary exploration of the parameter space.Then we allowed all parameters to vary and performed Markov Chain Monte Carlo (MCMC) simulations with emcee (Foreman-Mackey et al. 2013) to sample the LRT around its minimum value.We took a Gaussian ball around the minimum LRT value of model ( 1) and performed a MCMC run with 1000 samples and 64 walkers.Then we compute the LRT on a cube of parameters around the minimum LRT value given by the MCMC run to find the minimum LRT value and the confidence intervals.We found that the BPL that minimized the χ 2 and the LRT was around s 1 ∼ 1, s 2 ∼ 1.67 and E b ∼ 1 × 10 52 erg, with χ 2 5 /5 = 5.25/5.Confidence intervals were found by considering regions of the parameter space that satisfy p α > 1 − α, where α is the confidence level, and p α is the p-value of the statistical test (either KS, AD, or χ 2 test).We took α = 0.95.Fit goodnesses (estimated with the p-values) and confidences intervals for the different models, either for the full sample or for the GS sample, are reported in Table 2. Analogous conclusions were obtained for the GS, that is a PL distribution cannot account for the data, whereas a BPL with similar shape as the one obtained for the whole sample does.Hence, the need for a break in the intrinsic distribution does not depend on the quality of the modeling of the GRB light curves, at least as long as dN/dE iso is concerned.Figure 4 shows the best fit model for dN/dE iso obtained for the full sample.
We obtained similar results for dN/dL iso for the full sample, i.e a simple PL distribution is unable to describe the observed distribution (the best fit is obtained around κ ∼ 1.3 with χ 2 5 /5 ∼ 20/5), while a BPL works better (the best fit is obtained around κ 1 ∼ 1.1, κ 2 ∼ 1.5, L b ∼ 6 × 10 52 erg s −1 with χ 2 5 /5 ∼ 3.4/5).For the GS sample, the observed distribution is compatible with a PL distribution, the BPL distribution leading only to a marginal improvement.Results regarding L iso distribution, either for the full sample or for the GS sample, can be found in Table 3. Figure 5 shows the BPL that best fits dN/dL iso of the full sample.4, except that here we are considering only pulses with redshift z < 1.76.
The need for a break in the E iso distribution could be either real or due to our possible inability to model very accurately the selection effects.To verify the latter possibility, we carried out the same analysis after excluding the three redshift bins with the highest redshift, which are most severely affected by the selection effects.As a result, we limited to pulses from GRBs with z max < 1.76, which make up ∼60% of the entire sample.Our conclusions remained essentially unchanged: the PL model is rejected (p-values from KS, AD, and χ 2 tests were 8 × 10 −3 , < 10 −3 , 10 −11 respectively), whereas the BPL model was acceptable, with best-fit values fully compatible with those obtained for the full sample.Figure 6 shows the BPL that best fits dN/dE iso of the sample of pulses with z < 1.76.
It is worth noting that we did not add pulses missed by mepsa in the simulated light curves, as we did in Section 2.4 in the case of the real ones.Consequently, the estimated detection efficiency in the low tail of the distribution is somewhat lower than the effective one acting on real data.This further supports the evidence for a break "a fortiori".

Comparison between pulse-rich and pulse-poor GRBs
We studied separately pulse-rich (N p > 6) and pulsepoor (N p ≤ 6) GRBs.The separation value was chosen so that the number of pulses belonging to pulse-poor GRBs is large enough to enable a statistical analysis.To this aim, we chose 6 pulses as the boundary, which corresponds to ∼ 20% of the whole sample, that is 220 pulses, belonging to the pulse-poor GRBs.In the case of dN/dE iso , a two-population KS (AD) test between the two subsamples yielded a p-value of 0.3 (> 0.25), thus providing no evidence for a different parent population.On the contrary, for dN/dL iso the same tests yielded a p-value of 10 −15 (< 0.001), thus rejecting a common parent distribution for the two groups.As shown in Figure 7, least luminous pulses in pulse-poor GRBs are relatively more abundant than in pulse-rich GRBs, although being comparably energetic, suggesting that they are on average less luminous and longer.As was noted in Section 3.2.1, the SS is predominantly populated by pulserich GRBs; thus, the difference in the luminosity distribution that was pointed out in the comparison between GS and SS is due to the pulse richness here considered.

Populations of pulses within individual GRBs
Thus far we studied the population of pulses as a whole, regardless of which GRB a given pulse may belong to.We explored more in detail this aspect, by testing whether pulses distribute among different GRBs completely randomly, starting from the overall sample.To this aim, we considered the distribution of peak luminosity of GRBs (that is, the luminosity of the most prominent pulse within each GRB) and focused on individual redshift bins, to limit the impact of selection effects.
For the 1.08 < z < 1.76 redshift bin, which includes most GRBs (N grb = 41, collecting N p = 287 pulses), we computed the GRB peak luminosity distribution.We then compared it with a distribution obtained by mixing together all the pulses belonging to the GRBs in the same redshift bin in the following way: for each GRB (i = 1, . . ., N grb ) we drew a random set of N p,i pulses, where N p,i is the real number of pulses of GRB i and took the peak luminosity of simulated GRB i. Repeating this for all of the N grb GRBs, we obtained a fake peak luminosity distribution, which appears to be significantly different from the real one (AD test p-value < 10 −3 that both sets share a common parent population): in particular, the fake distribution has a more populated hard tail (see Figure 8).This indicates that pulses belonging to a given GRB are not completely independent of one another, or, at least, that they are not the result of randomly assembling from the whole population of observed pulses.High luminosity pulses, in particular, tend to cluster within a few GRBs, rather than being randomly distributed among all the GRBs.
To further test this possibility, we performed a multinomial test.Firstly, we counted N grb,lum = 4 GRBs that contain the top 10% most luminous pulses (n = 28 out of 287 pulses).We then compared this number with the one that would be obtained by randomly distributing the pulses among all the GRBs, keeping the same distribution of number of pulses per each GRB.This problem is similar to simulate the outcomes of n throws of a dice with m faces and can be addressed by sampling from the multinomial distribution, described by its probability mass function:16 where n = (n 1 , n 2 , . . ., n m ) is the array of number of luminous pulses in each GRB ( m i=1 n i = n = 28), m = N grb = 41, and p = (p 1 , p 2 , . . ., p m ), with p i = N p,i /N p being the probability for a pulse to belong to GRB i, calculated from all N p = 287 pulses in the chosen redshift bin.We performed 10 7 simulations and in all cases the number of GRBs that contained the luminous pulses was higher (four times in average) than the real one, N grb,lum = 4.The p-value that luminous pulses randomly distribute among GRBs within this redshift bin is < 5 × 10 −5 .
We replicated the same analysis on the next redshift bin, 1.76 < z < 2.90, which contains N grb = 25 and N p = 253 pulses.The top 10% luminous pulses (n = 25) belong to 5 GRBs.The same test yielded a p-value of 0.9 × 10 −4 that luminous pulses are distributed over ≤ 5 GRBs.We can therefore conclude that the most luminous pulses are more likely to belong to fewer GRBs than what is expected from a pure random distribution.

Waiting time and duration distributions
We obtained the distribution of the rest-frame WT calculated as ∆t i = (t 0,i+1 −t 0,i )/(1+z) where t 0,i is the pulse peak time of the i th pulse.We also calculated the pulse duration distribution, where duration was defined as the full width at half maximum (FWHM) of the modeled pulse, T i = τ ri (1+ξ i ) log(2) 1/νi .Durations were not corrected for cosmological dilation because of the combination of different effects that mostly cancel each other (see Camisasca et al. 2023 and references therein for a detailed explanation).To ease the comparison with previous results, we modeled the distributions as in Guidorzi et al. (2015), using the following model: consisting of a PL hard-tail with a PL index 3 − γ, plus a characteristic break time δ below which the distribution flattens (γ and δ corresponds to α and β in Equation 8 of Guidorzi et al. 2015).The first bin of the two distributions was not considered (∆t ≥ 0.05 s) since mepsa efficiency drops at low WTs/durations.We sampled the posterior distribution by running a MCMC simulation with 10 4 steps and 32 walkers.We discarded the 10 3 first steps of the posterior distribution and we used only every 15 steps from the chain.For the duration and the WTs distribution, we obtained a PL index 3−γ T = 2.08 +0.19  −0.16 and 3−γ W T = 2.04 +0.14 −0.12 , respectively.Best fit parameters are reported in Table 4 while the distributions along with the best fit models are shown in Figures 9 and 10.

Effects of the pulse model
We repeated the analysis by adopting the alternative pulse model of Norris et al. (2005), to explore to which extent our results are sensitive to the choice of the model   4, except that here it is the differential duration distribution dN/d∆t.The black area represents the first bin which is not considered in the computation of the posterior distribution.described in Section 2.4.Overall, the fit quality is worse than the one obtained with the former pulse model, since the number of GRB LCs that do not fulfill the condition |σ ϵ − 1| > 0.1 (see Section 3) is 57 vs. former 15.This was to be expected because of the smaller number of model parameters.We then computed the four different distributions considered in this paper and compared them with the ones obtained with the former pulse model (Figure 11).
We also performed two-population KS and AD tests.All but one test yielded p-values above 5%, except for 0.3% obtained for the time duration distribution.This is due to the fact that the FWHM computed with the second model is on average longer (by a factor ≲ 2) than the one computed with the first model.The distributions were modeled following the same procedure adopted for the first pulse model.Overall, the powerlaw indices do not differ significantly from the previous corresponding ones: α (2) W T = 2.14 +0.16 −0.14 .Compared with previous values reported in Table 5, they are all consistent.While the need for a break in energy/luminosity distributions is confirmed for both pulse models, in the luminosity case the value of the break point depends on the pulse model: for the second one, we found L (2) break ∼ 1 × 10 52 erg s −1 , which deviates by 2.5 σ from previous L (1) break ∼ 6 × 10 52 erg s −1 .As a result, we proved that our results do not depend significantly on the pulse shape model, except for the value of the break luminosity, which we conservatively estimate in the range L b ∼ 10 52−53 erg s −1 .

DISCUSSION
For the distributions of energy, peak flux (or the corresponding intrinsic quantity which is luminosity), and duration, the corresponding PL indices predicted by SOC models are given by the following (Aschwanden 2014): where β is the spreading exponent (β = 1/2 for subdiffusion, β = 1 for normal diffusion, β = 2 for linear expansion), d is the Euclidean space dimension, and D d is the fractal dimension of SOC avalanches, which is usually approximated by D d ≃ 1+d 2 .SOC theory also predicts α WT = α T , where α WT is the PL index of the WT distribution (this holds for short WTs: for WTs longer than the longest pulse duration, the WT distribution should break exponentially; see Aschwanden 2014).In the case of 3D avalanches (d = 3) and normal diffusion (β = 1), it is α E = 1.5, α P = 5/3, and α T = α WT = 2. Our estimates of the duration and of the WT duration distribution PL indices agree with SOC predictions: α T = 3 − γ T = 2.08 +0.19  −0.16 and α WT = 3 − γ W T = 2.04 +0.14  −0.12 .The PL indices of the high-value tails of the energy and luminosity distributions, respectively s 2 = 1.67 +0.23   −0.16   and κ 2 = 1.47 +0.60  −0.18 , are also roughly in agreement with SOC predictions (see Figure 12).However, our results show evidence for a break in both distributions, which is not predicted by SOC theory.Given the accuracy with which the detection efficiency was modeled, the evidence for a break seems to be hardly entirely ascribable to unaccounted selection effects.In our study, we have selected all Type-II GRB candidates.Our study implicitly assumes that all Type-II GRB engines are behaving in the same way.However, one cannot reject the possibility that the observed populations is actually the mixture of different kinds of engines or at least different behaviors.Even for a given type of central engine, say an hyper accreting BH, two or more jet formation scenarios (ν ν annihilation or BZ effect) producing different energy/luminosity distributions, can be at play and contribute to the observed population of GRBs.In principle, the observed deviations from the SOC predictions could be also ascribed to the physics and geometry of the prompt emission, i.e., the properties of the jet and its interactions with the stellar envelope.Finally, a PL behavior does not unavoidably imply SOC, is just a necessary condition, not a sufficient one, given that hard-tailed distributions modeled as PL or BPL can be the result from a number of different processes.
Table 5 summarizes and compares our findings with analogous studies.Interestingly, the energy index obtained is consistent with what was found for the burst activity observed in Galactic magnetars: α E = 1.43 − 1.76 (Göǧüş et al. 1999;Göǧüş et al. 2000), and in solar flares: α E = 1.62 +0.12 −0.12 (Aschwanden 2011).We obtained substantially lower values for α P than Lyu et al. (2020) and Li & Yang (2023).In these studies, α P is determined by analyzing the distribution of pulse count rates given by the instrument, as it was also done for solar flares.While the latter case is justified, as all flares come from the same source at a fixed distance from the detector, for GRBs the range of distances is so large that one should either consider a group of GRBs with similar redshift or better use luminosities and model the selection effects that inevitably affect the observed sample.
We compared dN/dL iso with the GRB (peak) luminosity function, as modeled for example in Ghirlanda & Salvaterra (2022).Even though the latter is in principle different from dN/dL iso , the two distributions are not completely independent of one another.Ghirlanda & Salvaterra (2022) modeled the distribution with a BPL with low/high luminosity indices a 1 ∼ 1 and a 2 ∼ 2.2 and a break luminosity L b ∼ 10 52 (1 + z) 0.64 erg s −1 slightly dependent on redshift.Our values are consistent with the low-luminosity index and with the break luminosity of the luminosity function.We find a flatter highluminosity index that can be understood since our distribution contains all pulses of all GRBs, while the GRB peak luminosity function is determined by the most luminous pulse within each GRB.As we already pointed out, high-luminosity pulses tend to cluster in relatively few luminous GRBs: consequently, the distribution resulting from the selection of the most luminous pulse of each GRB turns into a depletion of luminous pulses and hence into a smaller fraction of high-luminosity events in the GRB distribution compared with the pulse distribution.
Finally, our results do not seem to crucially depend on redshift, given that they did not change in essence when we ignored GRBs with redshift z > 1.76.Specifically, the evidence for a break in both distributions holds true in both cases, suggesting that the break is not likely to be entirely due to selection effects, but it is an intrinsic feature.

CONCLUSIONS
We determined the energy, luminosity, duration, and waiting time distributions of individual pulses identified with a well-calibrated algorithm from 142 Type-II GRBs with known redshift detected by the GBM.We then carried out a careful analysis of the selection effects through a suite of simulations, dividing our sample into 9 redshift bins, and modeled for each of them the detection efficiency as a function of energy (luminosity), η z (E iso ) (η ′ z (L iso )).For each redshift bin, we have then generated energy/luminosity samples of pulses accounting for the detection efficiency, under the assumption of some putative distribution models like simple PL or BPL.The plausibility of each assumed distribution was finally tested by comparing the resulting distribution with the observed one.
We found that a simple PL can reproduce neither E iso nor L iso distributions (especially in the case of E iso ).Rather, for the first time we found evidence for a break in E iso and L iso distributions (E b ∼ 10 52 erg and L b ∼ 10 52−53 erg s −1 ), which appears to be hardly ascribable to unaccounted selection effects.A possible interpretation for this break is that the underlying assumption of a unique stochastic process ruling the GRB dynamics is not valid, but instead different kinds of dynamics, possibly connected with different progenitors or different regimes, contribute to the observed population of GRBs.
Interpreting our results in the SOC framework with β = 1, our results are compatible with prompt emission being the result of 3D (d = 3) avalanches produced by a non-linear dynamical system driven slowly to a SOC state.The SOC interpration of GRB prompt emission could be compatible with utterly different theoretical models.Indeed, this scenario could be consistent either with the picture of a NDAF surrounding a hyper-accreting stellar mass BH becoming thermally unstable and cooled by neutrino emission, with rotational energy being extracted by BZ effect and converted into a Poynting-flux dominated outflow, or with cascades of magnetic reconnection events produced within a magnetically dominated relativistic outflow at the GRB emission site, as foreseen in the ICMART model (Zhang & Yan 2011).
We acknowledge the anonymous reviewer for a constructive report which improved the manuscript.R. We computed the isotropic-equivalent released energy using Equation ( 9), thus assuming a negligible impact of a possible spectral evolution.
GRB spectra are known to be temporally evolving.This is usually described in terms of two alternative behaviors: (i) a monotonic hard-to-soft evolution and (ii) pulse tracking (see Lu et al. 2012 and references therein).Therefore, the isotropic-equivalent energy should be computed using the time-resolved fluence Φ ′ k instead of  the time-average Φ k used in Equation ( 9), where Φ ′ k is calculated from modeling the spectrum extracted within the time interval that includes only the k-th pulse (or, at least, centered on it, given that pulses occasionally overlap and therefore their spectra cannot be completely separated).
We conveniently defined the time-resolved conversion factor f k = Φ ′ k /C ′ k , where C ′ k are the counts integrated over the same time interval used to compute Φ ′ k .In this way, the degree of approximation introduced in Equation ( 9) by assuming time-average instead of time-resolved fluences can be studied through the ratio f k /f tot .
To this aim, we considered a few cases of bright GRBs with numerous pulses, some of them being already known to show spectral evolution.
We modeled time-resolved spectra of these GRBs and we studied the impact of spectral evolution on the E iso distribution.For this study, we considered 5 GRBs among our sample17 having respectively 57, 34, 21, 20 and 13 pulses thus totaling 145, that is enough to build a statistically sound distribution.These GRBs are among the brightest of the GBM catalog, with three of them ranking in the 12 top fluence of the whole catalog. 18 Figure 13 illustrates the case of 180720B: for the timeaveraged spectrum, we took the fluence from the GRB catalog Φ tot = (2.9853± 0.0008) × 10 −4 erg cm −2 , integrated from −60.416 to 137.220 s.Our results are consistent with those by Chen et al. (2021).Our analysis shows a hard-to-soft evolution, modulated by f /f tot tracking the light curve peaks (Figure 13).As a consequence, the fluence of the brightest peaks is underestimated by up to 40%, whereas the fluence of the weakest ones is overestimated by a comparable amount.Among the examined GRBs, 180720B, 190114C, and 171010A  exhibited the larger deviations of f /f tot from 1, reaching 60-70%, indicating that the effect is probably stronger for brightest GRBs.We can therefore safely assume that neglecting the spectral evolution for the bulk of GRBs considered in this work has a milder effect than what we assessed for the brightest cases.We finally explored to what extent the E iso distribution was affected by ignoring the spectral evolution: we compared E iso calculated with both time-resolved and time-average analysis and derived the corresponding distributions (Figure 14).We ran statistical tests to assess the null hypothesis that the two E iso distributions are drawn from the same population.A KS test yielded a p-value of 0.87, while the AD test gave a lower-limit on the p-value > 0.25.We concluded that we cannot reject the hypothesis that the two distributions are drawn from the parent population, so we can safely ignore the bias introduced by neglecting the spectral evolution.

Figure 1 .
Figure 1.150314A background-subtracted light curve (red) along with the best-fit model (black).The typical error size is shown in the top left.The black horizontal line represents the interpolated background.Blue dashed curves represent the individual pulses detected by mepsa while orange ones represent pulses added by visual inspection to improve the fit of this complex GRB time profile.

Figure 9 .
Figure9.Same as Figure4, except that here it is the differential duration distribution dN/dT .The black area represents the first bin which is not considered in the computation of the posterior distribution.

Figure 10 .
Figure10.Same as Figure4, except that here it is the differential duration distribution dN/d∆t.The black area represents the first bin which is not considered in the computation of the posterior distribution.
M. and C.G. acknowledge the Dept. of Physics and Earth Science of the University of Ferrara for the financial support through the "FIRD 2022" grant.R.M. acknowledges the University of Ferrara for the financial support of his PhD scholarship.A.T. acknowledges financial support from "ASI-INAF Accordo Attuativo HERMES Pathinder operazioni n. 2022-25-HH.0".
Energy, luminosity, duration, and waiting times PL indices obtained in our study and in the literature in the case of GRB prompt emission, X-ray flares, precursors, solar flares, and SGRs.Values predicted by SOC theory are also indicated, in the case where β = 1.Special cases of d = 1 or d = 3 are also reported.The last two columns indicate whether the obtained values are in agreement with the values predicted by SOC theory in the two special cases.

Figure 12 .Figure 13 .
Figure 12.Top panels, left to right: violin plots of the posterior distributions for the PL indices of the energy, luminosity, duration, and waiting times distributions, respectively, obtained in this work on GRB prompt emission (for the energy and luminosity distributions, the post-break values were considered); the horizontal bars span 5 to 95% quantiles.Dashed lines show the SOC predictions for d = 3. Bottom panels, left to right: the same corresponding violin plots as in the top, including values reported from the literature.For GRB prompt emission: Lyu et al. 2020 (L20; brown), Li & Yang 2023 (LY23; cyan), and Guidorzi et al. (2015) (G15(GBM) and G15(BAT); in pink and purple).For GRB X-ray flares, values from Wang & Dai 2013 (WD13; blue), Wei 2023 (W23; green).The orange point refers to a study where prompt emission pulses and X-ray flares were considered as a unique process,(Guidorzi et al. 2015; G15(BAT-X)).The shaded areas show the results reported in the literature from analogous investigations about different classes of astrophysical sources: for solar flares, values fromWang &  Dai 2013 (WD13; black)  and Aschwanden 2011 (A11; red).For magnetars, values reported inGöǧüş et al. 1999; Göǧüş et al.  2000 are displayed (G99; yellow).Dashed (dotted) lines show SOC predictions for d = 3 (d = 1) and β = 1.

Table 1 .
The first 5 GRBs of the final sample.This table is available in its entirety in machine-readable form.

Table 2 .
Results of the modeling of the dN/dEiso distribution, once selection effects are accounted for.Values of frozen parameters are reported in square brackets.

Table 3 .
Results of the modeling of the dN/dLiso distribution, once selection effects are accounted for.Simulated distributions of Eiso, depending on whether selection effects due to detection efficiency are considered (orange) or not (blue).The nine panels refer to increasing redshift bins (left to right, top to bottom).The blue distribution was drawn assuming dN/dEiso ∝ E −2 iso .Vertical dashed lines mark the 50% detection efficiency: pulses with lower values of Eiso, highlighted by the gray area, are therefore hampered by a low detection efficiency.For the sake of clarity, 10 5 pulses were simulated in each redshift bin.The black solid line represents the best-fit BPL intrinsic differential energy distribution dN/dEiso obtained for the full sample.The gray shaded histogram is a binned version of the same distribution.The orange distribution was simulated from the intrinsic one after accounting for the selection effects, to be compared with the red distribution, which is the observed one.
Figure 6.Same as Figure

Table 4 .
Results of the modeling of the differential distributions dN/d∆t and dN/dT .Uncertainties are given with a 90 % confidence level.