Standard GRB Spectral Models"Misused"?

The standard model characterizing the gamma-ray burst (GRB) spectrum invokes a four-parameter empirical function, the so-called the BAND model. An alternative model named cutoff power law (COMP) implements a power law with an exponential cutoff. These functions achieve almost equally good fits on observed spectra, and are adopted in nearly all of the GRB literature. Here, we reanalyze the sample defined in Li.et al.,2021,ApJS,254,35 (39 bursts including 944 spectra). We classify the spectra by two methods: (1) checking their corner-corner plots of the posteriors to determine well-constrained $\beta$ (BAND-better) and unconstrained $\beta$ (COMP-better) categories; and (2) defining the four groups by difference of the deviance information criterion (DIC). We find inconsistent peaks of the parameter distributions between the BAND-better spectra ($\alpha=-0.64\pm0.28$ and $\rm log_{10}(E_{\rm p})=\rm log_{10}(191)\pm0.41$) and the COMP-better spectra ($\alpha=-0.96\pm0.33$ and $\rm log_{10}(E_{\rm p})=\rm log_{10}(249)\pm0.40$). With the statistically preferred model and vice versa the misused model defined based on DIC statistics, we also find that the fitted parameters obtained by the misused model (COMP) significantly deviate from those obtained by the statistically preferred model (BAND). This means that if a spectrum is statistically preferred, described as the BAND, applying COMP to derive the spectral parameters will prominently deviate from their intrinsic shape, therefore affecting the physical interpretation. Our analysis indicates that the better or statistically preferred model should be duly examined during GRB spectral analysis. In addition, the $\beta$ distribution exhibits a bimodal structure containing the BAND-better and COMP-better spectra, respectively, implying that BAND and COMP both may have physical origin.


