Multidimensional analysis of Fermi GBM gamma-ray bursts

TL;DR

This study uses principal component analysis and clustering on Fermi GBM gamma-ray burst data to explore classification, confirming the existence of short and two long-duration classes, but not supporting the intermediate class.

Contribution

The paper applies multivariate statistical techniques to a large gamma-ray burst dataset, providing insights into classification and questioning the intermediate class validity.

Findings

Identifies a short class and two distinct long-duration classes.

Spectral hardness is not a definitive property for class separation.

Fermi GBM spectral variables are complex and not straightforward for classification.

Abstract

The Fermi GBM catalog provides a large database with many measured variables that can be used to explore and verify gamma-ray burst classification results. We have used Principal Component Analysis and statistical clustering techniques to look for clustering in a sample of 801 gamma-ray bursts described by sixteen classification variables. The analysis recovers what appears to be the Short class and two long-duration classes that differ from one another via peak flux, with negligible variations in fluence, duration and spectral hardness. Neither class has properties entirely consistent with the Intermediate GRB class. Spectral hardness has been a critical Intermediate class property. Rather than providing spectral hardness, Fermi GBM provides a range of fitting variables for four different spectral models; it is not intuitive how these variables can be used to support or disprove previous GRB classification results.

Tables

Table 1. Table 1  : Numbers of bursts where the peak energy uncertainties are larger than 1500 keV, 800 keV, and 400 keV.

	$>1500$  $>$ 800 $>$ 400 keV
FlncCompEpeak	65	95	145
FlncBandEpeak	26	49	94
PflxCompEpeak	112	155	244
PflxBandEpeak	54	90	160
PflxSbplBrken	253	329	429
FlncSbplBrken	192	246	308

Table 2. Table 2  : Numbers of bursts for which flux / fluence uncertainties exceed 100% and 50% of the observed values.

	Number of bursts with
	errors exceeding
	100%	50%
T90	147	309
T50	127	329
Fluence	0	1
FluenceBATSE	0	2
Pflux1024	1	5
Pflux64	3	13
Pflux256	1	1
Pflux1024BATSE	0	7
Pflux64BATSE	1	54
Pflux256BATSE	0	9
PflxPlawAmp	0	0
PflxCompAmp	99	456
PflxBandAmp	377	711
PflxSbplAmp	6	11
FlncPlawAmp	0	0
FlncCompAmp	10	92
FlncBandAmp	238	440
FlncSbplAmp	1	1

Table 3. Table 3  : Numbers of bursts where the spectral index uncertainties exceed 3, 1, and 0.5.

	error $>$ 3	error $>$ 1	err $>$ 0.5
flncbandalpha	45	133	278
flncbandbeta	446	589	740
pflxplawindex	0	0	0
pflxcompindex	1	45	248
pflxbandalpha	79	241	516
pflxbandbeta	643	836	1050
pflxsbplindx1	91	208	361
pflxsbplindx2	439	690	942
flncplawindex	0	0	0
flnccompindex	2	12	48
flncsbplindx1	56	112	206
flncsbplindx2	239	436	657

Table 4. Table 4  : The number of GRBs having the best spectral model fit for fluence data (1802 GRBs) and for peak flux data (1792 GRBs).

	PLaw	Comp	Band	SBPL	Total
PLaw	408	492	12	16	928
Comp	22	595	95	64	776
Band	0	4	47	9	60
SBPL	0	5	12	11	28
no fit	10				10
Total	440	1096	166	100	1802

Table 5. Table 5  : Eigenvalues of the PCA from 16 variables.

	Eigenvalue	cumulative %
PC1	7.468	46.7
PC2	3.408	68
PC3	2.458	83.3
PC4	1.522	92.9
PC5	0.411	95.4
PC6	0.195	96.6

Table 6. Table 6  : The coefficients of the largest four PCA eigenvectors.

eigenvectors	1st PC	2nd PC	3rd PC	4th PC
pflxcompindex	0.109	-0.407	0.558	0.591
flnccompindex	0.060	-0.455	0.564	0.579
lgT90	-0.248	0.899	-0.033	0.257
lgT50	-0.288	0.870	0.012	0.292
lgFluence	0.549	0.753	0.163	0.228
lgFluenceBATSE	0.521	0.768	0.126	0.274
lgPflux1024	0.941	0.141	-0.208	-0.103
lgPflux64	0.949	-0.013	-0.089	-0.197
lgPflux256	0.961	0.028	-0.136	-0.181
lgPflux1024BATSE	0.963	0.112	0.031	-0.005
lgPflux64BATSE	0.961	-0.035	0.104	-0.109
lgPflux256BATSE	0.976	0.001	0.070	-0.082
lgPflxCompAmp	0.812	-0.234	-0.299	0.329
lgFlncCompAmp	0.695	-0.459	-0.166	0.384
lgFlncCompEpeak	0.189	0.025	0.886	-0.319
lgPflxCompEpeak	0.187	0.069	0.890	-0.310

Table 7. Table 7  : The communalities obtained using three and four PCs.

	3 PC	4 PC
