Self-Similarities and Power-laws in the Time-resolved Spectra of GRB 190114C, GRB 130427A, GRB 160509A, and GRB 160625B

[Shortened] CONTEXT: [...] AIMS: To identify and verify the BdHNe I properties in the additional sources GRB160509A, GRB160625B and GRB1340427A, and compare and contrast the results with the ones of a BdHN II source GRB180728A. We have also identified in all four sources, following the analysis GRB 130427A in the companion paper, the GeV radiation during and following the UPE phase. Also in all the four sources, we describe the spectral properties of their afterglow emission, including the mass estimate of the $\nu$NS, following the results presented in the companion paper. METHODS: [...] RESULTS: The results of the spectral analysis have validated the common properties in all BdHNe I: the three Episodes as well as the self-similar structures and the associated power-laws in the UPE phase. The profound similarities of the results have made a significant step forward in the taxonomy of GRBs and in evidencing a standard composition of the BdHN I. This opens the opportunity of a vaster inquire of the astrophysical nature of their components in the population synthesis approach: e.g., the BH formation in all BdHN I occurs due to accretion of the SN ejecta in a tight binary system with a neutron star companion which reaches its critical mass, leading to the formation of the BH. The SN-rise in all five BdHNe are compare and contrasted. CONCLUSIONS: The most far reaching discovery of self-similarities and power-laws here extensively confirmed, thanks also to the conclusions presented in the companion papers, leads to the existence of a discrete quantized repetitive polarized emission, both in the GeV and MeV observed by {\it Fermi}-GBM and {\it Fermi}-LAT, on a timescale as short as $10^{-14}$s. These results open new paths in the discovery of fundamental physical laws.


