A Tight Three-parameter Correlation and Related Classification on Gamma-Ray Bursts

Gamma-ray bursts (GRBs) are the most violent explosions from the deep universe and have been studied extensively ever since their discoveries in the late 1960s (Klebesadel et al. 1973). They are generally considered to be composed of short GRBs, with durations shorter than 2 s, and long GRBs, which last longer (Kouveliotou et al. 1993). The long GRBs are most likely powered by the collapse of a (rapidly rotating) massive star often associated with a bright supernova (MacFadyen & Woosley 1999; Woosley & Bloom 2006), while the short GRBs are expected to be from coalescence of compact object binaries involving at least one neutron star (Eichler et al. 1989; Paczynski 1991; Narayan et al. 1992).

There seem to be obvious differences in emission characteristics between long and short GRBs, for example, short GRBs have statistically larger hardness ratios than the long GRBs (e.g., Kulkarni et al. 2000). In addition, the host galaxies of long GRBs also have different features from those of short GRBs. Long GRBs whose progenitors are young massive stars occur in star-forming galaxies (Bloom et al. 2002; Fruchter et al. 2006; Li et al. 2011) while the host galaxies of short GRBs include both early-type and star-forming galaxies (Fox et al. 2005; Berger et al. 2005; Bloom et al. 2006; Fong et al. 2011; Fong & Berger 2013). The host galaxies of long GRBs also have more concentrated and smaller age distributions, smaller stellar mass (Savaglio et al. 2009; Leibler & Berger 2010), and lower metal abundance (Kobulnicky & Kewley 2004; Tremonti et al. 2004). The two subclasses of GRBs also have different offset (distance from the center of the host galaxy) distribution (Fong & Berger 2013). Zhang et al. (2017) have found a correlation between the luminosity of short GRB afterglows and the offset from host galaxies, while for long GRBs, neither prompt emission nor afterglow satisfied a similar correlation with its offset.

In spite of the many differences proposed between long and short GRBs, there is not a clear-cut observational standard to strictly distinguish the two types of bursts. In addition, it is not very certain whether GRBs can only be divided into two subclasses or there should be more subclasses. Over the years, researchers have proposed different classification scenarios, such as introducing a medium duration subclass between long and short GRBs based on T₉₀ (Horváth 1998; Mukherjee et al. 1998; Hakkila et al. 2000; Horváth et al. 2010) or classifying GRBs by defining and applying variable parameters (Lü et al. 2010). Li et al. (2016) conducted a comprehensive study of the comparison between long and short GRBs by analyzing 14 relevant parameters in detail and came up with three best parameters for the classification purpose (T₉₀, f_eff, F_light). However, no single parameter alone is yet good enough to accurately determine the category for a particular burst. Therefore, multiple-parameter criteria for physical classification seem to be unavoidable.

The development of machine learning and data mining technology provide us with new research tools. Traditional machine learning methods such as principle components analysis (PCA) and kth nearest neighbor, etc., can be used to perform spectral line certification in the spectrum and eliminate background radiation (Škoda et al. 2016; Marchetti et al. 2017). Ukwatta et al. (2016) estimated the redshift of GRBs by using the random forest algorithm.

In terms of GRBs classification, there have also been some related applications. With new methods adopted, however, some new and different results or conclusions about GRB classification have been proposed. Bhave et al. (2017) analyzed two-dimensional clustering of GRBs using duration and hardness and obtained different classification conclusions of either two or three components for different detectors or data sets. Svinkin et al. (2019) analyzed GRBs detected by Konus-Wind with similar methods and also suggested a three-component classification. Chattopadhyay & Maitra (2017, 2018) divided GRBs into five categories by Gaussian (and t-) mixture-model-based cluster analysis, and Horváth et al. (2019) proposed a three-subclass GRBs classification scheme. Meanwhile, some researches (Tarnopolski 2019a, 2019b) imply that the third class is unlikely to be a real phenomenon. Briefly, either the number of categories is different, or the characteristics of the subclasses are different. Thus the issues of GRB classification are far from being solved.