pflxcompindex	0.49	0.839
flnccompindex	0.53	0.865
lgT90	0.872	0.939
lgT50	0.842	0.927
lgFluence	0.897	0.949
lgFluenceBATSE	0.879	0.955
lgPflux1024	0.951	0.961
lgPflux64	0.91	0.949
lgPflux256	0.944	0.977
lgPflux1024BATSE	0.941	0.941
lgPflux64BATSE	0.937	0.949
lgPflux256BATSE	0.959	0.966
lgPflxCompAmp	0.805	0.913
lgFlncCompAmp	0.723	0.871
lgFlncCompEpeak	0.822	0.925
lgPflxCompEpeak	0.833	0.93

Table 8. Table 8  : Maximal distances between classes (Class 1 = C1, Class 2 = C2, Class 3 = C3).

	C1 v. C2	C1 v. C3	C2 v. C3
pflxcompindex	0.101	0.458	0.417
flnccompindex	0.055	0.542	0.57
lgT90	0.212	0.95	0.87
lgT50	0.226	0.941	0.9
lgFluence	0.382	0.486	0.723
lgFluenceBATSE	0.385	0.658	0.805
lgPflux1024	0.878	0.344	0.623
lgPflux64	0.872	0.772	0.2
lgPflux256	0.875	0.743	0.285
lgPflux1024BATSE	0.802	0.456	0.482
lgPflux64BATSE	0.805	0.843	0.191
lgPflux256BATSE	0.814	0.754	0.179
lgPflxCompAmp	0.692	0.618	0.152
lgFlncCompAmp	0.557	0.786	0.405
lgFlncCompEpeak	0.085	0.799	0.773
lgPflxCompEpeak	0.079	0.765	0.714

Table 9. Table 9  : The probability that the characteristics of two classes are similar. A probability greater than 0.01 indicates that two distributions are not significantly different.

	C1 v. C2	C1 v. C3	C2 v. C3
pflxcompindex	0.041	0	0
flnccompindex	0.608	0	0
lgT90	${10}^{-7}$	0	0
lgT50	${10}^{-7}$	0	0
lgFluence	0	${10}^{-6}$	0
lgFluenceBATSE	0	0	0
lgPflux1024	0	0.001	0
lgPflux64	0	0	0.168
lgPflux256	0	0	0.013
lgPflux1024BATSE	0	${10}^{-5}$	${10}^{-6}$
lgPflux64BATSE	0	0	0.209
lgPflux256BATSE	0	0	0.273
lgPflxCompAmp	0	0	0.464
lgFlncCompAmp	0	0	${10}^{-4}$
lgFlncCompEpeak	0.133	0	0
lgPflxCompEpeak	0.184	0	0

Table 10. Table 10  : Mean values of the sixteen GRB classification variables for each class. Within a row, black numbers indicate a negligible difference, while blue numbers differ significantly. Red numbers differ in a highly significant manner from the other two values.

	C1	C2	C3
pflxcompindex	-0.7164	-0.7264	-0.344
flnccompindex	-0.9761	-1.0082	-0.566
lgT90	1.5343	1.3345	-0.1418
lgT50	1.1082	0.8387	-0.5146
lgFluence	-5.3316	-4.9571	-5.7669
lgFluenceBATSE	-5.5928	-5.2326	-6.184
lgPflux1024	0.6154	1.1047	0.758
lgPflux64	0.8054	1.2518	1.2589
lgPflux256	0.7046	1.1928	1.0924
lgPflux1024BATSE	0.1862	0.7154	0.4185
lgPflux64BATSE	0.4251	0.8946	0.9353
lgPflux256BATSE	0.3026	0.8259	0.7733
lgPflxCompAmp	-1.7188	-1.1881	-1.2939
lgFlncCompAmp	-2.1025	-1.7452	-1.5708
lgFlncCompEpeak	2.3165	2.3008	2.9907
lgPflxCompEpeak	2.3663	2.351	2.9618

Table 11. Table 11  : Contingency table for the 16 and four variable cases with three and four groups.

	1	2	3	4	Total
1	241	68	94	24	427
2	22	115	170	33	340
3	0	0	0	34	34
Total	263	183	264	91	801

Table 12. Table 12  : Contingency table for the 16 and four variable cases with three groups.

	1	2	3	Total
1	391	13	23	427
2	35	237	68	340
3	2	0	32	34
Total	428	250	123	801

Full text

00footnotetext: National University of Public Service, Budapest, Hungary.     corresponding author: horvath.istvan@uni-nke.hu00footnotetext: College of Charleston, Charleston, SC, USA.00footnotetext: Eötvös University, Budapest, Hungary.

Keywords Astronomical databases: miscellaneous – Cosmology: miscellaneous – Cosmology: observations – Gamma-ray burst: general – Gamma-rays: general – Gamma-rays: observations – Gamma-rays: theory – Methods: data analysis – Methods: statistical

1 Introduction

Gamma-ray bursts (GRBs) are the universe’s most energetic electromagnetic events, with observed energies of between $10^{49}$ and $10^{54}$ ergs per burst occurring on timescales ranging from tens of milliseconds to hundreds of seconds (Zhang, 2011; Gehrels and Razzaque, 2013; Kumar and Zhang, 2015; Asano and Mészáros, 2016). Upon correcting for relativistic beaming and large Lorentz factors, these energies are narrowly distributed around $10^{51}$ ergs per burst (Nakar and Piran, 2017).