Introduction
As pointed out in the well documented book by Bing Zhang (Zhang 2018), the traditional approach in the spectroscopic data analysis of BATSE on-boarded the Compton Gamma-Ray Observatory (CGRO) (Preece et al. 2000) has typically addressed a time-integrated spectral analysis on the entire duration of T 90 and in finding commonalities in all GRBs. This approach has been continued all the way to the current Fermi-GBM observations and the observations of the BAT instrument on-board the Niels Gehrels Swift Observatory (see e.g., Abdo et al. 2009;Hamburg et al. 2019, by the Fermi team). The time-integrated spectrum has been traditionally fitted by a smoothly-connected, broken power-law function, named the "Band" function (Band et al. 1993). The Band function is based on four parameters whose values vary from source to source without reaching uni-Article number, page 1 of 19 arXiv:1910.12615v1 [astro-ph.HE] 28 Oct 2019 versal values. A complementary spectral analysis limited to the brightest time bin has been addressed by fitting with power-laws, smoothly broken power-laws, Comptonized and Band models (Gruber et al. 2014). Izzo et al. (2012); Ruffini et al. (2014) started a timeresolved spectral analysis approach with strong sources (e.g. GRB 090618). Since 2018, this approach has been improved starting with GRB 130427A  and other 16 GRBs (Ruffini et al. 2018c). This approach is now been adopted in GRB 190114C .
We have correspondingly defined our priorities: 1) to address only the brightest GRBs observed by Fermi-GBM, Fermi-LAT as well as by the Niels Gehrels Swfit Observatory, so addressing a more limited number of sources with high significance S and in a much wider range of spectral energies; 2) in view of the strongest significance S , to identify Episodes which present specific spectral structures and determine the duration ∆T of each Episode in the source rest-frame; 3) to perform an even more detailed time-resolved spectral analysis on ever decreasing time intervals, within the total duration ∆T , which has led to identify the presence of self-similar structures and associated power-laws. We have determined new statistical significant spectral distributions and evaluated the corresponding luminosity in the cosmological rest-frame.
On January 15, 2019, we indicated that GRB 190114C, discovered by Fermi-GBM on January 14, 2019 (Hamburg et al. 2019), with a redshift z = 0.424 observed by NOT (Selsing et al. 2019), had to be identified as a BdHN (Ruffini et al. 2019c). As a BdHN, see Fig. 1, within 18.8 ± 3.7 days, a SN should be expected to appear in the same location of the GRB. After an extended campaign involving tens of observatories worldwide, the expectation of the optical SN signal was confirmed (Melandri et al. 2019a;Wang et al. 2019b). This success and the detection of TeV radiation by MAGIC (Mirzoyan et al. 2019) make GRB 190114C one of the best example of multi-messenger astronomy.
For the first time, all the BdHN phases (see Fig. 1) have been observed in GRB 190114C (see companion paper Ruffini et al. 2019a). The main purpose of this article is to verity that the results obtained in GRB 190114C are not an isolated case, but on the contrary, they are verified to exist to an equal level of confidence in the other BdHN I: GRB 130427A, GRB 160509A, and GRB 160625B. Particular attention is given to the accuracy of the spectral analysis to identify the above three Episodes, as well as the much more complex iterative statistical analysis on the UPE to identify the self-similarities and the associated powerlaws.
In all these sources, starting from the GBM trigger and a well determined redshift, we have progressed to identify: in Episode 1, the precursor including the first appearance of the SN (the "SN-rise") and the accretion of the SN ejecta onto the companion neutron star (NS). For GRB 160625B, see Fig. 3 (upper left panel); for GRB 160509A, see Fig. 7 (upper left panel); for GRB 130427A, see Fig. 11 (upper left panel). In Episode 2, the moment of formation of the BH, the simultaneous onset of the GeV emission and the onset of of the UPE phase with its characteristic cutoff power-law plus blackbody spectra observed by Fermi; For GRB 160625B, see Fig. 4; for GRB 160509A, see Fig. 8. In Episode 3, the X-ray emission from the "cavity" recently modeled in the companion paper Ruffini et al. (2019b). For GRB 160625B, see  .
In an unprecedented detailed spectral analysis performed on ever decreasing time steps, we have here observed the selfsimilarities and power-laws in the GBM emission of the UPE phase in GRB 160509A and GRB 160625B and reported in the 4 and Fig. 8 . The details of the self similarities spectral features in the Table 2 and 4 , as well as the details of the numerical values of the fitting parameters for each GRB. We have extended the validity of of the corresponding results obtained for GRB 190114C in Ruffini et al. (2019a). In the case of GRB130427A this analysis has not been possible in view of the piled up data in the UPE phase. For all sources the spectral properties of the SN rise have been obtained quote Tables 1, 3, and 5, as well as the spectral properties of the GeV emission following the UPE phase and the spectral properties of the "cavity" and of their afterglows.
The results here presented indicate a clear progress in ascertaining the taxonomy of a standard BdHN I the comparison and contrast with the BdHN II (Wang et al. 2019b) has also greatly contributed in clarifying the role of the formation and/or absence of formation of the BH in BdHNe. There are three main new direction of research open: a) to submit to additional analysis each astrophysical component of each BdHN in the context of the different physical conditions characterising each Episode and harvest the physical novelties made possible by the observations of these previously unexplored regimes, not yet reproducible in earth based experimental facilities; b) to insert the BdHN observations within the larger context of population synthesis analysis, see e.g. discussion in Fryer et al. (2015) and c) address the micro-physical and physical origin of the self-similarities and the associated power-laws whose existence has been firmly confirmed in this article. In Ruffini et al. (2018b, in press), we have addressed the nature of the "inner engine" of BdHNe originating the structure of self-similarity: a Kerr BH embedded in a magnetic field aligned with its rotation axis and surrounded a very low density electron-ion plasma , and we have shown that the rotational energy extraction from the BH leads to the discrete and quantized MeV and GeV radiations. These results as well as the ones here presented open unprecedented opportunities for reaching new results in the determination of the physical laws in nature.
The structure of this article is as follows: In section 2 it is presented the detailed time-resolved data analysis procedure. We have fully considered the spectral contribution from thermal components. Our approach of the spectral analysis is based on fitting Bayesian models by using MCMC technique, which is superior compared to previous techniques (e.g., Frequency method). Based on the Bayesian analysis and MCMC technique, the preferred model is the one with the lowest DIC score.
In section 3, we derive our complete spectral analysis for all the episodes of GRB 160625B.
In section 4, we derive complete spectral analysis for GRB 160509A.
In section 5, the corresponding analysis for GRB 130427A, which is the only case that the UPE analysis is hampered by the pile-up problem.
In section 6, we recall the result of the BdHNe II GRB 180720A.
In section 7, we summarize the results on the analysis of the SN-rise of BdHNe I and II, and present the implications of these results in the physical and astrophysical scenario of BdHNe.
In section 8, we draw the general conclusions of this work.  Fryer et al. 2015), namely a carbon-oxygen star (CO core ) forming a tight (orbital period ∼ 5 min) binary with a neutron star (NS) companion. The CO core of mass 9-10 M undergoes core-collapse forming at its center a newborn NS (hereafter νNS) and, at the same time, ejecting the outermost layers in a type Ic SN explosion. The ejecta expand and their first observational appearance is what we call the "SN-rise". The ejecta reach the NS companion triggering a hypercritical accretion process onto it also thanks to a copious neutrino-antineutrino emission (Fryer et al. 2014;Becerra et al. 2018). Numerical simulations have shown that the NS companion, by accretion, reaches the critical mass for gravitational collapse, hence forming a BH. This was first shown by two-dimensional simulations in Becerra et al. (2015) and by threedimensional ones, first in Becerra et al. (2016), and more recently in improved smoothed-particle-hydrodynamics (SPH) simulations in Becerra et al. (2019), from which the simulated images shown in this figure have been taken. The fundamental contribution of these simulations has been to provide a visualisation of the SN morphology that is modified from its original sphericity. A low-density region, a "cavity", is carved by the NS companion and, once its collapses, it is further depleted to a density as low as ∼ 10 −14 g cm −3 by the BH formation process (see Ruffini et al. 2019b). The newborn Kerr BH, embedded in the magnetic field inherited from the collapsed NS, and aligned with the BH rotation axis, and surrounded by the low-density ionized plasma of the cavity, is what conform the "inner engine" of the GRB, see Ruffini et al. (2018b); . The "inner engine" leads to MeV emission due to the e + e − plasma created by vacuum polarization in the ultra-relativistic prompt emission (UPE) and to GeV emission by the synchrotron emission of accelerated electrons moving in the magnetic field. Details of these quantum and classic electrodynamics processes driven by the "inner engine" are given in companion papers; see Ruffini, Moradi, et al. (2019, submitted) and Ruffini, et al. (2019, submitted). The portion of the e + e − plasma that enters the high density region of the ejecta produce X-ray flares observed in the early afterglow (Ruffini et al. 2018c). The synchrotron emission by relativistic particles injected from the νNS into the expanding ejecta in the νNS magnetic field, explain the X-ray afterglow and its power-law luminosity (Ruffini et al. 2018a;Wang et al. 2019b). Finally, the optical emission from the ejecta due to the traditional nickel decay is observed in the optical bands few days after the GRB trigger.