In this work, we analyze the sample of 2989 GRBs detected by Fermi GBM detectors The parameters included in our analysis are duration T₉₀ (during which 90% of the burst fluence was accumulated, in units of seconds), fluence f (flux integrated over the full burst duration in 10–1000 keV energy band, in units of erg cm⁻²), and peak flux F (the peak flux in 10–1000 keV energy band on the timescale derived by Bayesian Block method (Scargle et al. 2013), in units of photon cm⁻² s⁻¹), which are publicly available from the Fermi GBM burst catalog (von Kienlin et al. 2020). We consider the data as a 2989 × 3 matrix X, which has three columns consisting of log T₉₀, $\mathrm{log}f$ and $\mathrm{log}F$, respectively. We will analyze the three-dimensional distribution of this data set.

All the algorithms used in this work are mature and complete tools in many fields. We will not go into details about mathematical theory, and just briefly introduce the characteristics of these algorithms related to the problems studied in this work.

2.1. Principle Component Analysis

2.2. Linear Discriminant Analysis

2.3. Gaussian Mixture Model

3.1. A Tight Three-parameter Correlation

As shown in the left panel of Figure 1, the three-dimensional distribution of the total fluence (f), peak flux (F), and duration T₉₀ concentrates very tightly in a certain plane of logarithmic scales. The best-fit function for this plane could be easily determined with a multivariate regression analysis as in:

$$\begin{eqnarray}&&\mathrm{log}\,f=0.75\,\mathrm{log}\,{T}_{90}+0.92\mathrm{log}F-7.14.\end{eqnarray} \tag{ 1 }$$

It is a very tight global correlation for all GRBs in the sample with a correlation coefficient of R = 0.92. A corresponding two-dimensional projection as shown in the right panel of Figure 1 appears be more intuitive.

Zoom In Zoom Out Reset image size

PCA is applied to the sample to compute the principal components and then transform the data into a new coordinate system that takes the principle components as the new axes. Accordingly, we get the principle components as in:

$$\begin{eqnarray}{P}=\left(\begin{array}{ccc}\,0.75 & -0.44 & -0.48\,\\ \,0.65 & 0.47 & 0.59\,\\ \,0.04 & 0.76 & -0.64\end{array}\right).\end{eqnarray} \tag{ 2 }$$

The initial sample variables X can be transformed into new coordinate system ${{X}}^{{\prime} }$ by:

$$\begin{eqnarray}&&{X}^{{\prime} }=(X-C)\cdot P\end{eqnarray} \tag{ 3 }$$

where C = (1.12, − 5.55, 0.81) is the center (judged by mean values) of three-dimensional distribution of initial sample variables X. Therefore, we can express the new three principle components as:

$$\begin{eqnarray}\begin{array}{ccl}{{X}}^{{\prime} } & = & \left({P}_{1}\ {P}_{2}\ {P}_{3}\right)\\ & = & {\left(\begin{array}{c}0.75\,\mathrm{log}\,{{\mathsf{T}}}_{90}+0.65\mathrm{log}\,f+0.04\mathrm{log}F+2.76\\ -0.44\,\mathrm{log}\,{{\mathsf{T}}}_{90}+0.47\mathrm{log}\,f+0.76\mathrm{log}F+2.47\\ -0.48\,\mathrm{log}\,{{\mathsf{T}}}_{90}+0.59\mathrm{log}\,f-0.64\mathrm{log}F+4.36\end{array}\right)}^{{\mathsf{T}}}\end{array}.\end{eqnarray} \tag{ 4 }$$

Three columns of the matrix ${{X}}^{{\prime} }$ correspond to three principle components P₁, P₂, and P₃, respectively, which have been reverse sorted according to the variance. We can also evaluate the variance ratios (Vr) of the three components as in:

$$\begin{eqnarray}&&{Vr}=\left({Vr}1\quad {Vr}2\quad {Vr}3\right)=\left(0.70\quad 0.27\quad 0.02\right).\end{eqnarray} \tag{ 5 }$$