In the 1980s a phenomenological GRB classification scheme was introduced to subdivide GRBs into Long and Short classes (Mazets et al., 1981; Norris et al., 1984) primarily on the basis of duration. Subsequent observations provided by BATSE (the Burst And Transient Source Experiment) supported this division (e.g. Kouveliotou et al. (1993); Koshut et al. (1996)); the two classes were modeled by overlapping lognormal duration distributions with a delineation occurring at roughly 2$s$. This simple taxonomy unfortunately became ensconced in the literature before more robust statistical and machine learning techniques could be applied, and the general validity of the two-class results overshadowed later concerns about the robustness of the classification approach.

Today it is widely accepted that Short and Long GRB classes have different progenitors, and that they represent dissimilar physical phenomena (Norris et al., 2001; Balázs et al., 2003; Zhang et al., 2009; Lu and Liang, 2010; Lü et al., 2010; Li et al., 2016). The large luminosities and short emission timescales of all GRBs have been explained by models involving black hole formation(Woosley, 2017; Feng et al., 2018; Fernández et al., 2018; Song et al., 2018). The stellar core collapse model used to explain long GRBs occurred on a timescale too long to explain short GRBs, so that models involving compact object mergers in binary systems were developed to explain these (Paczynski, 1986; Usov, 1992; Pérez-Ramírez et al., 2010; Berger, 2014).

The host galaxies and redshift distributions of Short GRBs and Long GRBs differ (Berger, 2014; Levan et al., 2016), with the more luminous Long GRBs being found in star-forming galaxies. Some low-luminosity Long GRBs have been associated with Type Ic supernovae (SN) (Hjorth et al., 2003; Campana et al., 2006; Pian et al., 2006; Blanchard et al., 2016), supporting the idea that the Long GRBs in general are related to deaths of massive stars (Woosley, 1993; Paczyński, 1998; Woosley and Bloom, 2006; Blanchard et al., 2016). For Short GRBs the absence of SN association, the location of these events in metal-poor regions, and their lower luminosities disfavor a massive star origin and point to compact binary mergers (Paczynski, 1986; Usov, 1992; Berger, 2014).

Despite this, the 2$s$ boundary separating Long GRBs and Short GRBs is not robust: the position of the ”boundary“ between Long and Short GRBs depends not only on the methodology used and the size of the database, but also on instrumental characteristics such as the detector threshold and spectral response (e.g. Zhang et al. (2012); Qin et al. (2013)). Furthermore, some GRBs are not amenable to the simple classification scheme. For example, GRB060614 is a long duration burst that shares enough properties with Short GRBs for Zhang et al. (2007) to suggest a compact star merger origin. Similarly, the observed properties of the GRB 090426 indicate that this 1.24$s$ duration burst had a collapsar origin (Lü et al., 2014; Li et al., 2016).

Perhaps because of the lack of a clear demarcation between Short GRBs and Long GRBs, Zhang (2006) proposed classifying GRBs based on their expected behaviors, as predicted from models involving compact objects (Type I) or massive stars (Type II). Subsequently, Zhang et al. (2009) suggested that multi-wavelength properties could be used in conjunction with theoretical models to characterize real physical GRB classes. Li et al. (2016) carried out a multi-wavelength study between duration-defined long-short and the Type I-Type II GRBs. They found several observables to be useful, but not fully reliable, for classification, due to overlap between the populations.

The simple two-class interpretation is further complicated by the statistical existence of a third GRB class. Using multi- and uni-variate statistical analysis techniques, Mukherjee et al. (1998) and Horváth (1998) found evidence for a third GRB class in data from the Third BATSE Catalog (Meegan et al., 1996). Many authors (Hakkila et al., 2000; Balastegui et al., 2001; Rajaniemi and Mähönen, 2002; Horváth, 2002; Hakkila et al., 2003; Borgonovo, 2004; Horváth et al., 2006; Chattopadhyay et al., 2007; Zitouni et al., 2015) have since confirmed the existence of this Intermediate GRB class in the same database using statistical techniques and/or data mining algorithms. The Intermediate class has also been found in the Beppo-SAX (Horváth, 2009) and Swift data (Horváth et al., 2008; Huja et al., 2009; Horváth et al., 2010; Horváth and Tóth, 2016), even though Beppo-SAX had a smaller effective area than BATSE, and Swift works in a different energy range.

The growing number of bursts detected by the Fermi Gamma-ray Burst Monitor (GBM) provides additional data on which GRB classification schemes can be tested. GBM has a spectral energy response that is similar to, but broader, than BATSE, and a surface area that is smaller than that of BATSE. Given the complementary, yet different, characteristics of the Fermi GBM instrument to BATSE, Swift, and Beppo-SAX, the time has come to apply statistical clustering techniques to explore GRB classification using this instrument.

The paper is organized as follows. Section 2 discusses the properties of the Fermi GBM catalog, Section 3 defines thirty six potential classification parameters and their structures, and Section 4 describes the classification process using sixteen GRB parameters following the elimination of sixteen duplicative parameters. Section 5 discusses the results and Section 6 provides the paper’s conclusions.

2 Classification Variables from the Fermi GBM Catalog