BATSE Data
In the BATSE era, the observation of GRB spectra is limited to 25-1800 keV (Preece et al. 2000). Band et al. (1993) studied 54 BATSE GRBs. It was there analyzed the time-integrated spectra without addressing the specific physical process operating at any given moment. The software used to perform the spectral analysis was the Burst Spectral Analysis Software (BSAS; Schaefer 1991), and the basic fitting algorithm is based on CURFIT of Bevington (1969) (see page 237 therein), which finds the optimum spectral parameters by minimizing χ 2 . For the fitting of the spectra, they proposed and adopted the Band function (Band et al. 1993), which they found well described all the BATSE spectra. The parameters of the Band function are the low-energy spectral index α, the high-energy spectral index β, and the peak energy E p . They pointed out that these parameters vary for each GRB without reaching universal values, so they were not able to open up the GRB taxonomy 1 .

Fermi Data
The Fermi satellite, launched in 2008, provides a wider observational window in energy (Fermi-GBM: 8 keV to 40 MeV, Fermi-LAT: 100 MeV to 100 GeV), as well as a higher time resolution (as low as 5 µs for time-tagged event data) (Meegan et al. 2009). Gruber et al. (2014) presented the catalog of spectral analysis of GRBs by Fermi-GBM during its first four years of operation. They studied two types of spectra, the time-integrated spectrum and the spectrum of the brightest time bin. The software of RMfit (version 4.0rcl) was employed, which applies a modified, forward-folding Levenberg-Marquardt algorithm for spectral fitting. Four different spectral models were adopted: Band, Comptonized cut-off power-law (CPL), power-law (PL), and smoothly broken power-law models (SBPL). The PL and CPL models are preferred for most GRBs, the popularity of the simple PL model was interpreted as an observational effect.
In our approach since 2018, we implement the data from the Neil Gehrels Swift and Fermi satellite, our priority of having bright GRBs has already been stated in the introduction.

Spectral Models
There are several basic spectral components have been proposal in previous literature (e.g., Zhang et al. 2011). The observed GRB spectrum in keV-Mev band usually can be fitted by a nonthermal component, namely, the Band (or CPL) function (Band et al. 1993). The Band function defined as has two power-law photon indices: the low-energy power-law photon spectral index α (typically ∼ -1.0), and the high-energy 1 Quote from Band et al. (1993): "BATSE spectra are well described by this spectral form, but that α, β and E 0 all vary; there are no universal values. Such diversity must be addressed by physical models of the burst process ... we do not find any striking characteristics upon which to base a classification taxonomy." power-law photon spectral index β (typically ∼ -2.2), they are connected at the peak energy E p (typically ∼ 300 keV) in the νF ν space (e. g., Preece et al. 2000;Kaneko et al. 2006), A is the normalization factor at 100 keV in units of ph cm −2 keV −1 s −1 , E piv is the pivot energy fixed at 100 keV, the break energy E 0 in units of keV. For the UPE phase, we mainly adopt the CPL model, or the so-called Comptonized model (COMP), which is given by where A is the normalization factor at 100 keV in units of ph cm −2 keV −1 s −1 , E piv is the pivot energy fixed at 100 keV, α is the low-energy power-law photon spectral index, and E 0 is the break energy in units of keV. Some bursts have an additional thermal component, and generally fitted with Planck blackbody (BB) function. The Planck function, which is given by where A(t) is the normalization, k is the Boltzmann constant and kT (t) the blackbody temperature.
For the high-energy Fermi-LAT emission, the best-fit spectral model is usually a power-law model (e.g., Abdo et al. 2010;Zhang et al. 2011;Ajello et al. 2019) in the 0.1-100 GeV energy band, i.e., where A is the normalization, and Γ is the power-law index.
In the spectral fitting for the MeV-UPE phase, we adopt a Bayesian analysis and model comparison using ∆DIC value (e.g., Burgess et al. 2017;Li 2019a,b,c). For the GeV emission, a Maximum Likelihood Estimate (MLE) analysis is used to obtain the best fitting (e.g., Goldstein et al. 2012;Ackermann et al. 2013;Gruber et al. 2014;Narayana Bhat et al. 2016;Ajello et al. 2019;Li et al. 2019).

Bayesian Analysis for Fermi-GBM Data
The temporal and spectral analysis of Fermi-GBM data is applied by the Bayesian approach package, namely, the Multi-Mission Maximum Likelihood Framework (3ML, Vianello et al. 2015). The GBM (Meegan et al. 2009) carries 14 detectors that includes 12 sodium iodide (NaI, 8 keV-1 MeV) and 2 bismuth germinate (BGO, 200 keV-40 Mev) scintillation detectors. The pre-source and the post-source data are used to fit the background by a 0-4 order polynomial function. The time interval of source is selected longer than the duration of bursts (T 90 ), in order to cover the entire background subtracted emission. During the fitting procedure, the likelihood-based statistics, the so-called Pgstat is used, given by a Poisson . We use the typical spectral parameters from Fermi-GBM catalogue as the informative priors: α ∼ N(µ = −1., σ = 0.5); E c ∼ N(µ = 200, σ = 300); β ∼ N(µ = −2.2, σ = 0.5). Each time we perform 20 chains, each chain includes 10, 000 time iterations. The final value and its uncertainty (68% (1σ) Bayesian credible level) are calculated from the last 80% of the iterations. In this paper, we adopt the deviance information criterion (DIC) to select the best one from two different models, defined as DIC = −2 log[p(data |θ)] + 2p DIC , whereθ is the posterior mean of the parameters, and p DIC is the effective number of parameters. The preferred model is the model with the lowest DIC score. We define ∆DIC = DIC(CPL + BB) − DIC(CPL), for instance, if ∆DIC is negative, indicating the CPL+BB is better than CPL. These methods have be applied in each Episode. The best-fit parameters for each spectrum (α, E c ), along with its time interval, ∆DIC, blackbody temperature kT , blackbody flux (F BB ), total flux (F total ), thermal to total flux ratio, and the total energy are summarized in Table 2 for GRB 160625B and in Table 4 for GRB 160509A.