The variance ratio of the third principle component is only 0.02, which is almost negligible in comparison with that of the other two principle components. So the three-dimensional data can be reduced to two dimensions by dropping the third component. The updated distribution of our sample data in the new coordinate system ${{X}}^{{\prime} }$ is shown in Figure 2.

Zoom In Zoom Out Reset image size

From another perspective, the dropped component is not useless. Instead, it carries important consistency information of the data sample. In other words, the third principle component indicates a correlation between data parameters of $\mathrm{log}{T}_{90},\mathrm{log}f$, and $\mathrm{log}F$. Equation (3) indicates that the data sample should be centered at 0 in the new coordinate system. So, from the third principle component, the mean value of data is:

$$\begin{eqnarray}&&-0.48\,\mathrm{log}\,{T}_{90}+0.59\mathrm{log}\,f-0.64\mathrm{log}F+4.36=0.\end{eqnarray} \tag{ 6 }$$

This is essentially a constraint on the data sample, which is a plane in the three-parameter space. Through equivalent transformation of the formula, we can see that Equations (6) and (1) are almost the same.

To sum up, even though the dropped component in PCA rarely contains the difference information, it does contain the most consistent information of the sample. Later in this work, we will carry out analyses on the classification of GRBs based on the remaining components, i.e., the first and second principle components, which contain most of the difference information.

3.2. Improved Classification Criteria

Based on the result of PCA, we consider first and second principle components as a new data matrix:

$$\begin{eqnarray}\begin{array}{ccc}{X}_{2} & = & \left({P}_{1}\ {P}_{2}\right)\\ & = & {\left(\begin{array}{c}0.75\,\mathrm{log}\,{{\mathsf{T}}}_{90}+0.65\mathrm{log}\,f+0.04\mathrm{log}F+2.76\\ -0.44\,\mathrm{log}\,{{\mathsf{T}}}_{90}+0.47\mathrm{log}\,f+0.76\mathrm{log}F+2.47\end{array}\right)}^{{\mathsf{T}}}\end{array}.\end{eqnarray} \tag{ 7 }$$

This is obtained by removing the third component P₃ from ${{X}}^{{\prime} }$. We use GMM to divide it into two subclasses. The result is shown in the upper left panel of Figure 3 where the two subclasses are represented in red and blue colors, respectively. In fact, the red and blue clusters are roughly composed of long and short GRBs, respectively, though the clusters are not single-handed chosen based on T₉₀. For convenience, we still use LGRB and SGRB to represent the red and blue data clusters, respectively.

Zoom In Zoom Out Reset image size

We then use the GMM classification result as the known training set, and use the LDA algorithm to find the optimal line that can distinguish the two subclasses, which is shown by an oblique line in the upper right panel of Figure 3. For a small number of GRBs near the line, LDA gives a different class prediction from GMM. We can also reduce the result to one-dimensional space by transforming the diagonal line to a vertical one and projecting all data into the horizontal axis, as shown in lower left panel of Figure 3. The conversion formula is:

$$\begin{eqnarray}&&L=\left({P}_{1}\ {P}_{2}\right)\cdot \left(\begin{array}{c}\,1.80\,\\ \,-0.70\,\end{array}\right)\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&=1.66\,\mathrm{log}\,{T}_{90}+0.84\mathrm{log}\,f-0.46\mathrm{log}F+3.24,\end{eqnarray} \tag{ 9 }$$

where the constant 1.66 just means a horizontal shift for the data points and has no meaning for the sample distribution characteristics. As for the classification, the important thing is the relative value of L, not the absolute value.

We use the new parameter L defined by Equation (9), which is combination of three observational parameters, as the improved criteria to distinguish short and long GRBs. As shown in the lower left panel of Figure 3, there is a significant bimodal structure. For comparison, we also fit T₉₀ in the lower right panel of Figure 3 where the separation of two Gaussian components is smaller. The degree of improvement can be quantitatively described by resolution, which is defined by:

$$\begin{eqnarray}&&R=\displaystyle \frac{2\times ({t}_{2}-{t}_{1})}{1.70\times ({W}_{1,h/2}+{W}_{2,h/2})}\end{eqnarray} \tag{ 10 }$$

where t₁, t₂ are the peak positions of the two components, and W_1,h/2, W_2,h/2 are the full widths at half maximum (FWHM) of the two components, respectively. We get these parameters by fitting these distributions by two Gaussian components. The resolutions of SGRBs and LGRBs are R_T ≈ 0.49 based on logT₉₀ and R_L ≈ 0.58 based on L, respectively. The classification resolution has been improved with the new criterion.

3.3. Another Possible Classification Scheme

Much research on GRB classification indicated different models with different numbers of subclasses. The model selection can be performed with multicomponent GMMs using the Bayesian information criterion (BIC), which can effectively prevent the high model complexity caused by the high accuracy of the model. In order to improve the confidence of the results, we used the bootstrap with replacement sampling method to extract 10,000 samples, and then performed the distribution of BIC scores. As shown in Figure 4, we performed from one- to eight-component models and found that a three-component model with the lowest BIC score is the optimal result.

Zoom In Zoom Out Reset image size

The three-component classification scheme is shown in the right panel of Figure 4. From the equal-probability lines, two protrusions extending from the highest density center cluster to two different directions that are caused by the other two clusters. LDA attempts to model the difference between the classes of data and it could project data to direction on which the different classes separated as much as possible while the same classes concentrated as much as possible. By this method the data could be transferred into new coordinate system, as shown in Figure 5, whose axes are obtained by LDA according Fisher criteria. The axes of new coordinate system are:

$$\begin{eqnarray}\begin{array}{ccc}({L}_{1}^{\mbox{'}}\ {L}_{2}^{\mbox{'}}) & = & \left({P}_{1}\ {P}_{2}\right)\cdot \left(\begin{array}{c}\,1.877\,-0.005\,\\ \,-0.010\,\,2.399\end{array}\right)\\ & \approx & \left(1.9{P}_{1}\ \ 2.4{P}_{2}\right)\end{array}.\end{eqnarray} \tag{ 11 }$$

Notice that the nondiagonal elements of the coordinate transfer matrix in Equation (11) are almost negligible; this means that the coordinate transfer hardly rotates the original coordinate system. It is an interesting coincidence that the first principal component is almost identical to the first direction of LDA, even though they are obtained through different criteria (maximum projection variance for PCA and Fisher criteria for LDA). This may be caused by physical conditions of different origin of GRBs.

Zoom In Zoom Out Reset image size

Indeed, the highest density center gathers most of the long GRBs with lower luminosity (the red cluster) and one of the protrusions points to the short GRBs area (the green cluster); the other protrusions of blue cluster points to the long GRBs with higher luminosity. For convenience, we use short GRBs (SGRBs), long-bright GRBs (LBGRBs), and long-faint GRBs (LFGRBs) to refer to the three subclasses, although here short and long is no longer strictly identified by T₉₀.

The three subclasses correspond to three clusters as shown in lower left panel of Figure 5 whose axes are ${L}_{1}^{{\prime} }$ and ${L}_{2}^{{\prime} }$. The upper and right panels show the histogram of ${L}_{1}^{{\prime} }$ and ${L}_{2}^{{\prime} }$, respectively. The distribution of each subclass is fitted with a Gaussian function, the histogram of total sample is shown fitted by a black line, which is the sum of the three Gaussian function. One can easily see that the resolution from projection distribution on ${L}_{1}^{{\prime} }$ is better than the projection distribution on ${L}_{2}^{{\prime} }$.

For further detail of the statistical characters of the three subclasses, we do the comparison of distribution of $\mathrm{log}{T}_{90},\mathrm{log}f$ and $\mathrm{log}F$ of the three subclasses; the result is shown in Figure 6. Obviously, LBGRBs have greater statistical fluence and peak flux than LFGRBs, and T₉₀ of SGRBs is statistically shorter than others.

Zoom In Zoom Out Reset image size

