Bayesian Time-resolved Spectroscopy of Multipulse GRBs: Variations of Emission Properties among Pulses

We use data obtained by the Fermi Gamma-ray Space Telescope, which was launched in 2008, and carries two instruments: the GBM (Meegan et al. 2009) and the Large Area Telescope (LAT; Atwood et al. 2009). Together, they cover an energy range from a few keV to a few hundred GeV. The GBM harbors 14 detectors, of which 12 are sodium iodide (NaI; 8 keV–1 MeV) and two are bismuth germanate (BGO; 200 keV–40 MeV) scintillators. By 2019 June, Fermi had completed 11 yr of operation, and at least 2388 GRBs had been observed.

Among these bursts, we want to identify the ones that have at least two individual emission episodes, which can be assumed to be independent from each other. Since there is no theoretical prediction of the shape and form of individual emission episodes, we choose to make a general definition for our selection, which does not limit the sample to smooth pulses only. Indeed, numerical simulations of jet emission propagating through the progenitor star show that not only smooth pulses are expected but also more complex morphologies (Lazzati et al. 2009; López-Cámara et al. 2014). This choice also avoids omitting too many bursts. Following Yu et al. (2019), we thus allow subdominant variations on top of the main pulse structures. In addition, we also include emission activities with many smaller spikes but whose heights are limited by an approximate pulse-shaped envelope. The pulse-shaped envelope indicates that such emission could be connected, in spite of its large variability. Examples of the light curves that are selected for the sample are shown in Section 3.2 and the Appendix.

We first visually inspected the time-tagged event (TTE) light curves from each of the 2388 GRBs obtained from NASA/HEASARC. ⁹ We identified more than 120 bursts that had at least two clear emission episodes. The light curves exhibit diverse temporal properties. Some have a precursor-like pulse followed by a main pulse, while others exhibit a main pulse followed by a small pulse, and yet others consist of several pulses with similar strength. Many of the pulses are well separated by quiescent intervals, while others have a slight temporal overlap.

We use the standard selection criteria adopted in GBM catalogs (Goldstein et al. 2012; Gruber et al. 2014; Yu et al. 2016). Following the common practice, we choose the triggered detectors that have a viewing angle of less than 60° (Goldstein et al. 2012). Typically, one to three NaI and one BGO detector are selected. The period of GRB emission that is considered is in most cases somewhat longer than the T₉₀ reported in the HEASARC database. This is done so that all relevant features in the light curve can be incorporated. In order to determine the background emission, we select one interval located tens of seconds before the triggered time and one interval located tens of seconds after the emission has ended. We fit the background photon count level with a polynomial function, which typically has an order lower than 4. The optimal order is determined through a likelihood ratio test. The polynomial is applied to fit all 128 energy channels and then interpolated into the pulse interval to yield the background photon count estimate. In a few bursts, there are several pulses that are separated by long quiescent intervals. In these cases, we select three background intervals to better constrain the background levels (e.g., GRB 140810782). The spectral energy range is set from 10 to 900 keV for the NaI detectors and 300 keV to 30 MeV for the BGO detectors. In order to avoid the K-edge at 33.17 keV, we ignore the range 30–40 keV.

In order to follow and study spectral evolution in individual bursts, detailed time-resolved spectroscopy is needed (see, e.g., Crider et al. 1997; Kaneko et al. 2006; Goldstein et al. 2012). The question is then how to divide the light curve into time bins, since this choice might affect the results. Foremost, it is important to minimize the amount of variation of the emission during a time bin, since such variations will obscure the intrinsic spectral shape. A method to account for this is to identify Bayesian blocks (BBlocks) in the light curve (Scargle et al. 2013). Therefore, we apply the BBlock method with a false-alarm probability p₀ = 0.01 to the TTE light curve of the brightest NaI detector (see Li 2019a; Yu et al. 2019, for further details). This binning is then used for the other detectors as well.

The BBlock method will create time bins that have varied signal-to-noise ratios. Therefore, some time bins will not have enough signal for a fit to be reliable. A suitable measure of the signal-to-noise ratio is the statistical significance S. ¹⁰ Dereli-Bégué et al. (2020) demonstrated that to fully determine the spectral shape, a value of S = 15–20 is needed. In particular, this ensures that the model parameters converge properly (see also, e.g., Vianello et al. 2018; Li 2019a, 2019b; Ryde et al. 2019). We therefore follow this recommendation and only select the BBlock time bins that have S ≥ 20.

This specific selection thus provides time bins during which (i) there is no significant spectral evolution and (ii) there is highly significant data. Such spectra are Required in order to make firm inferences on the intrinsic emission spectrum. However, a consequence of this selection is that intervals with rapid variations will be dismissed. This is often the case for the rise phase of pulses, which could carry particular information about the emission process (see further discussion in Section 5.2). A possibility of incorporating these dismissed intervals would be to perform joint fits of many intervals at the same time by assuming a prescription of the spectral evolution. This would increase the statistical significance while maintaining the temporal resolution. However, since further assumptions need to be made on the spectral evolution (e.g., Ryde & Svensson 1999), such studies are deferred to other publications.

