Classification of Photospheric Emission in Short GRBs

In order to better understand the physical origin of short duration gamma-ray bursts (GRBs), we perform time-resolved spectral analysis on a sample of 70 pulses in 68 short GRBs with burst duration $T_{90}\lesssim2$~s detected by the \textit{Fermi}/GBM during the first 11 years of its mission. We apply Bayesian analysis to all spectra that have statistical significance $S\ge15$ within each pulse and apply a cut-off power law (CPL) model. For our analysis, we identify the maximal values of the low energy spectral index, $\alpha_{\rm max}$. Under the assumption that the main emission mechanism is the same throughout each pulse, such an analysis is indicative of pulse emission. We find that $\sim$1/3 of short GRBs are consistent with a pure, non-dissipative photospheric model, at least, around the peak of the pulse. This fraction is larger compare to the corresponding one (1/4) obtained for long GRBs. For those bursts which are compatible with a pure photospheric origin, we find (i) a bi-modal distribution in the values of the Lorentz factors and the hardness ratios; (ii) an anti-correlation between $T_{90}$ and the peak energy, $E_{\rm pk}$: $T_{90} \propto E_{\rm pk}^{-0.35\pm0.01}$. This correlation disappears when we consider the entire short GRB sample. Our results thus imply that short GRB population may in fact be composed of two separate populations: a continuation of the long GRB population to shorter duration, and a separate population with distinct physical properties. Furthermore, thermal emission is initially ubiquitous, but is accompanied at longer times by additional radiation (likely synchrotron), making time-resolved spectroscopy absolutely necessary in analyzing GRB pulses.