On April 18, 2017, the Fermi GBM Catalog contained 2060 GRBs.111https://heasarc.gsfc.nasa.gov/W3Browse/fermi/fermigbrst.html The most recent GRB at that time was GRB170416583, although the last GRB for which spectral fits were available was GRB170131969 (the 2011th GRB). More than three hundred variables that might be used in GRB classification are published for each GRB, including negative and positive uncertainties. Upon excluding the uncertainties, the number of available variables is less than two hundred. Most of these are spectral fit parameters, including fitting parameters for each of four different spectral models: Band, Comptonized (Comp), Power Law, and Smoothly Broken Power Law (SBPL); for details see Narayana Bhat et al. (2016). All four fits are applied to both the set of photons observed in the one second-peak (the 1024 ms peak flux) and the set of all photons in the burst (the fluence).

Published analyses have used a variety of different variables in their classification approaches. Some analyses use only the duration information (Horváth, 1998; Balastegui et al., 2001; Rajaniemi and Mähönen, 2002; Horváth, 2002; Horváth et al., 2008; Huja et al., 2009; Horváth, 2009; Zitouni et al., 2015; Tarnopolski, 2015a, b; Horváth and Tóth, 2016; Tarnopolski, 2016; Kulkarni and Desai, 2017), others primarily use the duration-hardness plane (Horváth et al., 2004, 2006; Veres et al., 2010; Horváth et al., 2010; Koen and Bere, 2012; Qin et al., 2013; Tsutsui and Shigeyama, 2014; Shahmoradi and Nemiroff, 2015; Řípa and Mészáros, 2016; Yang et al., 2016; Zhang et al., 2016), while still others use more than two variables (Mukherjee et al., 1998; Hakkila et al., 2003; Chattopadhyay et al., 2007; Kann et al., 2011; Lü et al., 2014; Li et al., 2016; Modak et al., 2017; Chattopadhyay and Maitra, 2017).

For our analysis we have chosen to use thirty six potentially useful classification variables, even if they overlap in content, with the assumption that any deemed to be nonsense parameters can be removed later. There are two duration measures (T90 and T50), two fluence measures (one based on Fermi GBM energy channels and the other based on BATSE energy channels), six peak flux measures (on the 64-ms, 256-ms, and 1024-ms timescales measured in both the Fermi GBM and BATSE energy channels), eight Band spectral fit parameters ($\alpha$, $\beta$, $E_{\rm peak}$, and the fit amplitude obtained from both the peak flux and fluence spectra), eight broken power law spectral fit parameters (low-energy index, high-energy index, break energy, and fit amplitude obtained from both the peak flux and fluence spectra), four single power law spectral fit parameters (power law index and the fit amplitude obtained from both the peak flux and fluence spectra), and six Compton spectral fit parameters (spectral index, peak energy, and fit amplitude obtained from both the peak flux and fluence spectra). Because most of these parameters span large dynamical ranges, we use the base 10 logarithmic measures of these variables (lg) instead of the measures themselves. There is no need to take the logarithmic values of the spectral indices since these have already been obtained as power law functions.

3 Statistical Clustering of GBM Data

3.1 Errors and correlations of pre-selected classification variables

Large numbers of objects and many classification variables do not guarantee that a GRB classification will be successful. Because GRBs are observed in low signal-to-noise regimes where measurements are difficult to make with accuracy, large uncertainties often accompany potential classification variables. Furthermore, the information contained in each classification variable is not always independent; many variables have content that overlaps with that found in others (e.g., T90 and T50).

Before applying classification tools to this GRB sample, we reduce the scatter in the data by removing GRBs having large measurement uncertainties; the numbers of such GRBs are found in Tables  1 -  3. Among the 2060 GRBs in our sample, 1729 have usable spectral fits of one form or another and 1487 have values measured for all 36 variables.

We have calculated the correlation matrix in order to determine the amount of information overlap found in the 36 classification variables. The correlation matrix is found with the built-in R-environment (R Core Team, 2015) cor() function using the available Pearson, Kendall, and Spearman methods. Seriation is performed on each of these matrices with the package seriation (Hahsler et al., 2008) using the function seriate() and the methods PCA and PCA_angle. The results obtained from the different techniques are very similar to one another, and the correlation matrix constructed with the Spearman method and seriated with the PCA_angle method is shown in Figure  1.

Clusters of variables are observed along the main diagonal where content overlap is significant: The two fluences are related to T90 and T50 and to the peak fluxes and the spectral amplitudes (the large 14x14 square in the lower right side of the matrix). The content overlap between fluence, peak flux, and duration make intuitive sense as fluence is the time-integrated flux, measured over the burst’s duration (e.g.,Hakkila et al. (2003)). Spectral amplitudes are also types of peak fluxes, but that are binned spectrally rather than temporally. Low energy spectral indices form another cluster in the matrix, as these indices all measure low-energy burst behaviors. The power law spectral indices overlap with break energies rather than with the other spectral indices because power-law indices identify single-component spectral models that increase through measured energies rather than decreasing at high energies (as more complex spectral models do). Finally, on the upper left side of the matrix, the four high energy spectral indices form a distinct group.

3.2 Spectral fits for GRBs in the Fermi GBM catalog

The Fermi GBM catalog contains information about which of the four spectral fits is best for each burst. Spectral fits have been obtained at the time of the peak flux for 1792 of the 2060 GRBs, and using the fluence for 1802 of the 2060 GRBs. Among the peak flux spectral fits, the power-law function provides the best fit for 928 GRBs, the Compton model for 776 GRBs, the Band function for 60 GRBs, and the SBPL function for the remaining 28 GRBs. Among the fluence spectral fits, the power law model provides the best fit for 440 GRBs, the Compton model for 1096 GRBs, the Band model for 166 GRBs, and the SBPL model for the remaining 100 GRBs.