The purpose of this study is to relate the emission properties between pulses within a burst. We therefore only consider bursts that have multiple distinguishable pulses. One of the properties that we will relate is the spectral evolution during the individual pulses. Hence, in order to include a pulse in the study, it needs to have at least four time bins with S ≥ 20 (N_(S≥20) ≥ 4). This allows one to determine the spectral evolution and the correlation between spectral parameters. Our final selection criterion is, therefore, that the burst should have at least two such pulses. This criterion strongly reduces the sample size. However, the pulses that pass the criterion will provide well-determined spectra and evolution characteristics. This final selection yields a sample of 39 bursts. From these, we obtain 117 pulses and 1228 time-resolved spectra, of which 103 pulses have N_(S≥20) ≥ 4, consisting of 944 spectra.

The properties of our sample are summarized in Table 1 and include the Fermi/GBM ID (column 1), redshift (column 2), and observed duration t₉₀ (column 3), together with the used detectors (column 4), the selected source and background intervals (columns 5 and 6), the number of total (column 7) and used (S ≥ 20; parentheses in column 7) BBlock time bins, and the number of total (column 8) and used (parentheses in column 8) pulses for each burst. The detector in parentheses is the brightest one, used for the BBlocks and background determinations.

Inevitably, there will be some overlap between many of the pulses. During such periods, emission from both pulses will be present, and the observed spectrum will be a superposition of two spectra. Depending on their mutual strengths, a fit with a single spectral component might give misleading results. We therefore group the pulses in the 39 bursts into three subsamples.

• 1.  
  Gold. No overlap, and the pulses are completely independent. There are 13 such pulses, which are listed in Table 2 (column 5).
• 2.  
  Silver. Slight overlap. There are 90 such pulses (column 5 in Table 2).
• 3.  
  Bronze. Pulses that have less than four time bins S ≥ 20. There are 14 such pulses (column 5 in Table 2). Note that these pulses are not used in our final statistical analysis.

The physical model of the GRB spectra is not yet firmly established. Likewise, the possible regimes of the suggested emission models are not clarified. Therefore, empirical models are typically used in the analysis. This offers the flexibility to explore many different models and, at the same time, to explore the models without limiting certain spectral regimes. The typically used functions are the Band function, which is a broken power-law function, and the cutoff power-law (CPL) function (Band et al. 1993; Gruber et al. 2014). Another advantage of these empirical models is that they are used in all catalogs mapping the behavior and characteristics of observed bursts.

However, the empirical models do not necessarily correspond to the underlying physical model (e.g., Lloyd-Ronning & Petrosian 2002; Burgess 2014; Acuner et al. 2019). Therefore, the empirical parameter values that are inferred from the data cannot be directly translated into physical model parameters. However, the way the physical model parameters map onto the Band parameters has been established for a few emission models (e.g., Lloyd & Petrosian 2000; Acuner et al. 2019; see also Section 4). Moreover, Acuner et al. (2020) showed that the values of, e.g., α can be used to decisively distinguish between various types of emission models, as long as the data have a high signal strength (as is the case in the present study).

We will therefore make use of these empirical functions. The Band function is defined by the low-energy power-law index (α) and the high-energy index (β), which are connected smoothly around the break energy (E₀).

The photon number spectrum is defined as