A potentially relevant factor may be the number of GRB pulses. The normalized cumulative distribution of the number of pulses for three subclasses are shown in Figure 7. According to our classification results, more than 96.6% of SGRBs have only one pulse, which for LFGRBs and LBGRBs are 85.1% and 64.3%, respectively. The KS-test derived very small P values, which indicates that three subclasses of GRBs have distinct distributions of the pulses numbers. LBGRBs statistically have more pulses than LGFRBs which may be the direct factor leading to the higher f and F. However, it is not yet known whether the number of GRBs pulses involves different physical mechanisms. More in-depth theoretical research is still needed.

Zoom In Zoom Out Reset image size

The connection between these parameters points to the similarity of each energy release process of GRBs. Under such rough assumptions, we can describe the process of GRB prompt emission by:

$$\begin{eqnarray}&&{f}^{{\prime} }={\bar{e}}_{\mathrm{ph}}^{{\prime} }\cdot {\bar{F}}^{{\prime} }\cdot {T}_{\mathrm{eff}}^{{\prime} }\end{eqnarray} \tag{ 12 }$$

where ${T}_{\mathrm{eff}}^{{\prime} }$ represents the effective duration of each GRB, ${\bar{F}}^{{\prime} }$ and ${f}^{{\prime} }$ represent the average photon flux during the burst and the fluence (flux integrated over the burst duration), ${\bar{e}}_{\mathrm{ph}}^{{\prime} }$ represents average energy of each photon, which can be considered equal for different bursts. All the prime marks indicate the parameters valued in local coordinate systems.

Considering the expansion effect of the universe, in the observer coordinate system, Equation (12) can be written as:

$$\begin{eqnarray}&&\frac{{f}^{{\prime} }}{(1+z)}=\frac{{\bar{e}}_{\mathrm{ph}}^{\mbox{'}}}{(1+z)}\cdot \frac{{\bar{F}}^{{\prime} }}{(1+z)}\cdot {T}_{\mathrm{eff}}^{\mbox{'}}(1+z).\end{eqnarray} \tag{ 13 }$$

Substitute the local variables to observer variables and Equation (13) can be written as:

$$\begin{eqnarray}&&f=\frac{{\bar{e}}_{\mathrm{ph}}^{\mbox{'}}}{(1+z)}\cdot \bar{F}\cdot {T}_{\mathrm{eff}}.\end{eqnarray} \tag{ 14 }$$

In view of the complexity of parameter definitions and measurements, we cannot directly measure $\bar{F}$ and T_eff. However, we can substitute them by common parameters according to $\bar{F}\propto {F}_{p}^{\alpha }$ and ${T}_{\mathrm{eff}}\propto {T}_{90}^{\beta }$ and then we have:

$$\begin{eqnarray}&&f\propto {\left(1+z\right)}^{-1}\cdot {F}_{p}^{\alpha }\cdot {T}_{90}^{\beta }.\end{eqnarray} \tag{ 15 }$$

Since the measurement of redshift is a strong bias on brightness, and most measured redshift located in 1 ∼ 3, ignoring the redshift factor can greatly increase the statistical sample size. So the upper correlation could be roughly expressed as:

$$\begin{eqnarray}&&{\rm{l}}{\rm{o}}{\rm{g}}f\,\text{\unicode{x0223C}}\,\alpha {\rm{l}}{\rm{o}}{\rm{g}}{F}_{p}+\beta \,{\rm{l}}{\rm{o}}{\rm{g}}\,{T}_{90}+const,\end{eqnarray} \tag{ 16 }$$

which is formally consistent with Equation (1).

Based on the above discussion, we use the Monte Carlo (MC) method to simulate this correlation and classification by generating 3000 GRB light curves. The single-pulse profile of GRBs can be expressed by the function (Norris et al. 2005; Hakkila & Nemiroff 2009; Nemiroff 2012):

$$\begin{eqnarray}&&I(t)=A\lambda {e}^{[-{\tau }_{1}/(t-{t}_{{\rm{s}}})-(t-{t}_{{\rm{s}}})/{\tau }_{2}]}\end{eqnarray} \tag{ 17 }$$

where A is the amplitude, t is the time since trigger, t_s is the start time of the pulse, τ₁ and τ₂ are characteristics of the pulse rise and decay, and $\lambda ={\exp }[2{({\tau }_{1}/{\tau }_{2})}^{1/2}]$. Two parameters equivalent to τ₁ and τ₂ are μ = τ₁/τ₂ and ${\tau }_{p}=\sqrt{{\tau }_{1}{\tau }_{2}}$.

Equation (17) gives the GRB prompt emission containing a single pulse. Based on this, we can get the multi-pulse GRB light curves:

$$\begin{eqnarray}&&{I}_{n}(t)=\displaystyle \sum _{i=1}^{n}{A}_{i}{\lambda }_{i}{e}^{[-{\tau }_{1,i}/(t-{t}_{{\rm{s}},i})-(t-{t}_{{\rm{s}},i})/{\tau }_{2,i}]},\,\,n\geqslant 1.\end{eqnarray} \tag{ 18 }$$

The simulation procedure is summarized as follows: First, simulate a GRB prompt light curve by randomly inputting A_i , τ_1,i , τ_2,i , t_s,i and n into Equation (18); second, calculate parameters (T₉₀, f, and F) from simulated light curve; finally, repeat the above steps by adjusting the initial parameter distributions and make the result consistent with the real data. To this end, we need two sets of independently distributed parameters, one set is used to generate SGRBs, the other set is used to generate LGRBs which is divided into two categories according to the number of pulses. Finally, the following parameter distribution assumptions (not from real statistics) is adopted:

• 1.  
  τ₁ and τ₂: Since μ = τ₁/τ₂ and ${\tau }_{p}=\sqrt{{\tau }_{1}{\tau }_{2}}$ more intuitively reflect the characteristics of light curve profiles, we generate random μ and τ_p instead of τ₁ and τ₂ directly. For LGRBs, μ ∼ uniform(0.1, 0.5) and $\mathrm{log}{\tau }_{p}\sim N(0.7,0.5)$. For SGRBs, μ ∼ uniform(0.1, 0.5) and $\mathrm{log}{\tau }_{p}\sim N(-1,0.5)$.
• 2.  
  A: Intuitively, peak of a single pulse. For SGRBs, $\mathrm{log}A\sim N(0.6,0.5)$ and for SGRBs $\mathrm{log}A\sim N(1,0.5)$.
• 3.  
  t_s : For SGRBs, t_s ∼ N(2, 12) and, for SGRBs, t_s ∼ N(0.5, 1.5).
• 4.  
  n: Number of pulses for GRBs modeled by integer of Erlang distribution Erlang(1, 1).

In order to make the GBM data consistent with the simulated data, two other parameters should be adopted according to Equation (14), i.e., average energy of each photon(${\bar{e}}_{\mathrm{ph}}^{{\prime} }$) and redshift (z). The average energy of each photon for LGRBs and SGRBs are 1.2 × 10⁻⁸ erg photon⁻¹ and 3 × 10⁻⁸ erg photon⁻¹. Redshift is generated for each GRB randomly according to the probability density function Erlang(2, 1) (Shao et al. 2011).

In the left panel of Figure 8, as comparison, the black circles represent real GRB data from GBM, green dots represent simulated SGRBs and red (blue) dots represent simulated long GRBs with single (multiple) pulse(s), simulated data perfectly satisfies the same correlation as real data. From the perspective perpendicular to the left panel, the right panel of Figure 8 shows the most discrete distribution of all GRBs. In addition to SGRB as an independent cluster, the cluster of LGRBs with multiple pulses also deviates from the LGRBs with a single pulse. This result confirms the GMM classification scheme in Figure 4 and the pulse statistics of each subclass reflected in Figure 7.

Zoom In Zoom Out Reset image size

There are a small number of simulated samples that exceed the data distribution at the location of low peak fluxes and fluences, which are most likely to correspond to GRBs that were not observed by the bias effect.

In summary, the difference between LFGRBs and LBGRBs is ostensibly derived from the light curve model. Statistically, those GRBs with multiple pulses have higher peak fluxes and fluences. As for the essential difference between the two subclasses of LGRBs, more dimensional data is needed for research.

In this work, we carry out some analysis on various distributions of GRBs by the data mining method. We first extract a tight correlation by using PCA on three-dimensional data. The third principle component contributes the least and has ignorable variance ratio, which indicates the intrinsic consistency of GRBs.

Based on the first and second principle (P₁ and P₂) components, which contain most of the variance, we composed classification analysis of two- and three-subclass schemes.

In order to better distinguish SGRBs and LGRBs, we reduced the dimensionality of the data by LDA and found a compound parameter (L) that provided a higher resolution between SGRBs and LGRBs. The resolution of the two subclasses measured by L is 0.58, it is obviously improved from 0.49, which is given by T₉₀. For the three-component model, we found the first principle component P₁ is almost identical to the first direction ${L}_{1}^{{\prime} }$ of LDA, which is obtained through Fisher criteria with the purpose of maximizing the covariance between classes and the minimum covariance within classes. This may be purely a coincidence in the data distribution or may be caused by different physical conditions where GRBs originated.

In general, the only uncontroversial and conclusive subclass of GRBs is SGRBs. Although the distinguish at the parameter boundary is controversial, there is a lot of evidence for its existence, which is widely recognized. Regarding LGRBs and their subclasses, we still need more information introduced by other independent parameters to help solve this problem.

Regardless of the number of parameters, PCA and LDA can be used to reduce the dimensionality to one or two if the number of subclasses is not more than three. It is easier to visualize data distribution in low-dimensional space.

The three subclasses classification scheme has been discussed in different aspects (Virgili et al. 2009; Bhave et al. 2017; Horváth et al. 2019; Svinkin et al. 2019). Classifying LGRBs based on hardness ratio, peak flux, or luminosity have similar conclusions as our analysis. However, such a classification scheme is still not sufficiently supported by the combination of theory and data.

So far, the existence and origin of SGRBs is widely recognized in that the classification based on T₉₀ points to the difference in the physical scales of the GRB progenitors. The influence of other physical parameters, such as Lorentz factor (Γ), the fraction of energy of shocked electrons (_e), or magnetic field (_B) on radiation is more complicated. Complex physical information is coupled in many observed parameters; this indicates that it may be possible to derive the physical information of the source through the combination of observed parameters. It may be helpful to introduce more independent parameters containing more physical information, and the analysis method of this work is easy to extend by introducing new independent variables. However, many parameters of SGRBs or LFGRBs are difficult to measure, such as the redshift and related quantities whose measurements require higher brightness.

As a supplement to pure data analysis, we generate a sample by the MC method based on an analytical light curve model (Norris et al. 2005; Hakkila & Nemiroff 2009; Nemiroff 2012), and then reproduce the strong correlation and clusters distribution. The simulation confirms that, under the unified light curve model, there must be a strong correlation as in Equation (1). And finally we found that the light curves composed by multiple pulses lead to statistically higher fluence and peak flux. This also explains why the LGRBs are composed of two subclusters, and consistent with the pulse distribution shown in Figure 7.

Another point worth noting is that, considering the definition of statistical parameters, if the research object is replaced with single GRB pulses, LBGRBs and LFGRBs may be unified into one category. And there are still other possibilities, such as the existence of fake LFGRBs that appeared purely due to background fluctuations. The fake triggers caused by background fluctuations are concentrated in low-brightness areas, which increased the proportion of single-pulse GRBs. If this part is removed, LFGRBs may have similar distributions of number of pulses as LBGRBs, i.e., in that case, the necessity to divide LGRBs into two subclasses still needs more in-depth discussion.

This work was funded by the National Natural Science Foundation of China (No. 11803006, 11903012); Science and Technology Project of Hebei Education Department (No. QN2018036) and Hebei NSF (No. A2019205166).