Table  4 summarizes the consistency of different spectral fit models obtained from fluence and peak flux spectral data, with fits to peak flux data listed in rows and fits to fluence data listed in columns. The peak fluxes of 10 bursts provide insufficient photon counts to generate any spectral fit and are shown in the first column; the fluences of these bursts have only marginally better statistics resulting in simple power law fits.

3.3 Revisiting the effectiveness of the spectral fit variables

The spectral variables described in Section 3.1 have been included in our tabulations regardless of the quality of these fits. To more accurately evaluate the potential effectiveness of the variables, we examine the reduced $\chi^{2}$s of each fit, as published in the Fermi GBM catalog.

Figures  2 and  3 demonstrate the reduced $\chi^{2}$ and associated degree of freedom (DOF) distributions for the power law spectral fits obtained from peak flux measurements. We note that the DOF distribution is identical for all eight spectral fits.

Most of the DOF values exceed 430, 96% are larger than 300, and only one is less than 200. For DOF = 400, a reduced $\chi^{2}$ value of 1.172 indicates a 99% significance (a value of 1.232 corresponds to 99.9%). For DOF = 300, the reduced $\chi^{2}$ value of 1.199 indicates a 99% significance (1.271 corresponds to 99.9%), and for DOF = 200 the reduced $\chi^{2}$ value of 1.247 indicates a 99% significance (1.337 corresponds to 99.9%). This means that, for the DOF values found in this distribution of Fermi GBM bursts, reduced $\chi^{2}$ values of 1.25 indicate poor spectral fits and values of 1.35 indicate very poor spectral fits. The reduced $\chi^{2}$ values of the different spectral fits range from being mostly poor ($>1.25$) in the case of the power law model fits (see Figure  2) to being generally good ($<1.25$) in the case of the Compton model fits (see Figure  4).

The signal-to-noise ratio ($S/N$ = value/error) is another quality indicator that should be examined when considering the value of potential classification variables. The effects of $S/N$ on spectral fit measurements can be seen in Figure  5, which shows the $S/N$ distribution of Band spectral fit amplitudes measured from fluence spectra. Several hundred GRBs in this distribution have ($S/N<1$), indicating that their fit amplitudes are less than their fit amplitude uncertainties. Another few hundred GRBs have signal-to-noise ratios smaller than two or three.

Very large uncertainties are also present in some of the high energy spectral indices measured for Fermi GBM bursts. For example, Figure  6 demonstrates the uncertainties in the Band $\beta$ spectral index obtained from peak flux measurements. Since $\beta$ typically has values in a narrow range ($-3\leq\beta\leq-2$), $\beta$ values having uncertainties ($\sigma_{\beta}>1$) are not likely to be useful for classification. Unfortunately, many Fermi GBM bursts fitted by the Band model fall into this range; several hundred have $\sigma_{\beta}>10$, and many have $\beta$ values that exceed 100 and even 1000.

Due to the uncertainties in measuring high energy power-law spectral indices, we have omitted the Band and SBPL fits (for all four cases, both fluence and peak flux fits) in this analysis.

We similarly analyze the $\chi^{2}$ distributions of the Compton and power law spectral models to determine whether the variables obtained from these fits have potential value as classification variables. Three hundred eighty one power law model and 60 Compton model $\chi^{2}$ values of the 1792 peak flux spectra are deemed unacceptable, as are 703 power law and 180 Compton model $\chi^{2}$ values. If we want to use for classification the acceptable spectral model variables from both peak flux and fluence spectral fits, then only the Compton model produces enough variables for a sufficiently large sample. This reduces the sample size by roughly half. We also exclude GRBs where the uncertainty in any variable exceeds 50% of that variable’s value. This excludes, for example, 184 GRBs on the basis of T90 and 288 GRBs on the basis of peak flux amplitude. Only 803 GRBs satisfy all of our requirements for a well-measured classification variable database. Finally, two more bursts have been excluded using the HDoutliers() function of the R-package HDoutliers (Fraley, 2016).

The 801 remaining GRBs are chosen to be the sample for the analysis described in Section 4. The list of these GRBs and the values of the variables can be found at http://itl7.elte.hu/ hoi/grb/f1704C801v16.csv

4 Cluster analysis with 16 variables

We have reduced the number of GRB classification variables by limiting our sample to the following sixteen well-measured characteristics of bursts in the Fermi GBM Catalog (note that all denoted by ’lg’ represent logarithmic values): lgT90 (T90), lgT50 (T50), lgfluence (fluence), lgPflux64 (64 ms peak flux), lgPflux256 (256 ms peak flux), lgPflux1024 (1024 ms peak flux), lgfluenceBATSE (fluence in the BATSE energy channels), lgPflux64BATSE (64 ms peak flux in the BATSE energy channels), lgPflux256BATSE, (256 ms peak flux in the BATSE energy channels), lgPflux1024BATSE (1024 ms peak flux in the BATSE energy channels), lgflncCompAmp (Compton amplitude from fluence spectral fit), lgflncCompEpeak (Compton Epeak from fluence spectral fit), flncCompIndex (Compton power law index from fluence spectral fit), lgpflxCompAmp (Compton amplitude from peak flux spectral fit), lgpPflxCompEpeak (Compton Epeak from peak flux spectral fit), and pflxCompIndex (Compton power law index from fluence spectral fit). Using the techniques described in Section 3.1, we obtain the correlation matrix using the built-in cor() function using the Spearman method and seriate the results with the seriate function of the FactoMineR package (Husson et al., 2016) via the PCA_angle method. The result can be seen in Figure  7.