$$\begin{eqnarray}\begin{array}{r}{f}_{{\rm{B}}{\rm{a}}{\rm{n}}{\rm{d}}}({E})=A\left\{\begin{array}{ll}{\left(\displaystyle \frac{{E}}{{{E}}_{{\rm{p}}{\rm{i}}{\rm{v}}}}\right)}^{\alpha }\exp \left(-\displaystyle \frac{{E}}{{{E}}_{0}}\right), & {E}\leqslant (\alpha -\beta ){{E}}_{0}\\ {\left[\displaystyle \frac{(\alpha -\beta ){{E}}_{0}}{{{E}}_{{\rm{p}}{\rm{i}}{\rm{v}}}}\right]}^{(\alpha -\beta )}\exp (\beta -\alpha ){\left(\displaystyle \frac{{E}}{{{E}}_{{\rm{p}}{\rm{i}}{\rm{v}}}}\right)}^{\beta }, & {E}\geqslant (\alpha -\beta ){{E}}_{0}\end{array}\right..\end{array}\end{eqnarray} \tag{ 1 }$$

The energy of the spectral peak of the ν F_ν spectrum (assuming β < −2) is in units of keV,

$$\begin{eqnarray}&&{E}_{{\rm{p}}}=(2+\alpha ){E}_{0},\end{eqnarray} \tag{ 2 }$$

A is the normalization factor at 100 keV in units of ph cm⁻² keV⁻¹ s⁻¹, E_piv is the pivot energy fixed at 100 keV, and α and β are the low- and high-energy power-law photon spectral indices, respectively.

The CPL function is given by

$$\begin{eqnarray}&&{f}_{\mathrm{CPL}}({E})=A{\left(\displaystyle \frac{{E}}{{{E}}_{\mathrm{piv}}}\right)}^{\alpha }\exp \left(-\displaystyle \frac{{E}}{{{E}}_{c}}\right).\end{eqnarray} \tag{ 3 }$$

Note that the CPL model approaches the Band model as β tends to −∞ . The peak energy E_p of the ν F_ν spectrum is related to E_c through E_p = (2+α)E_c.

The analysis in this work is performed with the Multi-Mission Maximum Likelihood Framework (3ML; Vianello et al. 2015), and we use a Bayesian approach, using the Markov Chain Monte Carlo (MCMC) method (e.g., Foreman-Mackey et al. 2013). We use the typical spectral parameters obtained from the previous Fermi/GBM catalog as the prior information of the Bayesian inference.

For the Band model,

$$\begin{eqnarray}\left\{\begin{array}{ll}A\sim \mathrm{log}{ \mathcal N }(\mu =0,\sigma =2) & {\mathrm{cm}}^{-2}\ {\mathrm{keV}}^{-1}\ {{\rm{s}}}^{-1}\\ \alpha \sim { \mathcal N }(\mu =-1,\sigma =0.5) & \\ \beta \sim { \mathcal N }(\mu =-2,\sigma =0.5) & \\ {E}_{{\rm{p}}}\sim \mathrm{log}{ \mathcal N }(\mu =2,\sigma =1) & \mathrm{keV}\end{array}\right.,\end{eqnarray} \tag{ 4 }$$

and for CPL model,

$$\begin{eqnarray}\left\{\begin{array}{ll}A\sim \mathrm{log}{ \mathcal N }(\mu =0,\sigma =2) & {\mathrm{cm}}^{-2}\ {\mathrm{keV}}^{-1}\ {{\rm{s}}}^{-1}\\ \alpha \sim { \mathcal N }(\mu =-1,\sigma =0.5) & \\ {E}_{{\rm{c}}}\sim \mathrm{log}{ \mathcal N }(\mu =2,\sigma =1) & \mathrm{keV}.\end{array}\right.\end{eqnarray} \tag{ 5 }$$

A posterior distribution is obtained from the prior distribution and the likelihood that combines the model and the observed data. The best model parameters are estimated from the posterior probability distribution obtained by MCMC sampling. When using MCMC sampling, in order to obtain the steady-state chains of the parameter distributions, the sampling needs to reach a certain number of times; therefore, the first part of the sample that has not reached the steady-state distribution is discarded. Therefore, for each parameter estimation, we take 10,000 MCMC samples and discard the initial 20%. The fitted parameters are estimated by the maximum a posteriori probability (MAP) with uncertainties at the 1σ (68%) Bayesian credible level and based the last 80% of the MCMC samples. We also provide all of the analysis results of each time-resolved spectrum, which includes the best parameter value estimates, covariance matrices, and the statistical information criteria. They can be retrieved at doi:10.5281/zenodo.4746267.

For a majority of time-resolved burst spectra, the CPL model is a sufficient model (e.g., Kaneko et al. 2006; Goldstein et al. 2012; Burgess et al. 2019; Yu et al. 2019). However, in some cases, a high-energy power law significantly improves the fit, indicating a significant flux contribution beyond the spectral peak. In the case of synchrotron emission, the high-energy power law provides information about the particle acceleration and the nature of the shocks, while in the case of photospheric emission, it provides information about the energy dissipation in the emitting region.

In order to determine whether a high-energy power law gives significant improvement in our analysis, we compare the information criteria of the Band and CPL fits. Acuner et al. (2020) showed that the information criteria capture everything important in the fits and that the difference in information criteria can be used in the model comparison. In particular, they showed that a significance of 99% of preferring one model over the other is found for a difference in log evidence greater than 5 (their Equation (6)) and that this corresponds approximately to a difference in information criteria of 10 (further discussion can be found in Acuner 2019). In this work, we therefore adopt the information criterion to compare models, and, following the literature in the field (e.g., Greiner et al. 2016; Yu et al. 2019), we particularly use the deviance information criterion (DIC; Spiegelhalter et al. 2002; Moreno et al. 2013), which is defined as DIC = −2log[p(data $| \hat{\theta }$)] + 2p_DIC, where $\hat{\theta }$ is the posterior mean of the parameters, and p_DIC is a term to penalize the more complex model for overfitting (Gelman et al. 2014). We consequently accept the Band function as the preferred model if the difference between the Band's DIC and the CPL's DIC ΔDIC = DIC_Band − DIC_CPL > −10.

For consistency, one needs to use the same empirical model throughout the whole pulse in order to avoid artificial fluctuation due to a change of spectral model (see Yu et al. 2019). Consequently, if one time bin has ΔDIC < −10, then we use the Band function throughout the pulse. Otherwise, we use the simpler CPL model. Therefore, the pulses that are fitted by a Band function have at least one time bin in which a high-energy power law is significantly detected. The model used for each pulse is listed in column 6 in Table 2.

Among all of the individual time bins, we find that in only 29% (274/944) is ΔDIC < −10; i.e., the Band function is preferred. For these time bins, a high-energy power law is required by the data (see also Gruber et al. 2014; Yu et al. 2019). For the sequence of pulses, we find the corresponding fractions to be 36% (122/338) for the first pulse P₁, 28% (101/361) for P₂, 22% (34/156) for P₃, 15% (7/46) for P₄, 15% (4/26) for P₅, and 29% (5/17) for P₆. There is only a weak trend, with the first pulse P₁ having the largest fraction.

Turning over to the pulses, we find that 66% (77/117) of them have at least one time bin with ΔDIC < −10; therefore, the Band function is used for the entire pulse. For the sequence of pulses, we find that the corresponding fractions are 77% (30/39) for P₁, 69% (27/39) for P₂, 63% (15/24) for P₃, 22% (2/9) for P₄, 50% (2/4) for P₅, and 50% (1/2) for P₆. Again, there is only a weak trend, with the first pulse P₁ having a larger fraction.

In addition to the ΔDIC criterion, we also examined the p_DIC values following Yu et al. (2019). In some instances, the values were found to be anomalously large (hundreds of thousands), which typically indicates that a local rather than a global minimum of the likelihood function is found. In these cases, we reran the Bayesian fits with new initial values to ensure proper convergence to the global minimum.

Figure 2 shows the parameter distributions for every group of pulses, including α, E_p, and the K-corrected energy flux F (erg cm⁻² s⁻¹), over the range 1–10⁴ keV, as well as the duration, Δt, of the pulses. In the following figures, the data points are orange (P₁), blue–magenta (P₂), violet (P₃), yellow (P₄ or P₄₊₅₊₆), ¹¹ cyan (P₅), and green (P₆). The best Gaussian fit for each distribution is presented, and the corresponding average values and standard deviations are presented in Table 3. The average value of the full sample, including all pulses, is α = −0.84 ± 0.35 and E_p = log₁₀(214) ± 0.42. These values are in agreement with previous catalogs (e.g., Kaneko et al. 2006).

Zoom In Zoom Out Reset image size

In order to assess whether the distributions change between the different pulse groups, we use the Kolmogorov–Smirnov (K-S) test. This test determines the chance probability, P, that two distributions are sampled from populations with identical distributions. For our purposes, we use P < 10⁻² to ensure that the distributions are truly different. In Table 4, we present the P-values for all of the α and E_p distributions for the pulse groups.

These results indicate that the groups of pulses could be divided into two categories with two different behaviors. The first category, with P₁ and P₂, is different from the rest of the pulse groups. In the first group of two pulses (P₁ and P₂), the peak energy E_p does not change, while the spectral slope α exhibits a clear softening. Compared to the second category, with the later pulse groups (P₃–P₆), there is a change, with both α and E_p decreasing. However, within this second category, the distributions are not significantly different. Based on this analysis alone, it can, therefore, be argued that these two categories represent different types of spectral characteristics and, therefore, possibly different emissions.

The comparison between the fluxes and pulse durations is shown in the bottom panels of Figure 2. The fluxes have a similar pattern as α. There is a steady decrease, except for the last two groups, which have a similar distribution. Finally, the pulse duration distributions are comparable, and the first three groups have increasing average values.

Another way to compare the variation of emission properties in between pulses is to study the maximal value of the parameters of α and E_p in each pulse and how they vary. Note that ${\alpha }_{\max }$ and ${E}_{{\rm{p}},\max }$ do not necessarily need to be at the same time bin. In the upper panels in Figure 3, we present ${\alpha }_{\max }={\alpha }_{\max }(t)$ and ${E}_{{\rm{p}},\max }={E}_{{\rm{p}},\max }(t)$ versus time. Here the trend of softening is again revealed for ${\alpha }_{\max }$. Later pulses typically have smaller ${\alpha }_{\max }$. There is also a trend that ${\alpha }_{\max }$ for the first pulse is softer the further it is delayed from the trigger. For ${E}_{{\rm{p}},\max }$, however, there is no trend, and the distinction between the two categories from the K-S tests is not as apparent. This means that, when it comes to ${E}_{{\rm{p}},\max }$ and ${\alpha }_{\max }$, it is mainly the spectral shape that changes and not the location of the peak energy when we compare the six pulse groups.

Zoom In Zoom Out Reset image size

In the lower panels in Figure 3, the dependency of ${\alpha }_{\max }$ and ${E}_{{\rm{p}},\max }$ on pulse duration Δt is shown. The colored stars correspond to the average values. For the first three groups, the increase in average pulse duration is correlated with the decrease in the average values of ${\alpha }_{\max }$ and ${E}_{{\rm{p}},\max }$ (see Figure 2). Again, the last group with P_(4+5+6) differs from the trend by having a shorter (average) pulse duration. However, instrumental effects might play a role here. The fluxes of later pulses are lower; hence, the full duration of the pulses might not be apparent above the background noise level.

Early investigations of GRB emission identified a significant correlation between spectral properties, such as a relation between the intensity and the shape of the spectrum (Wheaton et al. 1973), and a correlation between the intensity and the spectral peak (Golenetskii et al. 1983). In this section, we classify the evolution of α and E_p relative to the count light curves, as commonly done in the literature (e.g., Kargatis et al. 1994; Ford et al. 1995). In Section 3.3, we will quantify this further by investigating the actual correlations between the parameters.

Examples of the evolution of the spectral parameters E_p and α and the ν F_ν flux are provided in the upper panels of Figure 4. Here the parameter evolutions are overlaid on the count light curve. In the Appendix, we further provide the corresponding figures for all bursts (Figures A1–A3) as well as the figures for all bursts for the evolution of the spectral parameter β. Following the traditional classification, we use the notation that the parameter value is denoted as "hard" when referring to large values of both α and E_p, as opposed to soft values (low values of α and E_p). We categorize the evolution in the two main groups: those with a hard-to-soft (h.-t.-s.) pattern, that is, the parameters decrease independent of the rise and decay of the pulse, and those with a flux-tracking (f.-t.) pattern, that is, the parameters are correlated with the rise and decay of the flux with or without a time lag. In a handful of cases, other patterns are also identified. The classification of the pulses is given in Table 2 (columns 8–9).

Zoom In Zoom Out Reset image size

For the E_p evolution (see Figure A1 and column 8 in Table 2), we find that the hard-to-soft and flux-tracking patterns are the two dominant patterns in our multipulse sample. We find that about two-thirds (63/103 = 61%) of the pulses show a flux-tracking pattern and about one-third (32/103 = 31%) of the pulses exhibit a hard-to-soft pattern. Other evolution patterns are rarely observed. For instance, we only identify two cases showing the hard-to-soft-to-hard pattern (e.g., P₄ in GRB 101014175), one case displaying the soft-to-hard evolution (P₁ in GRB 170207906), and three cases exhibiting a hard-to-soft followed by a flux-tracking evolution within a pulse (e.g., P₁ in GRB 091127976). Here we note that previous investigations found that about two-thirds of cases have hard-to-soft behavior, while a smaller fraction has flux-tracking behavior (e.g., Ford et al. 1995; Crider et al. 1997; Lu et al. 2012; Yu et al. 2019). This fact is at odds with the observation in our sample. This is further discussed in Section 5.2.

For the α evolution (see Figure A2 and column 9 in Table 2), we find, similarly, that the hard-to-soft and flux-tracking patterns dominate. The flux-tracking pattern accounts for 65% (67/103) of the pulses, while the hard-to-soft pattern accounts for 25% (26/103) of the pulses. These two patterns thus account for 90% of the pulses. Among the rest of the pulses, we find that three cases have a soft-to-hard evolution (e.g., P₁ in GRB 110625881), and one case has a hard-to-soft-to-hard evolution (P₄ in GRB 101014175). Moreover, we also identify two cases showing a "flat" (or weak rise) behavior throughout the pulse (e.g., P₁ in GRB 120129580), and two cases have no clear trend at all (e.g., P₁ in GRB 120728434), which is not found in the E_p evolution.

We note that the E_p and α evolutionary patterns during a single pulse are not necessarily the same. In roughly half of the pulses (51% = 53/103), the E_p and α evolution are classified to have the same pattern, while 49% (50/103) of the pulses do not have the same pattern.

Finally, the patterns of the spectral evolution for E_p(t) and α(t) can vary from pulse to pulse within a burst. We find that in only about half of the pulses, the patterns of the E_p(t) evolution between two adjacent pulses are the same (55% = 35/64), and the rest 45% (29/64) are the inconsistent cases. Similarly, for α(t), we also find that for about half of the pulses, the spectral evolution among different pulses shares a similar pattern, accounting for 58% (37/64) of the pulses, while the inconsistent cases account for 42% (23/64) of the pulses. However, there is no significant variation in the fraction of pulses that are classified as having tracking behavior between the six pulse groups, which all have around a two-thirds fraction (Table 5 and Figure 5).

Zoom In Zoom Out Reset image size

We now turn to investigating the correlation between the following parameter pairs: ($\mathrm{log}F$, α), ($\mathrm{log}F$, $\mathrm{log}{E}_{{\rm{p}}}$), and (α, $\mathrm{log}{E}_{{\rm{p}}}$). Examples of these relations are provided in the lower panels of Figure 4, where functional fits are made. The leftmost panel shows the fit of F = F₀ e^{k α} , where k ∼ 3.37 ± 0.49 (a typical value; Ryde et al. 2019), the middle panel shows the power law between F and E_p with index 1.50 ± 0.15 (a typical value; Borgonovo & Ryde 2001), and the rightmost panel shows a fit to $\alpha ={k}_{2}\mathrm{ln}({E}_{{\rm{p}}}/{E}_{0})+{\alpha }_{0}$, where k₂ = − 2.01 ± 0.20. In the Appendix, we further provide the corresponding figures of these correlations for the full sample (Figures A5–A7).

We visually inspect the correlations and classify them according to the scheme in Yu et al. (2019), who classified them into three behaviors. The first behavior is a monotonic relation, defining type 1. It can be divided into three categories: type 1p, monotonic positive correlation; type 1n, monotonic negative correlation; and type 1f, flat relation. The second behavior has two piecewise monotonic relations combined at a break point. This behavior is defined as type 2 and is divided into two subcategories, either a concave (type 2p) or a convex (type 2n) function. No clear trend is classified as type 3. The classification is given in Table 2.

The strengths of the correlations are given by the Spearman's rank r. Strong correlations have r > 0.7, and weak correlations have r < 0.4. For both the F–α and F–E_p relations, around half of the pulses have strong correlations (47/103 and 54/103), and a quarter have weak correlations (27/103 and 26/103). In contrast, the α–E_p relations have the reverse properties: strong correlations in 21/103 and weak correlations in 57/103.

Among the pulses with strong correlations, we find that for the F–α relation, the vast majority has a positive monotonic relation (1p; 38/47), only a few have a negative relation (1n; 6/47). Similarly, for the F–E_p relation, a vast majority have a positive power law (1p; 45/51). Finally, for the α–E_p relation, negative and positive relations are equally common (1n; 9/21) and (1p; 8/21).

The three relations F–α, F–E_p, and E_p–α typically do not have strong correlations at the same time. However, in a few cases, they do: 13 out of the 103 pulses have r > 0.7 or r < −0.7 for all correlations at the same time. An example of a spectral evolution with simultaneous strong correlations is the single-pulse burst GRB 131231A, where all three relations have monotonic positive correlations (Li et al. 2019). In the following discussions, we summarize the correlations of all of the pulses independent of their r values.

3.3.1. Individual F–α Relation

3.3.2. Individual F–E_p Relation

3.3.3. Individual α–E_p Relation

In the upper panel of Figure 6, we plot the values of ${\alpha }_{\max }$ and the corresponding value of E_p, with one data point from every pulse. Comparing with the limiting lines for SCS and the NDP, we find that 67% (79/103) of all pulses have at least one data point above the SCS line and 21% (22/103) of all pulses have a data point above the NDP line.

The number of pulses that are above the NDP and SCS lines for the sequence of pulses is shown in the lower panel of Figure 6, together with the corresponding fractions. The first two pulses in a burst have a larger fraction above the SCS line. There is a clear decrease in later pulses. There are, however, a couple of bursts in which the hardest spectrum occurs at late times. This fact is reflected in the large fraction for P₅ in Figure 6. Examples of these bursts are GRB 121225 and GRB 140810, where the largest α occurs around 50 s after the trigger.

As mentioned above, the spectra with $({E}_{{\rm{p}}},{\alpha }_{\max })$ values lying above the SCS line can be expected to have a higher probability of having a photospheric origin. However, how big this probability is can only be answered by a model comparison of the physical models, for instance, through Bayesian evidence. In any case, for spectra close to the SCS line, a model comparison will be inconclusive, since both models can produce similar spectra (over the observed energy range). Acuner et al. (2020) investigated this point quantitatively and found that a photospheric preference can be claimed, with great confidence, only for spectra with α $\gtrsim$ −0.5, as long as the data have a high significance. To illustrate this point, we perform the Bayesian model comparison for two example time bins in GRB 150330 (at 1.4 and 137.0 s). We calculate the Bayesian evidence for the NDP spectrum, Z_NDP, and the evidence for the the slow-cooled synchrotron spectrum, Z_SCS. As shown in Acuner et al. (2020), a log-evidence difference of $\mathrm{ln}{Z}_{\mathrm{NDP}}/\mathrm{ln}{Z}_{\mathrm{SCS}}\gtrsim 2$ indicates that the NDP spectrum is preferred, while if the ratio is less than −2, the SCS spectrum is the preferred model. For the 1.4 s time bin, α = −0.24, and we find that $\mathrm{ln}{Z}_{\mathrm{NDP}}/\mathrm{ln}{Z}_{\mathrm{SCS}}=33.6$, strongly favoring a photospheric origin. Correspondingly, for the 137.0 s time bin, which has α = −0.9, we find that $\mathrm{ln}{Z}_{\mathrm{NDP}}/\mathrm{ln}{Z}_{\mathrm{SCS}}\,\lt -78.9$, strongly favoring an SCS origin (in the comparison between these two specific models). This shows again that the α value can be used to make an approximate model comparison in order to identify the preferred model, which is sufficient for the present study. A full model comparison investigation based on Bayesian evidence is beyond the scope of this paper.

We therefore identify spectra that have ${\alpha }_{\max }\gt -0.5$ according to the Acuner et al. (2020) criterion. These spectra are denoted by Ph in Table 2 and shown by the blue bars in the lower panel of Figure 6. The fraction of these spectra steadily decreases for subsequent pulses. For the first pulse in a burst, nearly half of all pulses (18/39) pass this very stringent criterion. The fraction decreases by approximately half for every following pulse.

We now examine the spectral evolution and correlations for the pulses that are compatible with photospheric emission (Ph in Table 2). The purpose is to assess whether they have different characteristics of their spectral properties.

For the photospheric pulses, the majority are classified as flux-tracking or hard-to-soft for both E_p and α. For the E_p evolution, the tracking pattern is found in about half, or 47% (20/43), which is only somewhat lower than in the full sample. However, for the α evolution, flux tracking is still dominant with 63% (27/43), similar to the full sample.

Turning to the parameter correlations, we find that 47% (20/43), 51% (22/43), and 26% (11/43) have a strong correlation (r > 0.7) for the F–α, F–E_p, and α–E_p relations, respectively. These fractions are similar to the ones for the full sample, which indicates that the strength of the correlation does not depend on whether the pulses are photospheric or not.

We again consider the two main types of correlations for each pair of parameters. For the F–α relation, the two main types in the full sample were positive (1p) and negative (1n) monotonic correlations. For the photospheric pulses, these fractions are 60% (26/43) type 1p and 21% (9/43) type 1n. For the F–E_p relation, the two main patterns are type 1p and the convex relation (type 2p), and the photospheric pulses have 60% (26/43) type 1p and 16% (7/43) type 2p. Finally, the α–E_p relation has 21% (9/43) type 1p and 35% (15/43) type 1n. Compared to the full sample (Section 3), these fractions are not significantly different.

The bursts in our sample were specifically selected, since they have multiple pulses. Here we address the question of whether or not the spectral properties of the pulses in our sample are similar to those of the single-pulse bursts. To do this, we compare the results in Yu et al. (2019), who studied a sample of 37 single-pulse bursts.

In Figure 7, we compare the distributions of duration T₉₀, peak flux, fluence, and the time-integrated spectral shape (E_p, α, and β from the Band function fit). Our sample is shown by magenta lines, and the Yu et al. (2019) sample is shown by the green lines. All of these parameters are collected from the Fermi/GBM catalog. Both the peak flux (1024 ms timescale) and the fluence are taken over the 10–1000 keV range. The average values and standard deviations (1σ error) are presented in Table 6. Both the peak flux and the fluence are approximately double for the multipulse bursts. Also, the α distribution is shifted to softer values for the multipulse bursts. On the other hand, the E_p and β values are similar between the samples. The corresponding probabilities from K-S tests are 0.02, <10⁻², <10⁻², 0.86, <10⁻², and 0.19 for t₉₀, peak flux, fluence, E_p, α, and β, respectively.

Zoom In Zoom Out Reset image size

In Figure 8, we compare the distributions of the time-resolved spectral parameters. In order to make a consistent comparison, we compare the results from using the CPL function throughout for both samples. This choice is based on the fact that in the Yu et al. (2019) sample, the CPL model is the best model. Likewise, we find that in only 29% of the time-resolved bins is the Band function a better choice (Section 2.6). We find that the average distribution of α is softer and the flux is larger for the multipulse sample. However, the E_p distributions are similar. These results are similar to the time-integrated spectra comparison. If we instead particularly identify the distribution of P₁, we find that α distributions have similar peak values. On the other hand, the flux of P₁ is double the corresponding flux for the single-pulse bursts.

Zoom In Zoom Out Reset image size

We also compare the frequency of the parameter relations between the samples. Most proportions between the patterns are similar. However, for the α–E_p relation, there is a notable difference. In our sample, the proportions are (72%) 74/103 for type 1 and (20%) 21/103 for type 2, while for the single-pulse sample, the proportions are (32%) 12/38 for type 1 and (45%) 17/38 for type 2. There is still a similarity in that the two pattern types are frequent in both samples; there is no strong dominance. However, in the multipulse sample, the type 1 pattern is relatively more frequent.

We find that the largest fraction of pulses in our sample follow a tracking pattern, which is at odds with earlier investigations that instead found that the hard-to-soft evolution accounts for about two-thirds and the flux-tracking evolution only accounts for one-third of the observations (e.g., Ford et al. 1995; Lu et al. 2012; Basak & Rao 2013; Yu et al. 2019).

Lu et al. (2012) argued that if there is a significant overlap between pulses, the classified pattern might be affected. For instance, an apparently flux-tracking pulse might be a consequence of two overlapping hard-to-soft pulses. In our sample, many pulses are slightly overlapping, which means that in the transition period between the pulses, there are contributions from both pulses. Therefore, this might affect our results. To investigate this point, we particularly study the "gold" sample (see column 5 in Table 2). These nine bursts ¹² have 22 pulses that are clearly separated by periods of no emission. Among these, the flux-tracking pattern accounts for 59% (13/22) of the pulses, while the hard-to-soft pattern only accounts for 36% (8/22) of the pulses. These fractions are similar to the ones from the full sample (see Section 3.2). We conclude, therefore, that our sample is not greatly affected by the overlap between pulses.

Another possible reason for the difference is the binning methods used. We use a combination of BBlocks and significance and only consider the most significant time bins. The advantage of our method is that the substantial variations in the light curve are captured and the spectral fits are reliable. On the other hand, we might miss some of the rising phases of the pulses. Including a larger fraction of the pulse duration by using less significant data (S < 20) could change the apparent pattern of evolution. The spectral evolution of the added periods might indeed be different, but more seriously, the spectral fits are less reliable and could therefore give misleading results.

A further uncertainty is the inherent problem in this type of classification. A subjective decision needs to be made on where a pulse starts and ends. This can affect the classification. Another point is the fact that the tracking behavior, in many cases, involves a time lag, both positive and negative. If this lag is large compared to the pulse size, the pulse could be classified as a hard-to-soft (or soft-to-hard) pulse instead. An example is the E_p pattern for GRB 120328, where the tracking pattern depends on one data point during the rise phase. The same issue happens for α in other cases. ¹³ These uncertainties can affect the classification of different samples. However, these cases are not sufficient to explain the whole discrepancy.

Our results thus indicate that the flux-tracking pattern is the prevalent pattern for α and E_p, at least around the peak of the pulse, where the significance is the largest. This is also consistent with the global F–E_p relation in the middle panels in Figure A9.

Based on the analysis of the change in α and E_p (Section 3.1 and Figure 2) alone, it was suggested that there are two categories of pulse groups. The first one, consisting of P₁ and P₂, has a constant E_p, while α softens. Within the second category, there is not much change, but they are all different from the first two pulses. Moreover, we found that the initial pulses had different spectral properties (frequency of different parameter correlation) as a group (Section 3.3). These distinctions lead to the possibility that they are due to different emission mechanisms.

In the simplest internal shock model, the late pulses (P_i with i > 1) occur above the photosphere and hence must have a synchrotron origin. We do not see such a clear distinction, indicating that at least the simplest version of the internal shock model does not represent what we see.

For the photospheric scenario, the α value reflects the energy dissipation and photon production below the photosphere (e.g., Pe'er et al. 2006b; Vurm & Beloborodov 2016). The E_p instead is mainly set by properties at high optical depths. Therefore, if the first two pulses are photospheric, the fact that E_p is similar while α varies could be due to similar properties close to the central engine but a varying amount of turbulence in the flow causing the dissipation. The variability and pulse structures of the photospheric emission have been reproduced by numerical simulations of a jet passing through the progenitor surrounding. For instance, López-Cámara et al. (2014) showed that even with a steady central engine, a light curve with a pulse structure and large variability arises. The main cause is the Rayleigh–Taylor instabilities that arise in the contact between the layers of jet and the surrounding progenitor material. This leads to a variable amount of mixing between the layers and thereby a variable baryon load of the jet, which has a direct influence on the radiative efficiency and spectral shape (see, e.g., Rees & Mészáros 2005; Gottlieb et al. 2019).

On the other hand, the pulse groups, P₃–P₆, all have α ∼ −1 (Figure 2). This could be explained by synchrotron emission from electrons that are marginally fast-cooled (Daigne et al. 2011; Yu et al. 2015; Geng et al. 2018). It can thus be argued that the first two pulses are typically photospheric, while the rest are due to synchrotron emission. If this interpretation is correct, one would expect that the two categories would have distinct and different spectral properties. There are several points in the analysis above that, therefore, do not support this interpretation. First, the values of ${E}_{{\rm{p}},\max }$ are similar for all of the pulse groups (Figure 3); there is no clear distinction between the pulse groups. Second, there is no clear distinction between pulse groups when it comes to the spectral evolution (Section 3.2) and types of correlations (apart from the initial pulses; Section 3.3). Third, the average pulse durations are also similar in all pulse groups. All of these points indicate that there is no drastic change in emission pattern between the pulses (apart from α).

On the contrary, the properties of the last four pulse groups can also be interpreted as photospheric emission. Previous studies have shown that pulses that, beyond any doubt, are photospheric all have significant spectral evolution (e.g., Ryde et al. 2019). In particular, at the end of such pulses, α values down to −1 and below are common, which is interpreted as the result of subphotospheric dissipation with varying jet properties (e.g., Ryde et al. 2010, 2011). Consequently, pulses with α ∼ −1 could be the result of dissipative photospheres throughout the pulse. Indeed, the observation that the last four pulse groups all peak at around α ∼ −1 is in line with the theoretical expectations of fully dissipative photospheres, i.e., flows where there is no strong limitation, not on the soft photon production deep in the flow or the amount of dissipation in the parts of the flow where the spectrum is formed (Beloborodov 2010; Vurm et al. 2013).

Even though this line of argument makes the subphotospheric dissipation model appealing for most of the pulses, some admixture of synchrotron pulses is expected. The most prominent example of this is GRB 190114C, which, apart from a very hard spectral component, showed a clear afterglow emission component already during the prompt phase (Ajello et al. 2020). This proved that a synchrotron component was present during the prompt phase (see also Axelsson et al. 2012; Iyyani et al. 2013). Further early examples of such a suggestion are given in Ryde (2005), Ryde & Pe'er (2009), and Guiriec et al. (2011) and more recently in Burgess et al. (2019 and Wang et al. (2019a, 2019b)).

Many of these synchrotron components are observed at late times. Combined with the fact that the fraction of (certain) photospheric pulses decreases with pulse number (Figure 6), emission compatible with synchrotron appears to prevail mainly at the end of the prompt phases. Different possibilities exist for synchrotron emission to arise toward the end of the prompt phase. One example is GRB 150330, for which Li (2019a) suggested that the jet composition changes from a baryon-dominated flow during the first pulse to a Poynting flux–dominated flow during the rest of the burst, producing synchrotron emission (see also Zhang et al. 2018). Alternatively, during the early phases of the afterglow, shocks due to interaction between the jet and slower-moving material ahead of the jet will also produce efficient synchrotron emission (e.g., Duffell & MacFadyen 2015).

However, from the GBM data alone, it is not possible to make a firm conclusion as to whether a soft pulse is due to synchrotron or photospheric emission subject to dissipation below the photosphere, since they are indistinguishable in many cases (e.g., Acuner et al. 2020). Simultaneous data from other wavelengths can be useful in some cases (e.g., Ravasio et al. 2018; Ahlgren et al. 2019; Ajello et al. 2020) or polarization measurements (e.g., Sharma et al. 2019).

In Figures A1–A4, we provide the evolution of the spectral parameters E_p, α, vF_v flux, and β for all the bursts based on the best models defined in Section 2.6.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Figures A5–A7, we provide the F–α, F–E_p, and α–E_p relations for all the bursts based on the best models defined in Section 2.6.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Figure A8, we provide the results of the parameter distributions, separated by two empirical photon models (Band and CPL). The overall parameter distributions we studied include four parameters (α, β, E_p, and F) in the Band model and three parameters (α, E_p, and F) in the CPL model. The corresponding average values and standard deviations from the best Gaussian fit for each distribution are presented in Table A1.

Zoom In Zoom Out Reset image size

The main results of parameter distributions (α, E_p, and F) among pulses obtained from a single model (Figure A8), either the Band or the CPL model, are consistent with the finding in the best model (Figure 2), as we discussed in the main text. For the β distribution, we find a bimodal distribution for each time-series pulse sample, as well as the global sample, where the harder peak is at ∼−2.3 and the softer peak is at ∼−6.1. The two peaks are basically the same for all different time-series pulses. We also find that the β indices are typically softer (about half of the Band spectra; the obtained β indices cannot be well converged) than some previous catalogs, whose analysis is based on the frequency analysis method, but consistent with the results found in Yu et al. (2019), who also adopted a fully Bayesian method but for a single-pulse sample.

We show global pulsewise parameter relations for our full sample in Figure A9. We present each parameter relation based on the two different photon models (the Band and CPL): F–α relation based on the Band model (first left panel), F–α relation based on the CPL model (first right panel), F–E_p relation based on the Band model (second left panel), F–E_p relation based on the CPL model (second right panel), α–E_p relation based on the Band model (third left panel), α–E_p relation based on the CPL model (third right panel), and α–β relation based on the Band model (bottom panel).

Zoom In Zoom Out Reset image size

We find that there is no significant difference in the global parameter relations between the Band and CPL models. In addition, the F–E_p Golenetskii relation based on the Band model shows a stronger monotonic positive correlation than the CPL model. For the F–α relation, a cluster of data points that significantly deviate from the peak of the distribution (probably mainly contributed by type 2 in the individual parameter relations) mainly come from the early (P₁) pulse. In short, the global F–α relation shows a monotone positive correlation following a break behavior. Such a break may be originated from thermal emission. The global F–E_p relation displays type 1p behavior, while the global α–E_p shows an anticorrelation type 1n behavior. Finally, we do not find a clear trend in the global α−β relation.