INTRODUCTION
After four decades of extensive research, the physical origin of gamma-ray burst (GRB) prompt spectra remains unclear and highly debated. The classification of GRBs is a tool which could help to understand the emission mechanisms at work. The main classification into short-hard and long-soft is based on their duration and spectral hardness. The short GRBs have a duration shorter than 2 seconds while the long GRBs have a duration longer than 2 seconds (Kouveliotou et al. 1993). The spectral peak energy of short bursts is, on the average, higher than that of long GRBs (e.g., Ghirlanda et al. 2011). However, both classes share many spectral characteristics. For instance, their spectra peak in the MeV range, with power law extensions below and above the peak. Both populations have a common inverse correlation between the duration and intensity for individual pulses (Hakkila & Preece 2011;Norris et al. 2011), and they follow a similar relation between the peak energy, E pk and the peak luminosity, L peak as well as the isotropic equivalent energy, E iso (Yonetoku et al. 2004;Amati 2006;Ghirlanda et al. 2009).
Despite these observed similarities, short and long bursts are thought to originate from different progenitors; the collapse of a very massive star for long GRBs (Woosley 1993) and a compact binary merger for short GRBs (Eichler et al. 1989). In fact, long GRBs are studied more than short ones. Indeed, they release more photons which allows more detailed spectral studies. In addition, more redshifts are known for long GRBs than for short ones since the afterglow after a few thousands seconds is brighter for long bursts. This allows the study of intrinsic properties (e.g., Howell & Coward 2013). The recent increase of interest in the study of short GRBs is mostly due to the detection of short GRB 170817B simultaneously with the first gravitational wave (GW) from a merger of binary neutron stars (Abbott et al. 2017;Goldstein et al. 2017).
Observationally, many GRB prompt spectra have too narrow νF ν peaks compared to what is expected from the synchrotron emission model (Axelsson & Borgonovo 2015;Yu et al. 2015). Yet, they are broader than a Planck spectrum (Goodman 1986;Paczynski 1986;Beloborodov 2011). Photospheric emission from highly relativistic outflows is often used to explain this observed spectral shape. Broadening of the spectrum by energy dissipation below the photosphere can be caused by shocks, dissipation of magnetic energy or collisional processes (Giannios & Spruit 2005;Pe'er et al. 2006;Beloborodov 2010). Moreover, broadening in a passively cooled jet without any energy dissipation can be due to geometrical broadening occurring during the coasting phase (Beloborodov 2011;Bégué et al. 2013;Lundman et al. 2013). In order for the emission to be detectable the outflow has to become transparent below or close to the saturation radius, r s , where the outflow saturates to its final outflow Lorentz factor (Mészáros 2006;Ryde et al. 2017).
The observed spectral shape of the prompt emission is commonly characterised by empirical models, such as the "Band" model (Band et al. 1993) or a cutoff power-law model (see, e.g., Yu et al. 2016). However, in making the link between observation and theory, the parameters of the empirical models should not be used directly for the comparison with the prediction of physical emission models. Indeed, attempt to make such a link leads to two main problems. The first one is known as an energy-window bias effect. When the empirical model does not match the true spectral shape (its curvature) then physical interpretation of the model parameters will be wrong; e.g, a positive correlation between parameters of empirical model at low energies (Preece et al. 1998;Lloyd & Petrosian 2000;Burgess et al. 2015;Ryde et al. 2019;Acuner et al. 2019). The second problem is the limitation due to the band-width of the detector which prevents the full spectrum to be detected (Burgess et al. 2015;Ryde et al. 2019).
There are two solutions to overcome these problems. The first one is to use a physically motivated model and fit it directly to the data (e.g., Lloyd & Petrosian 2000;Ahlgren et al. 2015). In this way, there is no need for an empirical function. However, it is computationally expensive due to the need to make a forward-folding of the theoretically generated spectra through the detector's response matrix, and the need to subtract the background -both vary from burst to burst. Thus, the claimed model has to be fitted individually to each burst. Furthermore, one has to assume knowledge of the physical model that should be used (e.g., Baring & Braby 2004;Burgess et al. 2016Burgess et al. , 2019bOganesyan et al. 2019). Due to these limitations, this direct method was applied, so far, only to a limited number of bursts (e.g., Vianello 2018; Ahlgren et al. 2019).
The second solution is to use an assumed physical model to generate a large number of synthetic spectra which are, in turn, fitted with empirical functions. This provides the distribution of the empirical model parameters that the given theoretical model corresponds to. The properties of the parameter distributions, for instance their widths, depend on how well the empirical model matches the theoretical model. These distributions can then be compared to the full GRB catalogue, in order to assess the theoretical model's ability to explain the data. This method was used by several authors (e.g., Burgess et al. 2015;Acuner et al. 2019) to make statistical claims about the ability of a theoretical model to fit the data.
In an attempt to fit a non-dissipative photospheric model (Beloborodov 2011;Lundman et al. 2013) to GRB spectra, Acuner et al. (2019) followed the second method and generated a series of synthetic spectra with a high signal-to-noise ratio (SNR) of 300 and peak energies at the range of 40-2000 keV. The simulated (synthetic) spectral data were fitted with a cutoff power-law model. It was found that the distribution of the low-energy photon indexes ranges from -0.4 to 0.0 and peaks at around -0.1. This was then compared with the distribution of the maximal, time-resolved value of the low energy spectral index, α in each long burst in the sample of Yu et al. (2016Yu et al. ( , 2019. They found that 1/4 of the long bursts have an α max which is consistent with a non-dissipative outflow, releasing its thermal energy at the photosphere. However, Acuner et al. (2019) did not consider short bursts since the selection criteria of Yu et al. (2016 is mainly based on bright burst with duration T 90 2 s. While the spectral properties of short GRBs are much less studied than that of long GRBs, evidence are accumulating that photospheric (thermal) emission could play an important role in these bursts as well. The main motivation for our current study is the large number of short GRBs seen in the cluster 5 in Acuner & Ryde (2018). This cluster was found to be consistent with a photospheric emission origin. Therefore, in this work, we are also using the fitted synthetic spectra from Acuner et al. (2019) to find the fraction of short GRBs compatible with a non-dissipative photosphere (NDP) model. As a first step, we apply time-resolved analysis to the spectra of individual pulses obtained from 68 short GRBs by using the full Bayesian analysis method. As a second step we study, in detail, the time bins with the hardest low-energy spectral index (maximal value of α) in each pulse. This paper is organised as follows. In Section 2, we define the sample of short GRBs and present the analysis methods. In Section 3, we present result of the spectral parameter relations, the observed α distributions, the Lorentz factor for the bursts consistent with thermal emission, and the hardness ratio. In Section 4, we then discuss our choice of spectral fitting model, and the correlations between temporal and spectral structures. Finally, in Section 5 we list our summary and conclusions.

Sample Selection
We select short GRBs, namely GRBs having duration T 90 2 seconds detected by the Gamma-ray Burst Monitor (GBM) onboard the Fermi Gamma-ray Space Telescope during the first 11 years of its mission. We scan all the short bursts for which automatic spectral fits are performed on the time-resolved data around the peak flux, within the time interval given in the GBM catalogue (Von Kienlin et al. 2014). We find a total of 147 short bursts for which spectral fits can be carried out and analysed. All the data is taken from the Fermi /GBM burst catalogue published at HEASARC 1 . We set a limit for at least one bin to have a signal-to-noise ratio (or statistical significance, see Section 2.2 for definition), S ≥ 15 in each pulse/burst; see appendix A for selection method of the signal-to-noise ratio. We end up having 70 pulses from 68 short GRBs as a final sample listed in Table 1.

Analysis Method
For the analysis, we follow the procedure of the Fermi /GBM GRB time-integrated (Goldstein et al. 2012;Gruber et al. 2014) and time-resolved catalogues (Yu et al. 2016. We select at most three NaIs and one BGO for the spectral analysis of each short GRB, see in Table 1 column 3. We use the response files RSP (except for two cases, GRB 090510, GRB 170127, in which the RSP2 files are used) for each short GRBs.
We select the source interval from the first few seconds of the burst light curve where the first pulse is most prominent, see in Table 1 column 4. Indeed, most of the bursts in the sample are single pulsed burst. We use a NaI detector in which the largest photon counts per second was recorded from the burst to define the background intervals before and after the pulse, see in Table 1 columns 3, 5 and 6 respectively. These intervals are then applied to all detectors. As a standard procedure in GRB background fitting of GBM data, we fit a polynomial with the order of between 0 and 4, to each energy channel (128 channels for TTE) of each detectors. The optimal order of the polynomial is determined by a likelihood ratio test independently for each energy channel. The polynomial is interpolated into the source and integrated over the source interval to obtain the background photon count flux and its corresponding errors in each channel.
To perform the time-resolved spectroscopy, we used the Bayesian spectral analysis package the Multi-Mission Maximum Likelihood 3ML (Vianello et al. 2015). We first applied the Bayesian Blocks method (Scargle et al. 2013) in order to rebin the light-curves into adequate intervals in which any evolution is small and to identify the statistical significance of each spectra. The statistical significance (S) adopted in the current work is a test statistic that incorporates the information of signal-tonoise ratio and suitable for Poisson sources with Gaussian backgrounds, see Vianello (2018); Yu et al. (2019) for further details. We find that the spectral parameters are typically well constrained for bins with statistical significance, S ≥ 15. Therefore, for each burst, we rebinned the TTE light curve of the brightest NaI detector into Bayesian Blocks with a correct detection rate for a single change point of p 0 = 0.01 (Scargle et al. 2013). The Bayesian Block binning is then transferred and applied to all other detectors. The total number of bin for each pulses is listed in Table 1 column 7.
We use the standard Fermi /GBM energy ranges: 8−30 keV and 40 keV to ∼850 keV for the NaI detectors (avoiding the K-edge at 33.17 keV) 2 , and ∼250 keV to 40 MeV for the BGO detectors. In the Bayesian spectral analysis, we fit the data using an empirical model with Poisson likelihood for the source and Gaussian likelihood for the background. This model is a power law with an exponential cut-off (hereafter CPL) that has been extensively used and discussed as a best model to fit GRB data in the Fermi GBM catalogues (Goldstein et al. 2012;Gruber et al. 2014;Yu et al. 2016Yu et al. , 2019. The CPL fit parameters are the normalization K(ph s −1 cm −2 keV −1 ), the low-energy power-law index α, the cut-off energy E c (keV). We further derive the CPL peak energy E pk (keV) and the CPL energy flux F (erg s −1 cm −2 ).
For the posterior simulations we use the informative priors for parameter intervals that can be detected with GBM and they are reasonable for capturing the shape of our seed function. While the normalization (K ∼ 10 −11 − 10 3 ph s −1 cm −2 keV −1 ) is assumed to have a log uniform prior, the low energy index (α ∼ −3 − +2) and the cut-off energy (E c ∼ 10 − 10000 keV) are assumed to have uniform priors. We visually checked the model fit variability by using the corner plots which display the posterior distributions of the resulting parameters. We point out that most of the data with a very low spectral index −2.1 < α < −1.1 belong to a single burst, the brightest short GRB 120323 (see Figure 1). As a result the values of α for the vast majority of bursts in our sample lie between -1.6 and +0.6 within 1-σ uncertainty. This range is narrower than that observed in long GRBs, −2 < α < 1 . However, the range of E pk (40 keV < E pk < 6 MeV) and flux (8 × 10 −6 erg s −1 cm −2 < F < 9 × 10 −3 erg s −1 cm −2 ) obtained from short GRBs are similar to the range obtained in long GRBs, 10 keV < E pk < 7 × 10 3 keV and 10 −7 erg s −1 cm −2 < F < 5 × 10 −5 erg s −1 cm −2 .
When comparing the results presented in Figure 1 to the results obtained by Yu et al. (2019), we conclude that most of the bins obtained from short GRBs generally have harder values of the spectral index α, higher energies and have higher fluxes than those in long GRBs.

Distribution of low energy spectral indexes
Time-resolved spectral catalogues present the evolution of the parameters over time bins (e.g., Kaneko et al. 2006;Yu et al. 2016Yu et al. , 2019. Their parameter relations for individual GRBs are interpreted by physical models (e.g., a qualitative photospheric emission scenario, Ryde et al. 2019). However, global parameter distributions, such as the α-distribution, contain varying number of spectra from each individual burst. This leads to a bias towards strong bursts with many time bins (e.g. GRB 120323; see the purple line in Figure 1).
We show the normalized histogram of the low energy spectral indexes, α, in Figure 2, left panel. The green curve is the kernel density estimation (KDE) of the distribution taken into account the errors (see appendix B for the KDE definition). The histogram contains 153 spectra. We compare the values of the low energy spectral index, α (considering 1σ lower limit) with the "synchrotron line of death" value (Preece et al. 1998;Kaneko et al. 2006), α = −2/3, which is shown by the red dashed line in Figure 2. We find that 56% of the analyzed spectra (within a 1σ error) violate the criteria set by the "synchrotron line of death".
Since the spectral index can vary between time bins of the same burst, it is difficult to draw firm conclusions on the emission mechanism based on the entire sample. The best way to constrain the emission mechanism during a pulse/burst is to select the time bin that contains the largest value of the spectral index α in each pulse/burst. The reason for this is that physical models typically have an upper limit on how hard the spectra can get. Therefore, it is enough that one single bin violates such a limit . Example of the temporal evolution (from light to dark blue) of the spectral index α, E pk and the calculated energy flux of each bin obtained from GRB 140209. The light curves (photon flux) are overlaid in gray color with Bayesian Block time bins. Data points with circle, square, diamond, and no shape indicate statistical significance S > 20, 20 > S > 15, 15 > S > 10, S < 10 respectively. In left panel, it is seen that the spectral index α peaks close to, or coincides with the peak of the light curve. This is typical to most pulses, see also Yu et al. (2019).
that the corresponding emission model is rejected by the data (under the assumption that a single emission mechanism is the source of the observed signal in the entire duration of the burst). This is discussed in several recent paper for long GRBs (e.g., Yu et al. 2019;Ryde et al. 2019;Acuner et al. 2019).
For each of the 70 pulses in the 68 short bursts in our sample we therefore select the time bin which contains the maximal (hardest) value of the low energy spectral index, which is denoted by α max . We find that for most pulses, α = α max close to, or at the peak of the light curve. As such, it contains the most valuable information of the spectra. This is demonstrated in Figure 3: the temporal evolution of α is shown in the left panel, with the corresponding temporal evolution of the peak energy (middle panel) and the temporal evolution of flux (right panel). The spectral index α max and its corresponding peak energy (E pk ), flux and statistical significance (S) are listed in the last four columns of Table 1.
We show in Figure 2 (right panel) the normalized histogram of α max in each of the 70 pulses from 68 short GRBs in the sample together with the KDE of the distribution. The red dashed line shows the "synchrotron line of death", α = −2/3. We find that 70% (within a 1σ error) of the pulses violate the "line of death" criterion, and are therefore better interpreted with a model that is not synchrotron, such as the photospheric model. This fraction is larger than the one obtained in the case of long GRBs, 60% within a 1σ error . We also note that the softest value is α max = −1.26 for GRB 160822.

On the consistency with non-dissipative photospheric model
We show in Figure 4 the relation of α max and E pk in all 70 pulses from 68 short GRBs in the sample. The light blue line corresponds to the values of α that are found when a CPL function is fitted to a non-dissipative photosphere (NDP) spectrum, peaking at different energies, E pk . The retrieved values of α are significantly smaller than the asymptotic slope of the theoretical NDP spectrum (α ∼ 0.4, Beloborodov 2010) due to (i) the limited energy band of the detector as well as (ii) the limitation of the CPL function to correctly model the shape of the true spectrum; see Acuner et al. (2019) for further details.
We find that 36% of the observed points (25/70) are consistent with being above the NDP line within 1σ error, and are therefore consistent with having a dominant quasi black body component. This fraction is significantly larger than the fraction found in the study of pulses from long GRBs, 26% (for single or multi pulses bursts) and 28% (for only single pulsed bursts) from both catalogues of Yu et al. (2016) and Yu et al. (2019) respectively. See Acuner et al. (2019) for further details.
Even though most of the short GRBs in the sample are single pulse bursts, two bursts in our sample (GRB 150819 and GRB 180703) are found to have two separate pulses (these are marked by blue and green points respectively in Figure 4). While the first pulse of GRB 150819 (blue color in Figure 4) shows a hard spectral index, α max = 0.25 +0.14 −0.17 , the spectral slope of the second pulse is much softer, α max = −1.02 +0.04 −0.05 . This might be an indication for a change in the leading emission mechanism, e.g., from photospheric emission to synchrotron (Zhang et al. 2018). On the other hand, both pulses of GRB 180703 (green color in Figure 4) are very hard, and are both compatible with NDP line. This might be an example of a burst in which a single emission mechanism is responsible for the full duration of a burst. In this case, the dominant emission mechanism throughout the full duration of the burst is likely a photospheric emission (Acuner & Ryde 2018).

Lorentz factor
If indeed the observed spectra above the NDP line have a photospheric origin, then one can use the data to calculate the coasting Lorentz factor, η. Here we estimate Table 1. A sample of 70 pulses from 68 short GRBs used in this study. Column 1: GRB names (bn), Column 2: T 90 , Column 3: detectors, the detector in brackets is the brightest one used for background and Bayesian block fitting. Column 4: source interval, Columns 5, 6: background intervals, Column 7: the number of time bins (N) using Bayesian blocks across the source interval. Column 8: hardest low energy spectral index α max . Column 9: significance (S) of the time bin with α max obtained for S ≥ 15. Column 10: corresponding peak energy (E pk ), Column 11: corresponding flux.  the Lorentz factor, η, for 25 pulses from 24 short GRBs above the NDP line by using equations 1 and 4 in Pe'er et al. (2007). As the redshift of most bursts in our sample are unknown, we have assumed redshift z = 1. We further assumed the flux in the analyzed time bin is equal to the black body flux, ∼ F BB and the observed temperature is related to the peak energy via E pk ∼ 1.48T ob . The flux and peak energy (E pk ) for each short GRB obtained in our analysis for the corresponding α max presented in Table 1. For comparison, we also compute the η for all 12 long GRBs found above the NDP line in the sample of Yu et al. (2019); Acuner et al. (2019). The normalized distributions of the Lorentz factor η for 25 pulses from 24 short GRBs and 12 long GRBs above the NDP line are presented in Figure 5. We find that the Lorentz factor (η mean = 775) of those short GRBs compatible with a photospheric origin are in the order of few hundred which is similar, though somewhat higher than that of the long GRBs (η mean = 416) (Pe'er et al. 2007;Racusin et al. 2011;Ghirlanda et al. 2018;Chen et al. 2018). Surprisingly, we find a bi-modal distribution in the values of η for short GRBs ( Figure 5, left panel). However, such a bi-modal distribution is not found in the analysis of long GRBs. The low peak coincides with the values obtained for long GRBs, while the higher peak is a factor of ∼ 3 higher. When cutting the sample at η c = 700, we find that the pulses with the high Lorentz factor, defined as "peak 2", have a corresponding higher E pk ( Figure 5, right panel). While this by itself may not be surprising, we point out that no such clear correlation is observed when analyzing the entire population of 70 pulses (see Figure 4).

Correlations in bursts above the NDP line
For these two groups from the bi-modal η distribution, a strong positive correlation between η and E pk is found (Figure 6, top left panel). This can be explained as due to the computational dependence, as η ∼ E 1/2 pk F 1/8 . However, interestingly, we also find an anti-correlation between T 90 and η (Figure 6, top right panel), as well as an anti-correlation between T 90 and E pk (Figure 6, bottom left panel), both unexpected. The Bayesian inference of the T 90 and E pk correlation shows that the slope is -0.35, the correlation coefficient is -0.59 and the fractional bias for the correlation coefficient is 0.01, within the 1σ errors of fitted parameters. We do not observe any correlation between T 90 and the spectral slope α max for bursts above the NDP line ( Figure 6, bottom right panel). Similarly, no such correlation is found when considering the entire sample of 70 pulses. However, when we consider all the sources having E pk > 800 keV (presented in Figure 4) we do observe an anticorrelation between T 90 and α max which is presented in Figure 7. This anti-correlation implies that sources with harder spectra have shorter T 90 . This suggests that there might always exist a thermal emission at short times, which is accompanied by other emission processes such as synchrotron at later times. If this interpretation is correct, the lack of thermal emission in a given GRB might be explained by the lack of observed or studied spectra at sufficiently short time.

Hardness ratio
The anti-correlation we found between E pk and T 90 of short GRBs with high E pk , motivated us to study a possible correlation between the hardness ratio (HR) and T 90 . Following Kouveliotou et al. (1993), we calculate the HR using the two typical energy bands, 100 − 300 keV and 50 − 100 keV. To integrate the spectra, we use the CPL fit parameters, α max and E c , for each given time bin in each of the 70 pulses from 68 short GRBs. For comparison, we added 38 pulses from 37 long GRBs (with T 90 2 s). These pulses were selected from single pulsed, long bursts that have at least 5 time bins in which the statistical significance is S ≥ 20. The time-resolved spectroscopy T 90 (s) peak 1 peak 2 Figure 6. Left top: η and E pk relations (the correlation coefficient is r = 0.92, the chance probability is p = 8.69 × 10 −11 ). Right top: η and T 90 relation (the correlation coefficient is r = −0.53, the chance probability is p = 0.006). Data points are the short GRBs (peak 1: black color and peak 2: green color) from the bi-modal η distribution. Left bottom: T 90 and E pk relations (The spearman correlation coefficient is r = −0.43, the chance probability is p = 0.033). The light blue line shows the mean of the posterior distribution from the Bayesian inference of the T 90 and E pk correlation. Within the 1σ errors of fitted parameters, the slope is −0.35, correlation coefficient is −0.59, fractional bias in correlation coefficient is 0.01. Right bottom: T 90 and α max relations. Data points are the short GRBs (peak 1: black color and peak 2: green color) from the bi-modal η distribution.
is performed in each time bin individually with the full Bayesian approach. Their spectra are fitted with the CPL model, similar to the analysis by Yu et al. (2019). The HR -T 90 relation is shown in Figure 8 (left panel), for both short and long GRBs. Short bursts with α max above the NDP model prediction are displayed in black while those below the prediction are in red. For the long bursts, those colors are blue and purple respectively. While the T 90 selection criteria enables to clearly discriminate the long and short GRB population, we do not observe any additional correlation in this plot.
The HR -T 90 relation for the 25 pulses from 24 short GRBs above the NDP line is shown in the right panel in Figure 8. The color coding (peak 1: black color and peak 2: green color) is the same as that used for the bi-modal η distribution in Figure 5. Now, a clear separation is observed: GRBs with lower peak energy, have low Lorentz factor, lower hardness ratio, and longer T 90 . . T 90 and α max relations (The spearman correlation coefficient is r = −0.37, the chance probability is p = 0.033) with a cut at peak energy, E pk =800 keV in Figure 4 .
Indeed, all the parameters of these GRBs in the first peak of the bi-modal η distribution (in Figure 5, black color) seem to form a continuous distribution of the parameters of the population of long GRBs. This is in contrast to GRBs in the second peak of the bi-modal η distribution (in Figure 5, green color), that have a higher HR than both the GRBs below the NDP line (45 pulses) as well as long GRBs (38 pulses). This result implies that the duration T 90 as a single criterion does not make a good separation between the two populations; rather as we show, the short GRBs may be composed of two separate populations -one which forms a continuation of the long GRB population, and another, separate population.

Spectral parameter correlations for the two groups
We find a weak positive correlations between the F and E pk (in Figure 9, left panel) and between the F and α max (in Figure 9, right panel) for bursts (in the two groups seen in the bi-modal η distribution) above the NDP line. However, no clear correlation is seen between these parameters for the bursts below the NDP line. This by itself is an interesting result. In literature, several publications claim that there is a strong correlation between the luminosity, L peak or the isotropic energy, E iso and peak energy, E pk (e.g., Yonetoku et al. 2004;Amati 2006;Ghirlanda et al. 2009). These claims are based on a large sample of long GRBs. In contrast, here we do not find any such  Flux (erg s 1 cm 2 ) peak 1 peak 2 Figure 9. Left: The flux and E pk relations (The spearman correlation coefficient is r = 0.40, the chance probability is p = 0.05). Right: The flux and α max relations (The spearman correlation coefficient is r = 0.39, the chance probability is p = 0.06). Data points are the short GRBs (peak 1: black color and peak 2: green color) from the bi-modal η distribution.
correlation when considering the entire sample of short GRB pulses, but we do find a correlation when we consider only those short GRBs whose spectral slope is above the NDP line. In this work, we fitted the Fermi /GBM data using the phenomenological cut-off power law (CPL) model (this is also known as the "Comptonized" model). Several empirical models are commonly used in the literature for spectral analysis of GRBs. In addition to the CPL model (Kaneko et al. 2006), these include the "Band" model (Band et al. 1993) as well as the smoothly broken power-law model (Ryde et al. 1999). Yu et al. (2016) showed that the CPL is the preferred model for the majority (70%) of bursts, according to the Castor C-Statistic (CSTAT; a modified version of the original Cash statistic derived by Cash (1979)). In addition, a consistent result was found by Yu et al. (2019) based on the Deviance Information Criterion (DIC) in Bayesian statistics (Spiegelhalter et al. 2002). This means that the results of the CPL model for these bursts have lower DIC and higher effective number of parameters, p DIC > 0 (Gelman et al. 2014) than those of the "Band" model. Additionally, the resulting parameters for the CPL fits are constrained within the prior ranges more often than those obtained from the BAND function fits; see also Burgess et al. (2019a,b).
It was recently argued by Burgess et al. (2019b) that a direct fit of the data to a synchrotron model enables to overcome the "line of death" criteria in many GRBs. However, this idea suffers several severe drawbacks. First, the bursts selected in this work are limited to long and smooth bursts, with known redshift, which are different than the short pulses considered here. Second, the values of the parameters found in their fits require unacceptable high ratio of explosion energy to ambient mass density, of more than 7 orders of magnitude (E/10 53 erg)/(n ism /1cm −3 ) 4 × 10 7 than the highest observed so far. In order to overcome this problem, Burgess et al. (2019b) suggested an additional acceleration of particles within the relativistically expanding jet ("jet within a jet"); however, no such mechanism that can lead to relativistic expansion within an already relativistically expanding jet is known. We thus find that this model is still incomplete, and an interpretation of an empirical fit still provides better insight. Another suggestion was given by Ghisellini et al. (2019) based on the low energy break in the prompt spectrum of GRBs. They argued that the emission process is still the synchrotron radiation but produced by protons and cannot be completely cool.

Correlation between temporal and spectral structures
When we consider the entire sample of 70 pulses analyzed in the short GRB population, we do not observe a correlation between the burst duration T 90 and the peak energy, E pk . Similarly, for the long GRBs that we considered (both below and above the NDP line) no clear T 90 −E pk correlation was detected. However, when we consider only those short bursts that have a hard value of the spectral index, α such that they are above the NDP line (Figures 4, 5) we do find an inverse correlation between T 90 and peak energy, T 90 ∝ E −β pk , where β = 0.35. A quantitative relationship between the temporal and spectral structure in gamma ray bursts has been considered by several authors in the past. However, these works treated only the long GRBs. Richardson et al. (1996); Bissaldi (2011) found a T 90 ∝ E −β pk correlation with β 0.4. Their samples contained only bright, long GRB population. When considering the entire set of long GRB population, Qin et al. (2013) report a similar correlation, but with weaker dependence, β = 0.2.
Here we report, for the first time, such a correlation in the sample of short GRBs. This could not have been done in the past, due to the small sample, or lack of suitable method to study the spectra. The similarity between the correlation found here for short GRBs above the NDP line and for the bright long GRBs [which have hard spectral index, α; see Ryde et al. (2019)], as well as the fact that we do not detect any correlation for bursts below the NDP line, suggest a possible correlation between the emission mechanism and the burst duration. Bursts above the NDP line are consistent with originating from the photosphere, hence the photons directly probe the inner engine. While bursts whose spectra are below the NDP line may have additional radiative mechanisms, such as synchrotron emission, which originates from the outer regions of the outflow (outside the photosphere) and as such do not necessarily follow directly the duration of the inner engine. If this interpretation is correct, it points to a possible correlation between the duration of the inner engine and the temperature -or total energy, of the released photons. This further points to the importance of spectral analysis in analyzing possible correlations in the GRB population. A second correlation we find is between T 90 and α max when we consider a cut at higher peak energy (E pk > 800 keV) (see Figure 7). This (anti-) correlation further suggests a dual emission mechanism: short duration GRBs might be dominated by a thermal component, while an additional emission process may both lead to shallower spectra and be characterized by a longer duration.
The results we find therefore strongly support the idea that the spectra of both short and long GRBs contain (at least) two separate components: a photospheric emission component that correlates directly with the inner engine activity, and a second component, possibly having a synchrotron origin, that is longer in nature, and less steep.

SUMMARY AND CONCLUSION
In this work, we have selected a sample of 70 pulses from 68 short GRBs with T 90 2s detected by Fermi/GBM. These GRBs have at least one time bin with statistical significance S ≥ 15. Time bins were selected using the Bayesian Block method. A total of 153 time-resolved spectra were obtained and fitted with empirical CPL photon model using a full Bayesian method.
We investigate the distribution of the maximal (hardest) value of the spectral index α in each of the pulses, denoted α max . Assuming that a single emission mechanism dominates throughout each pulse, the maximal value of the spectral index, α max provides a useful information on this emission mechanism. We find that 70% (within a 1σ error) of short GRBs have at least one interval in which the value of α is beyond the value allowed by the "synchrotron line of death" (see Figure 2). These values of α max are typically obtained when the flux is close to its peak (see Figure 3). Therefore the emission mechanisms in these pulses are inconsistent with being dominated by synchrotron emission.
When considering the intervals for which α = α max , we find that 36% (within a 1σ error) of the spectra are consistent with having a non-dissipative photospheric origin, namely are above the NDP line (Acuner et al. 2019). This is presented in Figure 4. These numbers are slightly higher than that of long bursts. Indeed, short bursts have been found earlier to be harder than long bursts (Kouveliotou et al. 1993;Tavani 1998).
For these bursts compatible with a non-dissipative photospheric origin, we calculate the coasting Lorentz factor, η, and find a bi-modal distribution in the values of η (see Figure 5), peaking around ∼ 300 and ∼ 1000. The first peak (η pk,1 ∼ 300) is compatible with the average Lorentz factor η found in long GRB population (Racusin et al. 2011) while the second peak (η pk,2 ∼ 1000) is larger by a factor of 3.
When considering those bursts above the NDP line, we further find a strong positive correlation between η − E pk and a negative correlation between T 90 − E pk with a slope −0.35 ± 0.01 (see Figure 6). A clear separation between bursts that belong to the two distinct distribution in η is further observed in their duration (T 90 ), peak energies (E pk ) and hardness ratio (see Figure 8).
These results thus suggest that short GRBs might in fact contain two separate populations. These results are consistent with earlier work of Fenimore et al. (1995), who studied pulses in only bright events and conclude that the pulses observed at higher energies are shorter and have narrower peaks.
We find an anti-correlation between T 90 and α max when we consider a cut at large peak energy (E pk > 800 keV), see Figure 7. We therefore conclude that here in our sample most pulses are compatible with thermal emission at short times but with some contamination from other emission processes such as synchrotron at later times.
Our results provide a direct indication that the GRB duration by itself is not sufficient to classify the nature of a GRB: T 90 2s or T 90 2s by itself is not enough to separate short and long GRBs. Rather, one needs to consider additional information, which includes a spectral information, such as the hardest value of α in each pulse/burst, and the corresponding Lorentz factor, η. Indeed, classification of GRBs is long discussed in the literature as a way of discriminating GRB progenitors (Kouveliotou et al. 1993;Tarnopolski 2015, and references there in). Here we show that the maximal (hardest) value of α in each pulse/burst can be used as an additional method for the classification of bursts, especially for the classification of short ones.
The bi-modal distribution we find in the values of the Lorentz factor, together with the differences in the harness ratio, provides a strong indication that what is currently classified as short GRBs, in fact is made of 2 separate populations. The first is an extension of the long GRB population to shorter duration, and the second is a truly separated population. A striking result is the difference in the Lorentz factors, by an average factor of 3, with this separate population having Lorentz factor of ∼ 1000, and in some cases higher. This implies that, on the average, the outflows of the separate population contain much less ejected material than that of long GRBs, which provides a further clue to the true nature of short GRB progenitors.
Our results prove the importance of using time resolved spectral analysis to access the physical information of GRBs. In order to do that, we use the new method of investigating the maximal (hardest) value of the spectral index α in each of the pulses. This is valid, under the assumption that the same emission mechanism dominates throughout an individual pulse. Doing so we find that 36% (within a 1σ error) of the spectra are consistent with having a non-dissipative photospheric origin. Our method lead us to find a correlation between the burst duration (T 90 ) and the peak energy (E pk ) of the spectra but only for those bursts consistent with a photospheric origin.
These surprising results lead us to conclude the following: (1) a thermal (photospheric) emission is ubiquitous among short GRBs, with ∼ 1/3 being consistent with having a pure thermal origin, and another large fraction may also have a thermal origin, which is distorted by sub-photospheric energy dissipation. However, this component is often accompanied by an additional emission mechanism (likely, synchrotron) which makes it hard to separate and clearly identify the dominant mechanism. (2) At early (short) times, the thermal component often dominates, but at longer times it is accompanied by a second mechanism, which makes it sub-dominant. (3) Only for those bursts in which the thermal component dominates, we find a correlation between pulse (burst) duration T 90 , and the peak energy, E pk which corresponds to the temperature: higher peak energy correspond to shorter burst duration. Since no corresponding correlation is found in the flux, this implies that a similar amount of energy is released in short time, which leads to higher temperature. This result may therefore provide a very strong hint towards a better understanding of the progenitor models and explosion mechanisms in short GRBs. (4) When considering only those bursts with high peak energy, E pk > 800 keV, we further find a correlation between the burst duration T 90 and the hardest spectral slope α max . This further supports the idea of a dual emission mechanisms: thermal and non-thermal (synchrotron).