Four distinct blocks of correlated variables are apparent in the matrix:

block1: lgPflxCompEpeak and lgflncCompEpeak

block2: pflxCompIndex and flncCompIndex

block3: lgPflux64, lgPflux256, lgPflux1024, lgPflux64BATSE, lgPflux256BATSE, lgPflux1024BATSE, lgflncCompAmp, lgpflxCompAmp

block4: lgT90, lgT50, lgfluence, lgfluenceBATSE

We have carried out PCA (principal component analysis) for the Fermi GBM bursts using these 16 variables. Table  5 identifies the six largest PCs, while Table  6 shows the coefficients of the four largest PCA eigenvectors, and Table  7 identifies the communalities of these PCs. To determine whether or not PCs are important, we use the criterion that significant PCs should be larger than unity and yield cumulative percentages larger than $80\%$ (see Jolliffe (1986)). Table  5 demonstrates that the top five PCs contain a combined 95.2% of the non-overlapping information contained in the sixteen variables, and the top six PCs contain 96.6% of that information.

Table  7 shows that the communalities obtained from three PCs are almost as robust as those obtained from four PCs (although the communalities of the two Compton indices are reduced to only 50% with three PCs), so we choose the three largest PCs for our analysis. The optimal number of classes contained within these PCs can be found with guidance from the Bayesian Information Criterion (BIC) (Schwarz, 1978). Figure  8 shows that the maximum BIC value (obtained from the mclust() function of the mclust (Charrad et al., 2014) R-package) suggests that the three PCs from the Fermi GBM 16-variable data optimally describe three GRB classes.

Cluster analysis is performed on the 801 GRB sample using the mclust() function with the aforementioned three PCs serving as classification variables. The function returns three classes: Class 1 (containing 427 GRBs), Class 2 (containing 340 GRBs) and Class 3 (containing 34 GRBs). The values returned by the mclust() function also contain probabilities $p1_{i}$,$p2_{i}$, and $p3_{i}$ that the $i$th GRBs belongs to Class 1, Class 2 or Class 3, respectively. We sum these probabilities to get $p1=400.1$ for Class 1, $p2=365.5$ for Class 2 and $p3=35.4$ for Class 3. These numbers differ from 427, 340 and 34 because each GRB has been placed into the group for which its probability is largest, and the non-integer counts of the group elements created in this way does not necessarily yield the same total as the some of the integer count probabilities.

We have calculated the $\log(T90)$ distributions of all three GRB classes by assuming that each GRB belongs to the class associated with it largest cluster probability (Max[$p1_{i}$, $p2_{i}$, or $p3_{i}$]). The three $\log(T90)$ distributions are shown in Figure  9, after having been normalized by the factors 427 for Class 1, 340 for Class 2 and 34 for Class 3. In addition to duration, we have also calculated the class distributions for the other fifteen variables. The fluence distributions of the GRB classes can be seen in Figure  10. One can easily see from these figures that a short duration GRB class (Class 3) is present, consistent with results obtained from other GRB experiments using a variety of classification techniques.

We have calculated the maximal differences between the three class distributions for each of the 16 observed parameters using the Kolmogorov-Smirnov statistic (implemented in the R environment by the function ks.test()); the results are listed in Table  8 and are shown in Figure  11. For the logT90 distributions, the maximal distance between Class 3 and Class 1 is 0.95 and between Class 3 and Class 2 is 0.87, indicating that the Class 3 durations are discernibly shorter than those of the other two classes. In contrast, the maximal distance between Class 1 and Class 2 is only 0.212.

The KS probabilities associated with these maximal differences are shown in Table  9. Using a significance requirement of at least 99.9%, probabilities larger than 0.001 indicate that there is no significant difference between two distributions. A near-zero probability therefore indicates a very significant difference (In our calculation zero means less than $10^{-7}$). Similar to Figure  11, which plots the numbers from Table  8, Figure  12 indicates the same D values colored according to their probabilities. Red indicates no significant difference ($p>0.1\%$), black indicates a very significant difference, and blue means that the two compared distributions may differ significantly.

Table  10 contains the mean values of each of the sixteen variables for all three GRB classes. Red indicates that a value differs significantly from the other two class values in the same row. Black numbers in a row do not differ significantly from one another. Blue numbers in the same row may be significantly different. Figure  13 shows the mean and the standard deviation of the sixteen GRB classification variables for each class. Figure  14 and Figure  15 show some 2-Dimensional distributions of our variables and also their correlation coefficients.

5 Robustness

It is important to check the robustness of our analysis. We do this by first excluding the variables lgFluenceBATSE, lgPflux1024BATSE, lgPflux64BATSE and lgPflux256BATSE. Second, we keep only the variables lgT90, lgFluence, lgPflux256 and flnccompindex which are independent of one another.

5.1 Robustness test with twelve variables

Excluding the BATSE-related variables (lgFluenceBATSE, lgPflux1024BATSE, lgPflux64BATSE, lgPflux256BATSE) the number of variables reduces to twelve. We performed cluster analysis with the same 801 GRBs using the mclust() function with these twelve variables. The BIC results can be seen in Figure 16.