Calculation of Luminosity and Energetics
The observed flux, Φ(E 1 , E 2 , z), integrated between the minimum energy E 1 and the maximum energy In principle, for different models and different energy bands the values of E 1 , E 2 and f obs would be different. For instance for GeV radiation E 1 = 0.1 GeV and E 2 = 100 GeV and f obs = f PL is obtained from Eq. (5) with typical value of Γ ≈ −2.5 (Ajello et al. 2019).
We adopt a flat FLRW universe model with Ω Λ = 0.714, Ω M = 0.286 and H 0 = 69.6 km s −1 Mpc −1 (Bennett et al. 2014;Planck Collaboration et al. 2016), and then the luminosity distance is given by (Weinberg 1972) The isotropic radiated luminosity is where d L is the luminosity distance, and z is the redshift. The observed fluence S is given by where ∆T i is the duration of the time interval in which the analysis is made, see Ajello et al. (2019) for details. The isotropic radiated energy which is assumed to be isotropically radiated is defined as

GRB 160625B
On 25  Swfit-XRT starts the observation at late time (> 10 4 s), a powerlaw behaviour with decaying index ∼ −1.25 (Melandri et al. 2016). GRB 160625B is one of the most energetic GRBs with an isotropic energy ≈ 3 × 10 54 erg Zhang et al. 2018). The redshift z = 1.406 is reported in Xu et al. (2016). GRB 160625B is a luminous GRB with the clear detected polarisation (Troja 2017). There is no supernova confirmation due to its high redshift; z > 1 (Woosley & Bloom 2006). GRB 160625B has been extensively analyzed in Troja et al. (2017) and Zhang et al. (2018). Both papers define the early emission as three episodes: a short precursor (G1), a main burst (G2), and a long lasting tail (G3).
Based on the temporal and spectral analysis, we confirm that the gamma-ray light curve of GRB 160625B has three different episodes, shown in Fig. 2 (see also Table 1). Three different physical episodes have been identified in the keV-MeV energy range (see Fig. 2, Fig. 3, and Table 1): (1) SN-rise, the timeinterval ranges from t rf = 0.00 s to t rf = 0.83 s. (2) UPE phase, the time-interval ranges from t rf = 77.72 s to t rf = 87.70 s.
(3) Cavity, the time-interval ranges from t rf = 87.70 s to t rf = 92.27 s.
In a BdHN I, the "inner engine" starts at the moment of formation of the BH, accelerating charged particles that radiate photons in a wide energy band, generating the UPE phase and the GeV photons. The onset of the UPE phase is indicated by the appearance of the thermal component since the plasma is originally optically thick. Since the count rate of GeV photons observed in the onset phase is a few per second, it requires to have the discrepancy of at most fractions of a second between the observed starting time of the UPE and the GeV. Indeed, as for GRB 160625B, the starting time of its thermal emission is just 0.38 s ahead of the observational time of the first GeV photon, which for the above reasons can be considered temporally coincident. This time coincidence is also observed in the other BdHNe I studied in this article.
-SN-rise. Figure 3 (upper left panel) shows the fit of the SN-rise spectrum during its rest-frame time interval of occurrence, i.e. from 0 to t rf 0.83 s. It is best fitted by PL+BB model with temperature 17.5 keV and power-law index −2.0.
-UPE phase. Similarly to GRB 190114C, we also find a self-similarity in the UPE phase for GRB 160625B after carrying out the detailed time-resolved spectral analysis, with a cutoff powerlaw + blackbody (CPL+BB) model, for five successive iteration process on shorter and shorter time scales (expressed in the laboratory and in the rest frame). For the first iteration, Fig. 4 (first layer) shows the best-fit of the spectrum of the UPE entire duration from t rf = 77.72 s to t rf = 87.70 s.
We then divide the rest-frame time interval in half and perform again the same spectral analysis for the two intervals, each of 4.99s, namely [77.72s-82.71s] and [82.71s-87.70s], obtaining the results shown in Fig. 4 (second layer). In the third iteration, we divide each of these half intervals again in half. We continue this procedure up to five iterations, i.e up to dividing the UPE in 16 time sub-intervals. For each iterative step, we give the duration and the spectral parameters of CPL+BB model, including: the low-energy photon index α, the peak energy E c , the BB temperature kT (k is the Stefan-Boltzmann constant), the model comparison parameter (DIC), the BB flux, the total flux, the BB to total flux ratio, and the total energy. The results are summarized in Fig. 4 and Table 2, which confirm the validity, also in GRB 160625B, of the self-similar structure first discovered in GRB 190114C.
A&A proofs: manuscript no. ms_arxiv  (Xu et al. 2016). The light-curve consists of two clear spikes, the isotropic energy in the first one is (1.09 ± 0.20) × 10 52 erg. The total energy is ≈ 3 × 10 54 erg. Lower panel: the rest-frame time and energy of Fermi-LAT photons in the energy band 0.1-100 GeV. The first photon of the GeV emission occurs at t rf = 78.1 s. The onset of the GeV radiation coincides with the onset of the UPE. The detailed information for each episode (SN-rise, UPE phase, Cavity, GeV, and afterglow emission), see Section 3 and Table 1, which include the typical starting time, the duration, the isotropic energy, and the preferred model. Table 1. The Episodes of GRB 160625B, including the starting time, the duration, the energy (isotropic), the preferred spectral model, and the references. For the starting time of GeV emission, we take the time of the first GeV photon form the BH. The GeV emission may last for a very long duration, but the observational time is limited due to Fermi-LAT is not capable to resolve the late-time low flux emission, therefore the ending time of GeV observation in the table gives the lower limit of the ending time of GeV emission. The starting time of X-ray afterglow in the table is taken from the starting time of Swift-XRT. The energy in the afterglow is integrated from 10 2 s to 10 6 s. All times are given in the rest frame.