INTRODUCTION
The standard approach to characterize the observational gamma-ray burst (GRB) spectral properties invokes a fourparameter empirical function known as the BAND model (Band et al. 1993). The photon number spectrum of BAND is defined as where A is the normalization factor in units of ph cm −2 keV −1 s −1 , E piv is the pivot energy always fixed at 100 keV, E 0 is the break energy correlated with the peak energy of νF ν spectrum (assuming β < −2) by E p = (2 + α)E 0 , α and β are the low-energy and high-energy asymptotic power-law photon indices, respectively. The spectral indices (α and β) and the peak energy 1 (E p ) are typically distributed around α=-0.8 (below the break energy), β=-2.5 (above the break energy), and E p =210 keV, respectively. An alternative empirical approach involves a simpler function called the cutoff power-law (COMP, aka the Comptonized model) model. This approach is valid when the power-law index β is poorly constrained (having fairly large absolute values and large uncertainties; see, e.g. in Fig 1). The COMP function is recovered from the BAND function Li as β tends to −∞. The COMP function is given by where the peak energy E p of the νF ν spectrum is related to the E c through E p =(2+α)E c . The physical origins of these empirical functions, however, have yet to be identified, although they have been the most widely used to fit GRB spectra. Neither BAND nor COMP functions correspond to an explicit emission mechanism. Whether these models are due to thermal or northermal emission is highly debated, depending on the slope values of their spectral parameters. Physically, the leading mechanisms for interpreting GRB prompt emission invoke either nonthermal photons originating from synchrotron emission (or inverse Compton scattering) (e.g., Meszaros et al. 1994;Rees & Meszaros 1994;Zhang & Yan 2011;Geng et al. 2018;Meng et al. 2018Meng et al. , 2019Li et al. 2019) or Comptonized quasi-thermal photons associated with photosphere emission (e.g., Thompson 1994;Pe'er et al. 2007;Ryde et al. 2010;Ruffini et al. 2013;Li 2019aLi ,b, 2020Xue 2021). The fast-cooling (α=-3/2) and slow-cooling (α=-2/3, so-called the line of death of synchrotron emission, Preece et al. 1998) synchrotron emission predicts two different values of α, whereas photosphere models predict much harder values of α (e.g., above α=-2/3). For instance, Acuner et al. (2020) argued that the spectra that prefer the photospheric model all have low-energy power-law indices α ∼>-0.5, as long as the data has a high significance. Therefore, applying these empirical models to the GRB spectral analysis plays an important role in identifying the GRB radiation mechanism, and investigation of spectral parameters, and therefore, will shed light on our understanding of GRB physics (e.g., Dai et al. 2006;Kaneko et al. 2006;Zhang et al. 2006;Gruber et al. 2014;Yu et al. 2016;Ruffini et al. 2018;Li 2019a,b;Li et al. 2019;Li & Zhang 2021;Xue 2021;Moradi et al. 2021, and references therein).
In general, we can apply either time-integrated or time-resolved spectral analysis to study the spectral properties of a GRB. Several spectral catalogs of GRBs exist in the literature based on either the time-integrated analysis (e.g., Goldstein et al. 2012;Gruber et al. 2014) or the time-resolved analysis (e.g., Yu et al. 2016;. The time-integrated spectrum represents the average spectral properties since the entire period of emission is treated as a single time bin. However, GRB prompt emission is well known to have strong spectral evolution (e.g., Kaneko et al. 2006;Yu et al. 2019;, and references therein), which requires the more detailed time-resolved spectral analysis (treating the whole period of emission as multiple timing bins, and spectral analysis is therefore performed on each timing event individually, e.g., Yu et al. 2016Yu et al. , 2019. Several early GRB spectral catalogs make use of the frequentist approach (e.g., Kaneko et al. 2006). In recent years, a fully Bayesian analysis method has been increasingly developed. For example, time-resolved spectral catalogs based on such a fully Bayesian analysis method for single-pulse bursts (Yu et al. 2019) and multi-pulse bursts  have been created. In the Bayesian analysis, Bayesian inference is used to account for relevant prior information and the resulting posterior probability distributions of parameters are obtained by the Markov Chain Monte Carlo (MCMC) iterations.
Phenomenologically, the BAND-like spectrum with well-constrained model parameters is typically observed in the time-integrated spectral analysis, while the simpler COMP-like spectrum is commonly observed in the time-resolved spectral analysis. This is because time-resolved spectral properties typically do not have good high-energy photon statistics, and therefore, the high-energy spectral index β for time-resolved spectra usually cannot be well evaluated due to the small number of photons available.
It is important to stress that the difference in fitting by BAND or COMP functions is not fully examined when performing the time-resolved analysis of large samples (e.g., Kaneko et al. 2006;Goldstein et al. 2012;Yu et al. 2019;. Moreover, in some articles COMP is applied throughout without a comparison with other models since the COMP is usually preferred for the majority of the time-resolved spectra. We have doubts about the statistical conclusions and the physical implications generated from possible misused model. Therefore, we dedicate this article to examining the deviation of spectral fittings between BAND and COMP. We wish to answer the question: Is the impact on parameters significant if misusing a model? Do BAND and COMP both have physical backgrounds? Here we reanalyze the sample (39 bursts including 944 spectra) defined in  to examine the spectral properties statistically of these two standard spectral models. This paper is organized as follows. The methods are presented in §2. The detailed results are summarized in §3. The discussion and conclusion are presented in §4 and §5, respectively. The convention Q = 10 x Q x is adopted in cgs units throughout the paper. The standard Λ cold dark matter (CDM) cosmology with the parameters H 0 = 67.4 kms −1 Mpc −1 , Ω M = 0.315, and Ω Λ = 0.685 are adopted (Planck Collaboration et al. 2018).

Sample Revisited
The Gamma-ray Burst Monitor (GBM; 8 KeV-40 MeV, Meegan et al. 2009), and the Large Area Telescope (LAT; 20 MeV-300 GeV, Atwood et al. 2009), are the two instruments on board Fermi providing unprecedented spectral coverage for seven orders of magnitude in energy. Fermi -GBM, together with Fermi -LAT, has been triggered by more than 2000 bursts since its launch in 2008. Here, we revisit the sample defined in . The sample is collected from the Fermi -GBM burst catalog published at HEASARC 2 , and it focuses on well-separated multi-pulse GBM-detected bursts. It consists of 39 bursts, 117 pulses, and 1228 spectra. There are two reasons that we included the sample in this task. First, all the spectra in the sample were selected to have a high statistical significance in order to allow us to perform a detailed time-resolved spectral analysis and ensure that the spectral fits are well determined, this is the key point. Second, the prompt-emission light curves of GRBs typically exhibit irregular, multi-pulse temporal profiles.
The sample selection in  includes the following main steps: (1) The first is to visually inspect the light curves for each burst that was observed by Fermi-GBM during its first 11 yr of mission (more than 2000 bursts), and about 120 bursts that have well-separated multi-pulse features are roughly identified; (2) The second is to capture the variations of the Time-Tagged Events (TTE) light curve and divide the light curve into time segments by following the Bayesian blocks (BBlocks; Scargle et al. 2013) algorithm, and the significance (S; Vianello et al. 2018) for each time bin was also calculated; (3) The third is to select at least two pulses in each burst whose individual pulse light curve has at least four time bins with high significance (S ≥ 20); 3 , hence, the final sample was defined (39 bursts, 103 pulses, and 944 spectra); (4) The final goal is to obtain the best spectral parameters by adopting a fully Bayesian analysis using the MCMC method and performing both the BAND and COMP functions to fit all the spectra, respectively. For information on the data procedure, including the burst, detector, source, and background selections, light curve binning method, sample definition, and Bayesian and MCMC spectral fitting approaches; please refer to ; Li (2019a); Li & Zhang (2021) for more details.
2.2. High energy power law, β, and the "better" models In reality, for a given spectrum, in order to determine which one (BAND or COMP) is better, one needs to check whether a well-constrained β can be determined. If β is not well constrained in some cases, there are two possibilities. Firstly, lack of photons in the analysed bins (e.g., S < 20), so that the spectral fit cannot be well determined. Secondly, the number of photons in the analyzed bins is sufficient (e.g., S ≥ 20), but the model that better characterizes the spectral shape is indeed the COMP. Our sample defined in  with S ≥ 20 rules out the first possibility. We therefore inspected all the posteriors of the BAND spectra to check their β indices. If a well constrained β is clearly identified by a certain spectrum, the BAND model is considered as better, otherwise the COMP model is better. Under these criteria, all the spectra can be identified into two categories: • BAND-better spectra: All the spectra selected in this category are identified with a well-constrained β, indicating that the BAND model is indeed better. It contains 35% of the total number of spectra.
• COMP-better spectra: All the spectra in this category are identified with an unconstrained β, implying that the COMP model is better. This is 65% of the total number of spectra.
In Figure 1, the left panel shows two-dimensional corner-corner plots of the spectral parameters using the Bayesian MCMC method used to perform the BAND fit. The spectral data is obtained from one time bin (between 24.215 and 25.597 s) from GRB 171227, and an unconstrained β is clearly identified from the posterior density map. While that of the COMP fit for the same spectral data is displayed in the right panel of Figure 1. For comparison, we also 2 https://heasarc.gsfc.nasa.gov/W3Browse/fermi/fermigbrst.html 3 Although the BBlocks method can better capture the intrinsic variation of light curve (e.g., , the time bins created by such a method usually have varied signal-to-noise ratios. This means that we cannot ensure that there are enough photons in each time bin in order to establish a reliably spectral fit. In order to ensure that the spectral fits can be well determined, a relatively high statistical significance for each selected time bin is needed. On the other hand, the threshold levels of statistical significance required by different spectral models may also be different. Practically, more complicated models (with more free parameters) require more signal photons in order to establish a reliable fit result. Therefore, a threshold level of S ≥ 20 is typically used for the BAND model while S ≥ 15 for the COMP model since the BAND model (A, α, β, Ep) has one more free parameter than the COMP model (A, α, Ec). The sample defined in  adopted S ≥ 20 to select the time bins, which is enough to study both two models.
present the same plots using another time bin (between 24.215 and 25.597 s) from GRB 171227 in Figure 2, where a well-constrained β is clearly identified in the BAND fit. In total, the fractions of the constrained-β and unconstrained-β spectra are 35% and 65%, respectively. This suggests that the majority of the spectra (two-thirds) can indeed be better fitted by the COMP, confirming the previous similar findings (Yu et al. 2019;).

Statistically Preferred Models Determined by Information Criteria
In practice, a more common approach for model comparison is by using information criteria, such as Akaike Information Criteria (AIC), Bayesian Information Criteria (BIC), and the Deviance Information Criterion. The Bayesian analysis and MCMC method are fully applied in this work, the deviance information criterion (DIC, Spiegelhalter et al. 2002;Moreno et al. 2013) is computed to compare models, it is defined as DIC=-2log[p(data|θ)]+2p DIC , whereθ 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). The values of the difference between the BAND's and the COMP's, defined as ∆DIC = DIC BAND -DIC COMP , can be used to indicate the preferred one. For each individual spectrum, a negative DIC value indicates that the observational data favors a Band-like spectrum. All of the spectra can be separated into the following groups using different threshold levels based on DIC statistics for Bayesian models (e.g., Gelman et al. 2014;Pooley & Marion 2018), in which the BAND-preferred or COMP-preferred spectra could also be determined and gathered.
• Group I: ∆DIC <=-10. BAND model is statistically preferred. This group contains 29% of the spectra, as shown in Figure 3.
• Group II: -10< ∆DIC <=-5. BAND model is still statistically preferred, but is not as strong as Group I. This group contains 11% of the spectra (Figure 3).
• Group IV: ∆DIC >0. COMP model is statistically preferred, and is stronger than Group III. This group contains 38% of the spectra (Figure 3).
Based on the BAND-better and COMP-better spectra defined in §2.2, we then check the fractions of the spectra with or without a constrained β for each DIC-defined group. Groupwisely, the corresponding fractions are [93%, 7%], [58%, 42%], [5%, 95%], and [1%, 99%] for Group I, Group II, Group III, and Group IV, respectively. In Group I, we find that for almost all the spectra (up to ∼93%) a well-constrained β can be clearly identified. However, there are very few spectra showing a well-constrained β in both Group III and Group IV. These results suggest that Band-like spectra dominate Group I while COMP-like spectra dominate Group III and Group IV. Interestingly, we also find that in Group-II, both Band-like (58%) and COMP-like (42%) spectra, are almost identical. The results also suggest that our two methods of classifying GRBs are consistent, but one needs to consider a relatively large DIC value (e.g. -5 is good and -10 is perfect, these values are in agreement with previous works (Acuner et al. 2020;.

RESULTS
Before we move forward, a few remarks need to be made here. Firstly, we focus on the two most widely used models for GRB spectra (BAND and COMP), and we miss several other models (e.g., power-law model, smooth broken power-law model, and the BETA model) that were used in the previous catalogs (e.g., Kaneko et al. 2006). In some cases, these models should be able to fit the spectra better than the BAND or COMP that we used in this task. For instance, if a break's energy lies outside the detector passband, or the source photon signal beyond the break energy is weak enough so that the break energy cannot be well determined. In such cases, the simpler power-law model (see one recent work, Tang et al. 2021, for instance) is superior to the other, more complicated models (having more parameters). As such, the models we used would be the better ones that characterize GRB intrinsic spectra, rather than the best ones. Secondly, our analysis is based on a sample of well-separated multi-pulse GBM-detected bursts and a criterion of statistical significance S > 20 was used to select the bright spectra for each individual burst. These may be causing some bias in our analysis results. Lastly, compared to previous catalogs (e.g., Kaneko et al. 2006;Yu et al. 2016) that used the χ 2 method to statistically compare the models and determine the best-fit model for each individual spectrum, our analysis is based on a fully Bayesian analysis approach using the MCMC method and we used the information criteria to compare the models. Unlike the χ 2 method involving a different degree of freedom in different models and resulting in the comparison not being straightforwardly performed, the information criteria (e.g., DIC statistics) that we used in this task may easily offer a straightforward comparison among different models because penalty factors for overfitting of more complex models have also been taken into account. Based on such a Bayesian analysis and MCMC spectral fit method, we may also be able to select the better models more straightforward by inspecting their posterior distributions from MCMC sampling, as compared to some previous studies that invoke a more complicated selection method to determine their good class of parameters (e.g., Kaneko et al. 2006;Yu et al. 2016).

Statistically Preferred Model Misused
In this task, our primary interest is to assess the effect of misuse of the models on fitting results. For instance, for a given spectrum that can be better (or statistically preferred) described by the BAND model, are the spectral parameters obtained from the simpler COMP model still consistent with that from the BAND model? Vice versa. To better address this question, we define (1) the BAND-to-BAND case (Column 6 in Table 1): a statistically preferred BAND model is used for a given Category (defined in §2.2) or Group (defined in §2.3); (2) the COMP-to-COMP case: a statistically preferred COMP model is used. Alternatively, if a statistically preferred model is not being used, in contrast, the model used is a statistically undesirable one. This may involve the better model being misused. We therefore also define (3) the BAND-to-COMP case: a statistically preferred model is BAND but COMP is misused, this invokes the case of underfitting; (4) the COMP-to-BAND" case: a statistically preferred model is COMP but BAND is misused, this invokes the case of overfitting.
We then investigate these misused cases using the following two typical examples, as shown in Figure 4. In the left panel of Figure 4, we present the spectral fit to a time bin (between 86.338 and 86.877 s in GRB 120728) using both BAND and COMP models. This spectrum can be statistically-preferred fitted by the COMP model, confirmed by the DIC statistics with a value of ∆DIC=1.2. The spectral parameters obtained by the COMP fit are α = −0.26 +0.23 −0.23 , and E p = 74 +16 −16 and obtained by the BAND fit are α = −0.17 +0.25 −0.27 , β = −6.46 +2.38 −2.38 , and E p = 71 +6 −5 . These results suggest that the fitted spectral parameters (α and E p ) between the COMP-to-BAND case all seem to agree. In the right panel in Figure 4, we present the spectral fit to another time bin (between 69.274 and 71.015 s in GRB 120728). The BAND model is the statistically preferred one that describes the spectral shape, which is also confirmed by the DIC statistics with a value of ∆DIC=-123.8. For comparison, we also fit the spectral data using the COMP model. For the BAND model: α = 0.09 +0.12 −0.12 , β = −2.52 +0.06 −0.06 , and E p = 76 +4 −4 . For the COMP model: α = −0.61 +0.05 −0.05 , and E p = 109 +7 −7 . We find a remarkable discrepancy in the spectral parameters between the BAND-to-COMP case (the COMP model with softer α values and higher E p , as compared with the BAND model). This is due to compensating for the lack of a high-energy spectral component in the model resulting in the underfitting.
We also test the misuse of models by fitting the simulated spectra, which are generated by the GBM Data Tools 4 . For the source spectra, we take the functional modeling of the BAND and set the initial model parameters of A=0.03, α=-0.55, E p =500, and β=-2.5, and those of the COMP with A=0.03, α=-0.55, and E c =500. For the background spectrum, we generate it using a phenomenological method that first to fit the background of GRB 210518A 5 then to simulate the background spectrum from the fitted parameters. The response matrix is taken from the first 10 seconds of GRB 210518A.
The simulated COMP-like spectrum fitted using BAND and COMP gives rise to a differential value of ∆DIC=0.4 (see the left panel in Figure 5). A small difference in DIC statistics suggests that adding a high-energy component to such a COMP-like spectrum does not significantly improve the fit. Whereas the simulated BAND-like spectrum results in large DIC statistics (∆DIC=-38.6) improvements (see the right panel of Figure 5), indicating a statistically significant high-energy power-law component. These simulated results are consistent with the findings using true spectral fittings as described above.

Comparisons of Parameter Distribution for β-Statistic-Based Categories
With the categories defined in §2.2, we present the distributions of the spectral parameters that are used to compare the BAND-preferred spectra and COMP-preferred spectra in Figure 6. The average values and the corresponding standard deviation obtained from the best Gaussian fits for parameter distributions are summarized in Table 1, including α (BAND and COMP), β (BAND only), and E p (BAND and COMP).
Before comparing the fitted parameters of BAND and COMP with different categories and groups, we caution that the low-energy spectral indices α obtained from BAND and COMP are asymptotic values rather than actual slopes, and therefore cannot be directly compared. In order to minimize the discrepancy, an effective α eff , computed at 25 keV (the BATSE 6 detector lower limit), was introduced by Preece et al. (1998). In the GBM observations, the lower limit of the detector is at 8 keV, which is much smaller than the BATSE, and the difference between the asymptotic values and the actual slopes can be negligible (Figure 7). The fit values of α, therefore, can be directly used for our further analyses.
For α distribution (the right panel of Figure 6), we find the BAND-better spectra and COMP-better spectra showing inconsistent peaks. The best fit gives α BAND = −0.64 ± 0.28 for the BAND-better spectra and α COMP = −0.96 ± 0.33 for the COMP-better spectra, with a difference between α of the BAND fits and the COMP fits of ∆α=0.32, where ∆α is defined as ∆α = α BAND − α COMP .
While if the BAND-better spectra are misused by the COMP fit, one has α COMP = −0.82±0.29, with a value of ∆α of 0.18. Likewise, if the COMP-better spectra are misused by the BAND fit, one has α BAND = −0.91 ± 0.31, with a value of ∆α of 0.05. Based on these results, several interesting results can be drawn: (1) a significantly statistical difference of the spectral parameters between BAND-better spectra and COMP-better spectra is found; (2) The deviation (the COMP model with higher E p and softer α indices than the BAND model) between the BAND-to-COMP case is much more significant than the COMP-to-BAND case. Similar results can also be found in the E p distribution (see the middle panels in Figure 6 and Column 9 in Table 1).
We use the Kolmogorov-Smirnov (K-S) test to assess whether the distributions change between the two distinct categories. The chance probability, P , determined by the K-S test, leads to a value of P K−S (α BAND , α COMP ) < 10 −4 for the α distributions and of P K−S (E BAND p , E COMP p ) = 1.06 × 10 −4 for the E p distributions between the BAND-better spectra and the COMP-better spectra, indicating that these distributions are indeed different from one another.
With separated well-constrained and unconstrained β categories, the β distributions show a single peak for each category, with the best fits giving β = −2.53 ± 0.39 for well-constrained β categories and β = −5.57 ± 0.90 for unconstrained β category, respectively (the lower panel in Figure 6 and Column 8 in Table 1).

Comparisons of Parameter Distribution for DIC-statistic-based Groups
We present the distributions of ∆DIC in Figure 3. We find that Group I, Group II, Group III, and Group IV can account for 29%, 11%, 22%, and 38% of the total number of the spectra, respectively. Based on these DIC-statisticbased groups, we then present the parameter distributions (α, E p , and β) by comparing BAND and COMP ( Figure  8) group-wisely. The parameter distributions obtained from the BAND model with the best Gaussian fit are shown by grey lines and those from the COMP model are shown by orange lines.
For α distribution (see the left panel in Figure 8 and Column 7 in Table 1), we find that α indices obtained from BAND are significantly harder than those obtained from COMP in each group 7 . More interestingly, such a statistically significant difference in parameters tends to be weaker during the transition from Group I (minimum-∆DIC) to Group IV (maximum-∆DIC).
Physically, we could diagnose the underlying physical mechanism through the distributions of α indices. This is because different theoretical models predict different distributions of α. The photosphere emission models usually associate with harder α indices while the synchrotron emission models typically relate to softer α indices. As pointed out by some previous works (e.g., Preece et al. 1998;Acuner et al. 2020), the low-energy index α is a good estimator for which model is preferred by the data. For example, the synchrotron emission explains the spectral indices with a limit, known as the "line-of-death", α = − 2 3 (Preece et al. 1998). Acuner et al. (2020) argued that the spectra that prefer the photospheric model all have low-energy power-law indices α −0.5. In Figure 8, the line-of-death of the synchrotron emission is indicated by the green lines for each group. As a result, we find that the fraction of the spectra with α beyond the synchrotron limit obtained from the BAND model is apparently greater than those obtained from the COMP model. Moreover, these fractions decrease for subsequent DIC-based groups. Groupwise, the corresponding fractions [BAND, COMP]  and Group IV, respectively. We also notice that the distribution of α has a smooth and well-defined Gaussian shape of at the "line-of-death", challenging the existence of the "line-of-death".
For E p distributions (the middle panel of Figure 8 and Column 9 in Table 1), unlike α distribution, we do not find a strong trend among the groups. Interestingly, the statistical significant difference in parameters tends to be weaker for subsequent DIC-based groups, resembling to the finding in the α distribution.
We present the groupwise β distributions in the right panel of Figure 8. Using the same data,  found a similar bimodal distribution based on the better model determined for each individual pulse, with the harder peak at ∼ -2.3 and the softer peak at ∼ -6.1. The results indicate that the BAND-better and COMP-better spectra should be mixed to compose the distributions, and the harder peak should be contributed by the BAND-better spectra while the softer peak should be contributed by the COMP-better spectra.
Separated by the DIC statistics, we group the spectra into four groups as defined in §2.3. Interestingly, we find that all of the groups show a single peak, but the peak is clearly shifted from Group I to Group IV with a hard-to-soft trend. The hardest peak is at ∼ -2.3 found in the Group I (minimum-∆DIC) (Figure 8). This value is the same as the harder peak of the bimodal distribution found in the pulse-wise categories , implying that the peak is dominated by the BAND-like spectra. Likewise, the softest peak is at ∼ -6.1 (Figure 8), which is found in Group IV (maximum-∆DIC). This value is the same as the softer peak of the bimodal distribution found in , suggesting the peak is dominated by the COMP-like spectra. However, the peak ( Figure 8) for Group II is β = −3.23 ± 0.74 while that for Group III is β = −5.01 ± 0.89, suggesting a mix of Band-like and COMP-like spectra.

Comparisons of Parameter Relations
In GRB physics, the study of parameter correlations is an open question, and it plays an important role in understanding the underlying physical processes and radiation mechanisms(e.g., Amati et al. 2002;Geng & Huang 2013;Srinivasaragavan et al. 2020;, and references therein).
We first compare the same spectral parameters between two models by plotting α COMP -α BAND (the left panels of Figure 9), E COMP p -E BAND p (the middle panels of Figure 9), and F COMP p -F BAND p (the right panels of Figure 9). We find that the α indices obtained from the BAND model are systematically harder than the ones obtained from the COMP model, particularly in the BAND-better spectra. However, this trend is weaker in the BAND-wise spectra as compared to the COMP-wise spectra, which is consistent with the finding based on parameter distributions as discussed in §3.3. Similar results are also found in the E COMP p -E BAND p plot. An interesting result is found in the F COMP p -F BAND p plot, where the energy flux obtained from BAND and COMP is similar, crossing different categories and groups.
It is even more interesting to see how these parameter relations are affected by the misused models. Based on the categories defined in §2.2, we therefore investigate the following pair parameter relations comparing BAND with COMP: (logF , α), (logF , logE p ), (α, logE p ). For each individual parameter relation, in order to ensure that the majority of spectra are BAND-like, we select three typical bursts (GRB 140206B for the F − α relation; GRB 130306B for the F − E p relation; and GRB 120827 for the α − E p relation), where the vast majority of spectra in these bursts satisfy ∆DIC <-10 (seven out of 10 from GRB 140206B, 13 out of 16 from GRB 130306B, and 30 out of 36 from GRB 120827). We use the following function to fit the data: (Golenetskii et al. 1983), and α − E p relation (e.g., ) are presented in the left, middle, and right panels in Figure 10, respectively. The results of our linear regression analysis for these parameter relations comparing BAND with COMP are reported in Table 2. We find that k Band 1 ∼ 2.87 ± 0.39 is significantly shallower than k CPL 1 ∼ 3.26±0.47, whereas k Band 2 ∼ 1.49±0.15 is clearly steeper than k CPL 2 ∼ 1.22±0.07, and likewise, k Band 3 ∼ −4.62 ± 0.47 is apparently steeper than k CPL 3 ∼ −1.14 ± 0.13.

DISCUSSION
BAND function and COMP are preferred respectively by a given group of GRBs, so a question is raised as to whether these two empirical functions have different physical origins? Or, is COMP just an approximation of BAND as demonstrated in Figure 4 when β << 0? We may find some clues in Figure 6, of which the histogram of β does not form the shape of the unimodal distribution; instead, it displays a clear bimodal structure which peaks at β −2.5 and β −6 respectively. Moreover, the two modes are separated at β −3.5, and each one is almost independently contributed by COMP-preferred spectra or BAND-preferred spectra. This separation is hardly due to that COMP is preferred by noisy data, because first the separation is distinct, and second all the spectra have S > 20 which is Li high enough to ensure data quality 8 . Figure 6 also exhibits the histograms of α and E p , for which the distributions of COMP and BAND preferred GRBs differ from each other, though not as distinguishable as the distribution of β. This statistical result infers that BAND and COMP may have different physical origins. One may propose that there exist two different mechanisms of prompt emission. One produces spectra consists of many power laws, for e.g., synchrotron emission of charged particles accelerated by kinetic shock waves (e.g., Sari et al. 1998); the other produces spectra of a power law with an exponential tail, for e.g., the convolution of blackbody spectra in photospheres, the cut of the highest temperature corresponds to the exponential tail (e.g., Ryde et al. 2010).

CONCLUSIONS
In this paper, we have revisited the catalog of the time-resolved spectrum of the multi-pulse Fermi-GBM bursts defined in . We used two methods to determine the better (or statistically preferred) spectra between two standard empirical spectral functions: BAND and COMP. Firstly, we grouped the spectra into the well-constrained β (BAND-better) and unconstrained β (COMP-better) categories by checking their two-dimensional corner-corner plots of the posteriors for each Bayesian MCMC spectral fit. Secondly, we also separated the spectra into four groups based on the DIC statistics: Group I with ∆DIC <-10 strongly suggests that BAND spectra are statistically preferred; Group II with -10< ∆DIC <-5 indicates that BAND spectra are still statistically preferred, but not as strong as Group I; Group III with -5< ∆DIC <0, indicating COMP spectra are statistically preferred; and Group IV with ∆DIC >0 significantly indicates that the COMP spectra are statistically preferred.
With these categories and groups defended, we therefore compared the spectral properties obtained by both BAND and COMP functions, including their spectral distributions, spectral relations, and spectral evolution.
In the categories defined by identifying well-constrained β and unconstrained β, we found inconsistent peaks of the parameter distributions (both α and E p ) showing between the BAND-better and COMP-better spectra.
These results were also independently confirmed by an analysis based on the DIC statistics. Moreover, such a statistical difference in parameters tends to be weaker when transitioning from Group I (minimum-∆DIC) to Group IV (maximum-∆DIC). The BAND-β distributions show a single peak for all the DIC-statistics-based groups, and the peaks obtained from the Group I and Group IV samples are the same as the harder and softer peaks found in , suggesting that these peaks are more likely dominated by the BAND-like spectra and COMP-like spectra, respectively.
We also discussed the effect of the misused model on the results for each category and group. We found that the apparent deviation from the parameters is found between the BAND-to-COMP cases, while the parameters between the COMP-to-BAND cases all seem to agree.
As a self-consistency test, we also compare the same spectral parameters between the BAND-better and COMPbetter categories and among the DIC-statistic-based groups by investigating the α COMP -α BAND , E COMP p -E BAND p , and F COMP p -F BAND p plotting. The greater dispersion for data points is still found between the BAND-to-COMP cases. We further investigated the F − α, F − E p , and α − E p relations for such the misused case using three example cases. The obtained power-law index (slope) between the misused model (COMP) and the better model (BAND) are significantly different. The index (k COMP ) derived from the "misused" model is shallower (F − E p relation and α − E p relation) and steeper (F -α relation) than that (k BAND ) derived from the better model.
We also discussed the bimodal distribution of β, which indicates that BAND and COMP may have different physical origins.
In conclusion, our analysis suggests that the choice between BAND and COMP spectral model for the GRB spectral analysis should be made with caution. The fit from the misused model deviates from the real spectral shape and then may lead to incorrect physical interpretation. ACKNOWLEDGMENTS I thank the anonymous referee for valuable comments and suggestions. I also thank Felix Ryde, M.G. Dainotti, Gregory Vereshchagin, Remo Ruffini, and ICRANet members for many discussions on GRB physics and phenomena. I particularly thank to Yu Wang for many useful discussions that greatly improved the paper. This research has made use of the High Energy Astrophysics Science Archive Research Center (HEASARC) Online Service at the NASA/Goddard Space Flight Center (GSFC).

APPENDIX
In this Appendix, we also present the spectral parameter evolution and relations (BAND versus COMP) for each individual burst in Figures A1-A2.