The BIC function reaches its maximum at four assumed groups. Thus, using the twelve Fermi parameters the mclust method prefers four groups rather than three. One can compare this result with the previous result by calculating the contingency table (Table 11) for cases of sixteen and twelve variables with three and four GRB groups. The table shows the stability of the group members. For the first group, 56% (241/427) remain in the same group. For the short bursts, all 34 GRBs remain in the same group (group4 in the twelve variable case). For the second group, 84% ((115+170)/340) remain in the new group2 and group3.

5.2 Robustness test with four variables

Keeping only four variables (lgT90, lgFluence, lgPflux256, flnccompindex) we performed cluster analysis with the same 801 GRBs using the mclust() function. The BIC function reaches its maximum at three assumed groups (see Figure 17).

One can compare this result with the previous result by calculating the contingency table (Table 12) for the sixteen and four variable cases with three assumed groups. The table shows the stabilities of the group members. For the first group, 92% (391/427) remain in the same group. For the second and the third groups, 70% (237/340) and 94% (32/34) of the GRBs do not change their group membership.

6 Discussion

Soon after the discovery of GRBs, it was suggested that they could be separated into two duration classes (Mazets et al., 1981; Norris et al., 1984). This hypothesis was later supported by BATSE data (Kouveliotou et al., 1993; McBreen et al., 1994; Koshut et al., 1996), and resulted in the identification of ‘Long’ and ‘Short’ GRB classes separated at around $T90=2s$. An additional ‘Intermediate’ class (with durations of roughly $2sT90<15s$) was subsequently found from analyses of the Third BATSE Catalog data Mukherjee et al. (1998); Horváth (1998); numerous authors (Hakkila et al., 2000; Balastegui et al., 2001; Rajaniemi and Mähönen, 2002; Horváth, 2002; Hakkila et al., 2003; Borgonovo, 2004; Horváth et al., 2006; Chattopadhyay et al., 2007; Zitouni et al., 2015) have verified this result using BATSE GRBs. The Intermediate class has also been found in Beppo-SAX (Horváth, 2009) and Swift data (Horváth et al., 2008; Huja et al., 2009; Horváth et al., 2010; Horváth and Tóth, 2016), as well as in a preliminary analysis of Fermi data Horvath et al. (2012)s.

The Intermediate GRB class is not as clearly delineated from the other classes as the Short class is; this is also true from recent analysis of the Swift data. Zhang and Choi (2008) analyzed 95 GRBs with measured redshifts and found that Swift bursts have a wider dynamic range in duration than pre-Swift and BATSE bursts. Koen and Bere (2012) subsequently analyzed the Swift BAT data and concluded that two classes sufficiently describe the spectral hardness distribution, whereas three components are needed to characterize the duration distribution. The Intermediate class identified by Koen and Bere (2012) has durations of around 3-20 seconds, which is in a good agreement with (Horváth and Tóth, 2016) who find the Intermediate class durations to be in the 4 to 30 second range.

In addition to Swift, instrumental effects might also be responsible for affecting Fermi classification results. Qin et al. (2013) analyzed the data of 315 Fermi GRBs, studying the dependence of the duration distribution on energy and on various instrumental and selection effects. They have suggested that the true durations of a GRB could be much longer than what is observed. They also suggested that the observed bimodal duration distribution might be due to an instrumental selection effect.

Analysis of data from a variety of orbital high-energy satellites continues to find evidence for three GRB classes. In a series of papers, Tsutsui and his coworkers have used data from several orbital instruments, as well as x-ray and optical afterglow data, to study GRB classes (Tsutsui et al., 2009, 2013; Tsutsui and Shigeyama, 2014) and have found a third group with durations of approximately five seconds. Zitouni et al. (2015) has analyzed the CGRO/BATSE and Swift/BAT GRB data to find a very similar class structure to (Horváth, 2002).

Although most rigorous GRB classification studies find three classes in the data, there have been exceptions. In one of his recent publications, Tarnopolski (2015b) proposes that the division between short and long bursts is at 3.4 seconds rather than at two seconds. In Tarnopolski (2015a) he analyzed the Fermi BAT duration data of 1566 GRBs. Although he found a third component in the distribution, the significance was not convincing. This may be due to methodology: he tested binned data with a $\chi^{2}$ fit rather than using a maximum likehood method with unbinned data. Tarnopolski also suggested that the log-normal fit may not be an adequate model for the duration distribution, which is an interesting suggestion because previous investigators have assumed that the underlying distribution is lognormal.

In this manuscript we have analyzed 801 Fermi GRBs observed by GBM using sixteen classification variables: two durations (T90 and T50), three Compton spectral parameters (amplitude, peak energy, and spectral index) for both fluence and peak flux fits, two fluences, and six peak fluxes. We find that content overlap in these variables can be reduced to a three-variable problem. These three main components are essentially the spectral (Comp) index, an amalgam of the peak fluxes, and a combination of fluences and durations (see Table  6). These variables can be used to identify three clusters (classes) of GRBs. As Table  10 and Figure  13 show, Class 3 (34 bursts) has short durations while the other two have long durations; Class 2 (340 bursts) is slightly shorter than Class 1 (427 bursts). Class 2 is, however, brighter as measured by fluence than Class 1, which is in turn brighter than Class 3. The peak fluxes of Class 2 and 3 are generally similar, and both of these are brighter than the peak fluxes of Class 1. There are no clear spectral differences between Class 1 and Class 2 bursts, although the Compton indices and peak energies of Class 3 are larger than those of the other two classes.