Episode
Starting    Table 1 that includes, for each Episode, the starting time, the duration, the isotropic energy, and the model that best fits the spectrum. Upper left: The spectrum of the SN-rise from 0 s to ≈ 2.0 s (t rf ≈ 0.83 s). The spectrum is fitted by a blackbody of temperature 17.5 keV (in the observer's frame) plus a power-law of index −2.0. Upper right: The cavity spectra, from ≈ 211 s (t rf = 87.70 s) to ≈ 222 s (t rf = 92.27 s), are well fitted by a CPL, with the photon index α is −1.67 and the cutoff energy 251 keV in the observer's frame. Lower left: Fermi-LAT rest-frame luminosity in the 100 MeV to 100 GeV energy band (the UPE region is marked with the grey shadow). Lower right: k-corrected X-ray afterglow luminosity observed by Swift-XRT in the 0.3-10 keV band, as a function of the rest-frame time. It is best fitted by a power-law with index 1.319 ± 0.021. Figure 5 shows the luminosity of the Fermi-GBM as a function of the rest-frame time, derived from the fifth iteration (see Table 2). We also show the corresponding time-evolution of the rest-frame temperature.
Article number, page 7 of 19 A&A proofs: manuscript no. ms_arxiv -Cavity. Figure 3 (upper right panel) shows the spectrum of cavity for GRB 160625B, from t rf = 87.70 s to t rf = 92.27 s, can be well fitted by a featureless CPL model, with the photon index α = 0.95 and the cutoff energy is at 239 keV.
-GeV emission. Figure 3 (lower left panel) shows the luminosity of the GeV emission in the rest-frame as a function of the rest-frame time.
-Afterglow. Figure 3 (lower right panel) shows the (kcorrected) afterglow luminosity (Swift/XRT data) in the restframe as a function of rest-frame time, and obtained a best fit parameters with the power-law index of 1.319 ± 0.021.

GRB 160509A
GRB 160509A was observed by the Fermi satellite on May 9, 2016, at 08:59:04.36 UT (Longo et al. 2016). It was a strong source of GeV photons detected by Fermi-LAT, including a photon of 52 GeV arrived at 77 s, and another one of 29 GeV, at ∼ 70 ks (Laskar et al. 2016). Swift has a late-time follow-up, with a total exposure time of 1700 s starting from 7278 s (Kangas et al. 2019). The redshift of 1.17 is measured by Gemini North telescope (Tanvir et al. 2016), inferring a high isotropic energy of 1.06 × 10 54 erg (Tam et al. 2017).
Pak-Hin Thomas Tam and collaborators (Tam et al. 2017) analyzed in great detail the bright multi-peaked pulse from −10 to 30 s, and a weaker emission period from 280 to 420 s. In Fig. 6, we show the result of our highly time-resolved analysis applied to GRB 160509A, which further extend the result of Tam et al. (2017).
Based on the temporal and spectral analysis, three different physical processes are identified in the keV-MeV energy range (see Fig. 6, Fig. 7, and Table 3): (1). SN-rise, the time interval ranges from t rf = 0.92 s to t rf = 1.84 s. (2). UPE phase, the time interval ranges from t rf = 4.84 s to t rf = 8.53 s. (3). Cavity, the time interval ranges from t rf = 10.14 s to t rf = 13.82 s.
-SN-rise. Figure 7 (upper left panel) shows the SN-rise and the spectral fitting of the cavity emission during its restframe time interval of occurrence. The spectrum of the SNrise of GRB 160509A are best fitted by CPL+BB model, from t rf 0.92s to t rf 1.84s. The spectrum contains a BB component of temperature 25.61 keV and a photon index α of −1.22, and E c = 1769.76 keV.
-UPE phase. We perform the corresponding time-resolved spectral analysis from which we can see that the self-similarity first discovered in GRB 190114C is confirmed in the case of GRB 160509A. For the first iteration, we present the best-fit of the spectrum of the UPE entire duration from t rf = 4.84 s to t rf = 8.53 s (see Fig. 8, first layer).
We then divide the rest-frame time interval in half and perform again the same spectral analysis for the two 1.85 second interval, namely [4.84s-6.68s] and [6.68s-8.53s], obtaining the results shown in Fig. 8. Iteration 3: we then divide each of these half intervals again in half, i.e., ∆t rf = 0.92s corresponding to [4.84s-5.76s], [5.76s-6.68s], [6.68s-7.60s] and [7.60s-8.53s] and redo the previous spectral analysis obtaining the results still in Fig. 8 -Cavity. Figure 7 (upper right panel) shows the spectral fitting of the cavity emission during the rest-frame time interval of its occurrence, i.e. from t rf = 10.14 s to 13.82 s. The best fit of the spectrum is a CPL model, with a photon index α is −1.20 and the cutoff energy is 314 keV.
-GeV emission. Figure 7 (lower left panel) shows the luminosity of the GeV emission and the luminosity in the afterglow are given as a function of the rest-frame time.
-Afterglow. Figure 7 (lower right panel) shows the (kcorrected) rest-frame afterglow luminosity (Swift/XRT data) as a function of rest-frame time, and obtained a best fit parameters with the power-law index of −1.259 ± 0.025.
The Fermi-GBM count rate of GRB 130427A is shown in Fig. 10. During the UPE phase the event count rate of n9 and n10 of Fermi-GBM surpasses ∼ 8 × 10 4 counts per second in the prompt radiation between rest-frame times T 0 + 3.4 s and T 0 + 8.6 s. The GRB is there affected by pile-up, which significantly deforms the spectrum; details in Ackermann et al. (2014); Ruffini et al. (2015), only the data between t rf = 0.0 and t rf = 1.49 s can be used for a spectral analysis in the prompt phase. As shown in Fig. 10 (see also Table 5) clearly identified parts are: -SN-rise. Figure 11 (upper left panel) shows the clear identification of the SN-rise, as also reported in Fig. 10. The spectrum of the SN-rise of GRB 130427A are best fitted by CPL+BB model, from 0.0s (t rf 0.0s) to 0.65 s (t rf 0.49s). The spectrum contains a BB component of temperature 42.63 keV and a photon index α = −0.58, and E c = 547.59 keV.
-Cavity. Figure 11 (upper right panel) shows the featureless spectrum of the cavity emission of GRB 130427A, fitted by a CPL model, from ≈ 15 s (t rf = 11.19 s) to ≈ 25.5 s (t rf = 19.03 s) fitted by a CPL, with the photon index α = −1.52 and the cutoff energy is at 496.13 keV.
-GeV emission. Figure 11 (lower left panel) shows the restframe luminosity of the GeV emission as a function of the restframe time.
-Afterglow. Figure 11 (lower right panel) shows the (kcorrected) luminosity of the afterglow (Swift/XRT data) as a function of the rest-frame time, and obtained as best-fit a powerlaw index of −1.276 ± 0.002.   We have converted to have their corresponding value in the rest-frame, see Table 2, where rest-frame time in column 2, and rest-frame temperature in column 6. Table 2. Results of the time-resolved spectral fits of GRB 160625B (CPL+BB model) from the t rf = 77.72 s to t rf = 87.70 s. This table reports: the time intervals both in rest-frame and observer frame, the significance (S ) for each time interval, the power-law index, cut-off energy, temperature, ∆DIC, BB flux, total flux, the BB to total flux ratio, F BB /F tot and finally the isotropic energy. To select the best model from two different given models, we adopt the deviance information criterion (DIC), defined as DIC=-2log[p(data|θ)]+2p DIC , whereθ is the posterior mean of the parameters, and p DIC is the effective number of parameters. The preferred model is the model with the lowest DIC score. Here we define ∆DIC=(CPL+BB)-CPL, if ∆DIC is negative, indicating the CPL+BB is better. After comparing the DIC, we find the CPL+BB model is the preferred model than the CPL and other model. The ∆DIC scores are reported in column 6. S r =-3.39±1.08  Table 2 are used to apply the k-correction and plot the rest-frame luminosity as a function of rest-frame time. The power-law index of the luminosity is at −19.89 ± 4.05. For more information about GeV luminosity behavior see (Wang et al. 2019a) and the companion paper (Ruffini, Moradi et al, 2019, submitted). Right: Corresponding rest-frame temperature of the UPE as a function of the rest-frame time. after t r f = 10.14 s, staring at t r f = 10.14 s and ending at t r f = 13.82 s. The redshift for GRB 160509A is 1.17 (Tanvir et al. 2016). The light-curve consists of two spikes, the isotropic energy in the first small one is ∼ 1.47 × 10 52 erg. The total energy is 1.06 × 10 54 erg (Tam et al. 2017). Lower panel: the energy and time of each Fermi-LAT photon of energy > 100 MeV. The first GeV photon occurs at 4.84 s in the rest-frame. The onset of the GeV radiation exactly coincides with the onset of the UPE. The detailed information for each episode (SN-rise, UPE phase, Cavity, GeV, and Afterglow) see Section 4 and Table 3, which include the starting time, the duration, the isotropic energy, and the preferred model.

BdHN II
A&A proofs: manuscript no. ms_arxiv    We have converted to have their corresponding value in the rest-frame, see Table 4, where rest-frame time in column 2, rest-frame cutoff energy in column 5 and rest-frame temperature in column 6.
Article number, page 13 of 19 A&A proofs: manuscript no. ms_arxiv Table 4. Results of the time-resolved spectral fits of GRB 160509A (CPL+BB model) from the t rf = 4.84 s to t rf = 8.53 s. The definitions of parameters are the same as in table 2.   Table 4 are used to apply the k-correction and measuring the luminosity as a function of time. The power law index of 0.51 ± 0.35 of the luminosity is similar to the one obtained in the GeV emission luminosity after the UPE phase with index of −0.22 ± 0.23. For more information about GeV luminosity behavior see (Wang et al. 2019a) and the companion paper (Ruffini, Moradi et al, 2019, submitted). Right: Time evolution of the rest-frame temperature of the UPE as derived from the fifth iteration with 16 sub-intervals, is reported in Table 4  A&A proofs: manuscript no. ms_arxiv   Table 5 that includes, for each Episode, the starting time, the duration, the isotropic energy, and the model that best fits the spectrum. Upper left: SN-rise spectrum, well fitted by a CPL+BB model, from 0 to 0.65s (t rf 0.49s); the spectral index α is -0.58, cutoff energy E c is 547.59 keV, and the BB temperature is 42.63 keV in the observer's frame. Upper right: Featureless spectrum of the cavity emission, from ≈ 15 s (t rf = 11.19 s) to ≈ 25.5 s (t rf = 19.03 s) fitted by a CPL, with the photon index α is −1.52 and the cutoff energy is at 496.13 keV in the observer's frame. Lower left: Fermi-LAT rest-frame light-curve in the 100 MeV to 100 GeV energy range. The UPE region is marked with the grey shadow. Lower right: k-corrected X-ray afterglow luminosity observed by Swift-XRT in the 0.3-10 keV energy range, as a function of the rest-frame time. It is best fitted by a power-law with index 1.276 ± 0.002. redshift z = 0.117 detected by VLT/X-shooter (Rossi 2018). On July 28, we made a prediction of the SN appearance in ∼ 15 days (Ruffini 2018;Wang et al. 2019b), and indeed the SN optical peak was confirmed then (Izzo 2018;Selsing 2018).
This GRB is composed of two pulses, see Fig. 12 and Table 6: -First pulse as SN-rise. The first spike, the precursor, shows a power-law spectrum with a power-law index −2.31 ± 0.08 in its 2.75 s duration. The averaged luminosity is 3.24 +0.78 −0.55 × 10 49 erg s −1 , and the integrated energy gives 7.98 +1.92 −1.34 × 10 49 erg in the energy range from 1 keV to 10 MeV. This energy emitted is in agreement the conversion of the SN-rise kinetic energy into electromagnetic emission.
-Second pulse as the hypercritical accretion of the SN ejecta onto the companion NS. This pulse starting from 8.72 s and lasts 13.82 s which contains 2.73 × 10 51 erg isotropic energy. The best-fit is which is a CPL+BB model of temperature ≈ 7 keV in the observer frame is shown in Fig. 12 Fig. 12. We identify the SN-rise from the CO core of a BdHN II in GRB 180728A (Wang et al. 2019b). This GRB is composed of two spikes. The first spike, the precursor, shows a power-law spectrum with a power-law index −2.31 ± 0.08 in its 2.75 s duration. The averaged luminosity is 3.24 +0.78 −0.55 × 10 49 erg s −1 , and the integrated energy gives 7.98 +1.92 −1.34 × 10 49 erg in the energy range from 1 keV to 10 MeV. This energy emitted is in agreement the conversion of the SN-rise kinetic energy into electromagnetic emission. We consider second pulse (prompt emission without self similarity) as due to the hypercritical accretion of the SN ejecta onto the companion NS, starting from 8.72 s and lasts 13.82 s. This pulse contains 2.73 × 10 51 erg isotropic energy. The best-fit is a CPL+BB model of temperature ≈ 7 keV in the observer frame. The BB component is interpreted as a matter outflow driven by the Rayleigh-Taylor convective instability developed in the accretion process. From the time of observing the SN-rise to the starting time of the hypercritical accretion, ∆t ≈ 10 s, it has been inferred a binary separation of ≈ 3 × 10 10 cm. From the binary separation, by angular momentum conservation, it has been inferred that the spin of the νNS left from the collapse of the CO core , is ≈ 2.5 ms (Wang et al. 2019b). This νNS powers the afterglow by dissipating its rotational energy. ponent is interpreted as a matter outflow driven by the Rayleigh-Taylor convective instability developed in the accretion process (see e.g. Izzo et al. 2012). From the time of observing the SNrise to the starting time of the hypercritical accretion, ∆t ≈ 10 s, it has been inferred a binary separation of ≈ 3 × 10 10 cm. The binary separation determines, by angular momentum conservation, the spin period of ≈ 2.5 ms of the νNS left from the collapse of CO core . This νNS powers the afterglow by dissipating its rotational energy (Wang et al. 2019b).

Discussion
In Table 7, we compare and contrast the duration, the fluxes, the energy, the temperature of the BB component associated with the SN-rise of the above BdHNe I and II, we give as well, for each GRB, the corresponding redshift and E iso . In the case of BdHN I, all of them have a similar SN-rise duration of nearly a second, consistent with the radius of the CO core of 10 10 cm, and energies of the order of 10 52 erg. These energies are much larger than the one we have here found in in the SN-rise of BdHNe II, ∼ 10 50 erg, comparable to the one of isolated SNe (see, e.g., Arnett 1982;Bethe 1990;Waxman & Katz 2017).

The SN-rise energetics of BdHN I
The larger energies of the SN-rise associated with BdHNe I here discovered can also be ascribed to a more energetic, rapidly rotating CO core . This can be the result of the binary nature of the progenitor with a short orbital period of the order of 4-5 min in which angular momentum transfer by tidal effects during the previous evolutionary stages has been at work very efficiently.
Let us estimate the rotational energy of the CO core assuming that the binary is tidally locked. In this case the CO core rotation period, P CO , equals the binary orbital period, P orb (see, e.g., Hurley et al. 2002), i.e.
which is related to the binary separation a orb and the total mass of the system M tot and G is the gravitational constant. Let us adopt a typical progenitor of a BdHN from Becerra et al. (2019): a CO core obtained from the evolution of a 30 M zero-age mainsequence (ZAMS) progenitor star, which have a total mass of M CO = 8.9 M and radius R CO = 7.83 × 10 9 cm, forming a binary with a NS companion of M NS = 2 M . As for the orbital period/separation, we constrain our systems by the condition that there is no Roche-lobe overflow at the moment of the supernova explosion of the CO core . The Roche lobe radius of the CO core can be estimated as (Eggleton 1983) R RL a orb = 0.49q 2/3 0.6q 2/3 + ln(1 + q 1/3 ) , where q = M CO /M NS . Therefore, the minimum orbital period of the binary, a orb,min , is obtained when R CO = R RL . For the above parameters, a orb,min ≈ 1.53 × 10 10 cm and correspondingly the minimum orbital period is P orb,min ≈ 5.23 min. The rotational energy for a CO core is where I CO is the CO core moment of inertia. So, adopting P CO = P orb,min (ω CO ≈ 0.03 rad s −1 ) and I CO ≈ (2/5)M CO R 2 CO , we obtain E rot,CO ≈ 8.7×10 49 erg. This is of course lower than the gravitational binding energy |W| ≈ (3/5)GM 2 CO /R CO ≈ 1.6 × 10 51 erg and lower than the internal thermal energy as from the virial theorem. If we adopt the CO core from the 25 M ZAMS progenitor (see Table 1 in Becerra et al. 2019), characterized by M CO = 6.85 M and R CO = 5.86 × 10 9 cm, and for the corresponding minimum orbital period P orb,min ≈ 4 min (ω CO ≈ 0.02 rad s −1 ), we obtain E rot,CO ≈ 6.3 × 10 49 erg.
Therefore, a much more energetic SN-rise can be the result of an exploding CO core which rotates at a much higher rotation rate with respect to the one set by tidal synchronization. In the above two examples, the rotational to gravitational energy ratio is E rot /|W| ≈ 0.05. However, from the stability point of view, it is known from the theory of Newtonian ellipsoids that secular axisymmetric instability sets in at E rot /|W| ≈ 0.14 and dynamical instability at E rot /|W| ≈ 0.25 (Chandrasekhar 1969).
Indeed, three-dimensional simulations of SN explosions confirm these stability limits and so explore SN explosions from pre-SN cores with high rotation rates of the order of 1 rad s −1 (see, e.g., Nakamura et al. 2014;Gilkis 2018;Fujisawa et al. 2019). These angular velocities are a factor 30-50 faster than the ones we have considered above. This implies that the rotational energy of the pre-SN core can be up to a factor 10 3 higher, namely E rot ∼ few × 10 52 erg.

The SN-rise energetics of BdHN II
In the case of BdHN II, the SN-rise has been shown to have a much smaller energy, 10 49 -10 50 erg. A similar case in the literature is represented by SN 2006aj, associated with GRB 060218 (Campana et al. 2006;Mirabal et al. 2006;Pian et al. 2006;Sollerman et al. 2006;Ferrero et al. 2006). The GRB 060218/SN 2006aj association was indeed interpreted in Becerra et al. (2016) as a BdHN II (at that time called "X-ray flash"). As we have mentioned, the energetics of these SN-rises are closer to the typical ones encountered in isolated SNe (see, e.g., Arnett 1982;Bethe 1990;Waxman & Katz 2017). This is consistent with the longer orbital periods of BdHNe II (Becerra et al. 2016) since, being farther apart, in the prior evolutionary stages binary interactions have been less effective in transferring angular momentum to the CO core . This explains why the SNe associated with BdHN II, even if they occur in a binary, are more similar to isolated SNe.
As a final remark, we recall that the occurrence of the SN is deduced from direct optical observations for GRB sources at z < 1, and for all cases the SN occurrence is also inferred, indirectly, by the observation of the afterglows. Indeed, the afterglow originates from the feedback of the emission of the νNS, originated in the SN event, into the expanding SN ejecta, given the proof of the SN occurrence (see Ruffini et al. 2018a;Wang et al. 2019b;, for details).

Conclusions
In a companion paper , we have introduced a novel time-resolved spectral analysis technique, adopting ever decreasing time steps, in the analysis of GRB 190114C. This has led to the discovery of the three Episodes and the self-similarity and power-laws in BdHN I.
In this paper, we have made a major effort of applying, such a time-resolved spectral analysis to BdHNe I: GRB 130427A, GRB 160509A, and GRB 160625B; see sections 3-5. We have proved that indeed, all the results obtained in GRB 190114C, far from being an exception, do characterize the physics of BdHNe I.
These results open new perspective of research: 1) to study the new physical process characterising each single Episode of a BdHN in the context of previously unexplored regimes: e.g. the analysis of the SN not following the traditional description as an isolated system and identifying their properties within a BdHN I and alternatively, in a BdHN II.
2) to insert the BdHN evolution in the framework a population synthesis analysis.
3) to address the new physical process underlying the existence of the observed self-similarities and power-laws revealing a discrete sequence of quantized events with quanta of 10 37 erg on new timescales of 10 −14 s; see the companion papers (Ruffini et al. 2018b;, Ruffini, Moradi, et al. (2019, submitted) and explore the vast new directions open to the identification of the fundamental new laws of our Universe.