These three GRB classes are partially but not entirely recognizable when compared to those obtained from previous analyses. Class 3 strongly resembles the Short class found in the 2- and 3-class cluster analyses of other GRB experiments. Like Short bursts, Class 3 GRBs are are less common, shorter, spectrally harder, and fainter than other burst classes. Classes 1 and 2 are not recognizable as the Intermediate and Long burst classes found in other 3-class cluster analyses. Classes 1 and 2 are both equally common and of similar spectral hardness. Class 2 is shorter yet brighter than Class 1. In contrast, Intermediate bursts are less common, shorter, spectrally softer, and fainter than Long GRBs.

We have identified several reasons why our analysis might not have delineated the Intermediate GRB class from the Long one. First, as indicated earlier in this section, the Intermediate class is not as well defined as the Long class in BATSE, Swift, and Beppo-SAX data; this may make it harder to find in the data collected by a different experiment. Second, Fermi GBM’s smaller surface area makes it less sensitive to detecting faint GRBs, such as those belonging to the Intermediate class (de Ugarte Postigo et al., 2011): the relatively small number of faint but clearly identified Short GRBs observed by Fermi GBM supports this statement. Third, and most important, the spectral information presented in the Fermi GBM catalog is in a form that requires more interpretation than the hardness ratios extracted from the data of other experiments.

Spectral fit parameters are published in the Fermi GBM catalog rather than simpler traditional data products such as hardness ratios. This is because burst spectral fitting is a complex process, and the Fermi GBM science team has chosen to publish the spectral fits so that burst intensity, localization, and spectral information have been correctly deconvolved from the detector response. This approach serves the uninitiated user, while also taking advantage of Fermi GBM’s ability to time- and energy-tag every photon. In contrast, other GRB experiments generally publish bulk spectral characteristics in the form of hardness ratios, which are simpler but not as intuitive or as potentially useful for theoretical modeling. Although the Fermi GBM approach has many advantages, it is disadvantageous when classifying GRBs because spectral fitting requires photons to be binned in many (16) rather than few (4) energy channels, and there are generally too few photons in a Fermi GBM burst for the spectrum to be fitted accurately. Since the underlying physical GRB spectral model is not known, approximations are obtained in the form of the four ‘standard’ spectral models. None of these models is accepted as representing the underlying burst physics, but each works well for burst spectra having certain characteristics at the cost of being less optimum for other types of spectra. For example, the Band and SBPL models provide better spectral fits when a significant amount of high energy flux is present, whereas the PL and Compton models provide better fits when high energy flux is lacking (Goldstein et al., 2012; Gruber et al., 2014).

The choice of a ‘best’ spectral model differs from burst to burst, depending on factors such as the true energy distribution of photons, the burst brightness, the burst redshift, and the detector response. Furthermore, GRB spectra evolve, and the time interval during which an evolving spectrum has been observed can impact the choice of which spectral model is optimal for fitting a burst: rapidly evolving GRB spectra are usually spectrally harder near the time of the peak flux, which explains why peak flux spectra might differ from fluence spectra. It is entirely likely and natural to think that a GRB’s optimal spectral model might change throughout its prompt emission.

However, in order to reduce measurement uncertainties in the difficult-to-measure spectral variables, we have limited our classification database to those GRBs for which the Compton model provides the ‘best’ spectral fits from both peak fluxes and fluences. In doing so we have likely biased our sample to bursts having few spectral differences that might not be representative of the larger distribution. Without having wide-ranging and well-defined burst spectral data, it is easy to see how our analysis might have had trouble clearly delineating the Intermediate class from the Long one, even as it finds evidence supportive of three GRB classes.

7 Summary and conclusion

We classify GRBs using cluster analysis and the attributes of the published Fermi GBM burst catalog. After initially selecting thirty-six potential classification variables along with 2016 prospective bursts, we show that measurement uncertainties limit our classification sample to 810 GRBs and only sixteen classification variables. Principal Component Analysis reduces the number of non-overlapping classification variables to three, constructed primarily of peak fluxes, a spectral index, and fluences coupled with durations.

Cluster analysis identifies three optimal GRB classes. The first appears to be the well-known Short GRB class, while the other two are types of Long classes. The peak flux distributions of these two Long classes differ, while their hardness distributions do not (duration and fluence distributions are slightly different). Neither one of these classes appears to be the previously-identified Intermediate class. We attribute this discordant result to weakly-delineated Intermediate class characteristics obtained from data in other GRB catalogs, coupled with the Fermi GBM catalog’s use of various model-dependent spectral fitting parameters as opposed to standard hardness ratios.

We are currently pursuing a more detailed analysis involving the various spectral models and associated fitting parameters and their effects on GRB classification. We hope to determine (1) whether or not systematic differences in the published GRB spectral fitting parameters can be overcome to yield understandable classification results, and (2) the extent to which the Intermediate class is present in the Fermi GBM catalog.

This research was supported by NASA EPSCoR grant NNX13AD28A and the Hungarian TIP and TKP grants and OTKA K131653 grant.