Constraining the Type of Central Engine of GRBs with Swift Data

The central engine of gamma-ray bursts (GRBs) is poorly constrained. There exist two main candidates: a fast-rotating black hole and a rapidly spinning magnetar. Furthermore, the X-ray plateaus are widely accepted by the energy injection into the external shock. In this paper, we systematically analyze the \emph{Swift}/XRT light curves of 101 GRBs having plateau phases and known redshifts (before May 2017). Since a maximum energy budget ($\sim2\times10^{52}$ erg) exists for magnetars but not for a black hole, this provides a good clue to identify the type of GRB central engine. We calculate the isotropic kinetic energy $E_{\rm K,iso}$ and the isotropic X-ray energy release $E_{\rm X,iso}$ for individual GRB. We identify three categories based on how likely a black hole harbor at central engine: 'Gold' (9 out of 101, both $E_{\rm X,iso}$ and $E_{\rm K,iso}$ exceed the energy budget), 'Silver' (69 out of 101, $E_{\rm X,iso}$ less than the limit but $E_{\rm K,iso}$ is greater than the limit), and 'Bronze' (23 out of 101, the energies are not above the limit). We then derive and test the black hole parameters with the Blandford Znajek mechanism, and find that the observations of the black hole candidate ('Gold'+'Silver') samples are consistent with the expectations of the black hole model. Furthermore, we also test the magnetar candidate ('Bronze') sample with the magnetar model, and find that the magnetar surface magnetic led ($B_{p}$) and initial spin period ($P_{0}$) are consistent with the expectations of the magnetar model. Our analysis indicates that, if the magnetar wind is isotropic, a magnetar central engine is possible for 20\% of the analyzed GRBs. For most GRBs a black hole is most likely operating.


Introduction
One of the most fundamental, yet unanswered questions in gamma-ray burst (GRB) physics is the nature of the central engine. Observationally, they have to fulfill several requirements. First, they should be able to release a huge amount of energy (∼ 10 49 − 10 55 ergs). Second, indirect es-timations of the jet Lorentz factor ∼100, together with the compactness problem, require the central engine to be as pure as possible of baryons. Finally, the rapid variability hints towards an intermittent central engine (e.g., Fishman 1996;Bloom et al. 2001;Liang et al. 2010).
There are two main progenitor models for the central engine. On one hand, the magnetar model involves a rapidly spinning and highly magnetized neutron star. In this scenario, the rotational energy of the neutron star powers the GRB (e.g., Usov 1992;Thompson 1994;Dai & Lu 1998;Wheeler et al. 2000;Zhang & Mészáros 2001;Metzger et al. 2008;Metzger et al. 2011;Bucciantini et al. 2012;Margutti et al. 2013;Lü & Zhang 2014;Lü et al. 2015). The limited energy budget provided by the magnetar model might, however, be at strain with observed GRB energies (Cenko et al. 2010). On the other hand, a newly formed blackhole (hereinafter BH) surrounded by an accretion 1 arXiv:1712.09390v3 [astro-ph.HE] 29 Mar 2018 disk (e.g., Paczynski 1991;Narayan et al. 1992;MacFadyen & Woosley 1999) can produce bipolar jets through the Blandford and Znajek mechanism (Blandford & Znajek 1977) or through the production and polar anni-hilation of neutrinoantineutrino pairs (Liu et al. 2017). The rapid rotation of the central engine required by GRB models might, however, prevent the collapse to a black hole (Burrows et al. 2007;Dessart et al. 2008). Distinguishing between a black-hole and a neutron star as the central engine is of the utmost importance as it would help to better constrain the progenitor of GRBs, and their emission mechanisms.
The long time activity and decay of the central engine can leave footprints on the X-ray and optical afterglows and can thus help distinguishing between the proposed progenitors. For example, the optical re-brightening of GRB 970228 was interpreted as the result of energy injection by a millisecond pulsar (Dai & Lu 1998). In addition, Zhang & Mészáros (2001) showed that a shallow decay followed by a normal decay phase can be naturally obtained for the energy ejection by a millisecond magnetar. The widely accepted scenario for the X-ray plateaus (also called shallow decay) is the energy injection into the external shock, either matter-dominated or Poyntingflux (electron-positron wind) dominated injection. The simplest way to deal with the energy injection process is by assuming the central engine is a millisecond magnetar and the injected energy is carried out by a Poynting-flux/electron-positron wind, extracting the spin energy of the magnetar through the magnetic dipole radiation (MDR). Of particular importance, the X-ray plateau is often interpreted as energy injection into the external forward shock, which is an evidence for extended central engine activity. The energy can be provided either in terms of Poynting flux from the central magnetar, or by the deposition of kinetic energy into the forward shock by a stratified jet. Both models are successful to explain the normal decay phase at the end of the plateau. However, it happens that the plateau is followed by a sharp flux drop, as steep as F ∼ t −9 . Such a sharp decay is very difficult to explain using the forward shock model, since it requires the radiation to be released from a small radius. In this case, the plateau has to be of internal origin (as opposite to external when it can be explained by energy injection). Such plateaus are easily explained in the magnetar framework, as residual rotational and magnetic energy that can be extracted after the main extraction event. However, no model for the emission mechanism exists todate. The spindown luminosity of energy injection from a central engine would evolve with time as (Zhang & Mészáros 2001, the magnetar injection corresponds to q=0 for t< τ and q=2 for t> τ ), Here q is the luminosity injection index. For black holes, they can also give energy injection, but q depends on density profiles of the star (e.g., Kumar et al. 2008a,b). An additional criterion which could be used to differentiate between the magnetar and the blackhole model is the total energy emitted during the bursts. Indeed, heuristic arguments lead to a maximum energy that can be released by a magnetar, corresponding to the maximum rotational energy of 2 × 10 52 ergs (Lü & Zhang 2014). Such a limit does not exist for black-holes. The emitted radiation energy will only be some fraction of the rotational energy depending on the radiative efficiency (see e.g., Beniamini et al. 2017). More important, such a limit requires an assumption to be made on whether the energy emitted from a magnetar is isotropic or not. The reason is that the inferred, total energy is given by integrating the observed flux over all solid angles. The angular distribution of the flux must therefore be assumed. For a millisecond magnetar, the energy is carried out by a Poynting-flux/electron-positron wind, extracting the spin energy through the magnetic dipole radiation. The wind is expected to be launched in a near isotropic wind (e.g. Usov 1992). For instance, for the pulsar wind nebular of Crab, the MDR of a neutron star is nearly isotropic. Moreover, in the case of a neutron star merger as the origin of the magnetar central engine in short GRBs, the total energy is also expected to be emitted quasi-spherical (e.g. Metzger & Piro 2014;Fan & Xu 2006). However, if the magnetar is born at the center of a collapsar, the progenitor star could cause a collimation of the emitted energy, both during the prompt phase (Bucciantini et al. 2007(Bucciantini et al. , 2009 and the afterglow phase (Lü & Zhang 2014). However, Mazzali et al. (2014) argues for the magnetar wind to be close to isotropic within a collapsar based on the energetics of the associated supernovae. Only a small fraction of the total energy is then channeled into the GRB jet, while most of the total energy is responsible for the supernova explosion. Therefore, the isotropic injection has a low efficiency. However, the energy injection from a fast-rotating black hole is naturally considered to be a collimated injection through the Blandford & Znajek (BZ) mechanism (Blandford & Znajek 1977). The collimated injection is high efficiency, since a high fraction of the energy can be injected into the GRB jet.
In this paper, we use the assumption that the magnetar wind is quasi-spherical during the energy injection phase of the X-ray afterglow. We analyze the plateau afterglow of GRBs in order to investigate the plausibility of having a magnetar as a central engine. We extensively search the Swift/XRT data for GRBs in which the light curve has a plateau phase and calculate the X-ray energy released during this phase and the kinetic energy. We then investigate whether the data are consistent with the black hole central engine model or not.
The paper is organized as follow. The XRT data reduction and the samples selection are presented in §2. The physical parameters of the GRBs and the hypothetical BH/magnetar parameters derived in §3. The statistical analyzes of the external shock model parameters are presented in §4. The discussions and conclusions are presented in §5. Throughout the paper, a concordance cosmology with parameters H 0 = 71 kms −1 Mpc −1 , Ω M = 0.30, and Ω Λ = 0.70 is adopted.

Data Analysis
We restricted our study to GRBs observed by Swift/XRT from December 2004 and May 2017, which show a transition from a plateau (also called shallow decay) to a normal decay phase (for the case of canonical light curve), or transition to a very steep decay (for the case of "internal plateaus"). We first identify such bursts by visual inspection of the X-ray light curves. To derive the properties of the GRBs in the rest frame, we only select the bursts with known redshifts.

Sample Definition
Two possible forms of light curves may confuse our sample selections performed by visual inspection. First, we may take an intrinsic single powerlaw decay as a broken power-law decay. Second, we may also take a normal decay followed by a post-jet break phase as a plateau followed by a normal decay phase.
Therefore, to characterize a plateau, we perform a temporal fit to the X-ray light curves within a time interval (t s , t e ) for each burst. Here t s is the time when the plateau begins, and t e when the segment after the plateau break ends. We employ either a smoothly broken power law (BKPL), or a single power-law 1 (SPL) to fit the light curves, In the above equations, α, α 1 , and α 2 are the temporal slope, t b is the break time, F b = F 0 2 −1/ω is the flux of the break time, and ω describes the sharpness of the break 2 . We perform the fits to the data with a non-linear, least-square fitting using the Levenberg-Marquardt algorithm, and an IDL routine called mpfitfun.pro 3 (Markwardt 2009).
To exclude an intrinsic single power-law decay, we first apply above two models, respectively. We then choose the best model according to the following two principles. (i) We compare the χ 2 r value for both models, and chose the model having the χ 2 r closest to 1 4 . For example, the value of the χ 2 r of GRB 070529 for the BKPL model fitting is given 29/30 (close to 1) and for the SPL model fitting gave 68/33 (much larger than 1), we then adopt the BKPL for the best model for GRB 070529. (ii) In case both models have similar χ 2 r value (close to 1), we always chose the simplest model (SPL).
For instance, GRB 101219B has χ 2 r 16/18 (close to 1) for the BKPL model and 21/21 (close to 1) for the SPL model, in this case, SPL model is the best model for our choice. We then remove these cases that have the best fit of the SPL model 5 .
To properly select only plateau phases and remove the normal decay followed by a jet break phase, we only keep the temporal decays which are too fast to be explained by an external shock. For this, we use the closure relations to obtain an estimate of what should be the spectral index α(β) in the case of different spectral regimes (ISM/wind), here α 1 is the temporal decay index before the break time, Γ 1 is the photon spectral index, β 1 is the spectral index (β = Γ − 1), ν m is the minimum frequency and ν c is the cooling frequency.
Comparing the theoretical spectral indexes to the measured ones, we only keep the bursts which satisfy 6 α s 1 < α(β). Finally, 101 GRBs are selected. These are shown in Figs. 1, 2 and 3.

Type of Plateau
In order to identify the types of the plateau (internal or external), we compare the temporal (α 1 and α 2 ) and spectral (Γ 1 and Γ 2 ) indices 7 before and after the plateau break, and combine them through the closure relations. Since the external plateau is purely dynamical, no change of the spectral index is expected. As such we require Γ 1 =Γ 2 (within the errors) associated to the equality between the decay index measured and those obtained by closure relations from Γ 1 and Γ 2 . The pre-break slope α 1 (correspond to q=0) is given by: 5 We excluded GRB 080210, GRB 130514A, GRB 140430A from our samples. 6 With this selection, the burst removed from the study are: GRB 081208, GRB 081203A, GRB 080413B, GRB 090313, GRB 111123A, GRB 130427B, GRB 140213A. 7 Throughout the paper, the convention Fν ∝ ν −β t −α is adopted.
while the post-break slope α 2 (correspond to q = 1, Zhang & Mészáros 2001), according to Zhang et al. (2006), where β 2 is the spectral index during the normal decay phase, and p is the electron's spectral distribution index. These criteria gave 89 bursts in this sample. Fig.4 shows all the GRBs external plateau in the α 1 − α 2 plane, with three theoretically favored lines (Equation 5) plotted. Fig.5 shows the β − α plane based on the ISM or the wind model. The GRBs satisfying the "closure relations" are identified as our sample GRBs with colored data points. The internal plateau (Troja et al. 2007;Lyons et al. 2010;Rowlinson et al. 2010Rowlinson et al. , 2013 cases are selected based on Γ 1 different Γ 2 and α 2 > 2 + β. This sample is made of only 12 bursts. The bursts which cannot be put in those two categories are discarded from further analysis 8 .

X-ray Plateau and Afterglow Kinetic Energies
We derive the X-ray isotropic energy released (E X,iso ) during the plateau phase, by integrating the fluence from t s to t b : where k = (1 + z) β−1 , β is the spectral index, S X is the X-ray fluence 9 integrated over the plateau phase, D L is the luminosity distance and z is the redshift. Another important parameter is the isotropic kinetic energy (E K,iso ), which is measured from the afterglow flux during the normal decay. Here, E K,iso is calculated following the method discussed in Zhang et al. (2007b). Since E K,iso depends on the afterglow models (we remind it here briefly for completeness), it is important to identify the afterglow spectral regimes. We first use the observed quantity (α,β) in the normal decay phase to constrain/determine the afterglow spectral regimes of each individual burst . We find that 18 of 101 GRBs belong to spectral regime I (ν x > max(ν m , ν c )), and 83 out or 101 GRBs belong to spectral regime II (ν m < ν x < ν c ). Here ν x is the frequency of the X-ray band. In addition, 42 out of the 83 bursts satisfy the ISM model and 31 out of 83 satisfy the wind model.
We consider the relations obtained by Zhang et al. (2007b) to compute the kinetic energy of the external plateau sample using the normal decay segment. Since the normal decay phase is always observed at late time (t∼ 10 4 seconds) in our samples, we consider the slow cooling scenario (ν m < ν x < ν c ) (see also Beniamini et al. 2015). In the case (ν x > max(ν m , ν c )), the same express expression can be applied for both the wind-like and the constant interstellar medium since the kinetic energy does not depend on the medium profile in this spectral regime. E K,iso,52 = νF ν (ν = 10 18 Hz) 5.2 × 10 −14 ergs −1 cm −2 where νF ν (ν = 10 18 Hz) is the energy flux 10 at frequency 10 18 Hz in units of ergs −1 cm −2 , t d is the time in the observer frame in units of days, D 28 = D/10 28 , is the luminosity distance in units 11 of 10 28 cm, n is the number density in the constant ambient medium, A * is the parameter of stellar wind, Y is the Compton parameter, and e and B are the energy fraction assigned to electron and magnetic field, respectively. f p is a function of the 10 Note that since the νFν is not an observed quantity, and it can be converted to the observed flux by integral frequencies over X-ray energy range (0.3 keV-10 keV). We note that for a few cases, the photon index Γ X,2 = 2.0 during the normal decay phase, we adopt Γ X,2 2.01 (within the errors) make a approximate treatment in order to avoid the 1 − β = 0 (β = Γ − 1) for the calculations. 11 The convention Q = 10 x Qx is adopted in cgs units for all other parameters throughout the paper.
power law electron distribution p index 12 (Zhang et al. 2007b).
When ν m < ν < ν c , the kinetic energy depends on the type of external medium. For a wind medium, one has while for a constant ISM, one has E K,iso,52 = νF ν (ν = 10 18 Hz) 6.5 × 10 −14 ergs −1 cm −2 After deriving E K,iso at the break time t b , we next estimate the upper limit of E K,iso (t 0 ) by E K,iso (t b ). If t s > t 0 , the upper limit of E K,iso (t 0 ) can be estimated by The error on E K,iso (t b ) is derived by bootstrapping. For other parameters with simple mathematical expressions, the errors are calculated by error propagation.
After determining the bursts with an external plateau, we further split them into three groups.
The classification is made depending on whether the X-ray isotropic energy E X,iso and or the kinetic isotropic energy E K,iso are larger than 2 × 10 52 erg, corresponding to the maximum energy released by a magnetar. This is illustrated in Fig.6.
• Gold Sample: bursts for which E X,iso of the plateau phase is greater than 2 × 10 52 erg. The sample is made of 9 GRBs and the power law fits the light curves are shown in Fig.1. This corresponds to 9% of all the bursts.
• Silver Sample: bursts for which E X,iso of the plateau phase is less than 2 × 10 52 erg, but E K,iso is greater than 2 × 10 52 erg. It is also consistent with the theoretical prediction by a black hole central engine model. We note that because the computation of E K,iso relies on assumed values for the microphysical parameters, it is not possible to completely rule out the magnetar progenitor for the bursts. This sample is made of 69 bursts, and their light curve is presented in Fig.2. This corresponds to 68% of all the bursts.
• Bronze Sample: bursts for which neither E X,iso nor E K,iso are above 2 × 10 52 erg. For the 23 bursts in this sample, it is not possible to completely rule out the BH progenitor. Fig.3 shows the fitted light curves. This corresponds to 23% of all the bursts.
This classification indicates that 77% of the bursts have a black hole central engine, under the assumption that a magnetar emits its energy isotropically, see further discussion in §5.1.

γ-ray Energetics
The k-corrected isotropic energy released in the gamma-ray band (E γ,iso ) is retrieved from the published paper. For bursts which did not have published values, we derived them from the observed fluence (S γ ) in the detector band, and kcorrect to the rest-frame (1 − 10 4 ) keV using the spectral parameters.
We obtained the fluence and spectral parameters (including the α and β spectral indexes, peak energy E p ) of the Band function, from published papers (e.g., De Pasquale et al. 2016) or GCN Cir-culars Archive 13 . First, if a burst was detected by Swift/BAT and also by Fermi/GBM and/or Konus Wind, we consider the spectral parameters from Fermi /GBM or Konus Wind instruments. Second, if a burst was only detected by Swift/BAT instrument, and fit with a power-law or a cutoff power-law model, we make the following steps: (1). For the cutoff power law model, the fit parameters include α and E p (Sakamoto et al. 2011), so we adopt the typical value of β=-2.3. (2). For the single power law model, we derive E p by the empirical correlation, log E p = (2.76 ± 0.07) − (3.61 ± 0.26) log Γ BAT (Zhang et al. 2007a), Γ is the photon index of the Swift/BAT band, which was found from the Swift data archive, and adopt β with a typical value -2.3. (3). We finally test the results against the Amati relation (Amati et al. 2002), as shown in Fig.7. Our results are consistent with the Amati relation, which indicates that our method is reasonable.
We derive the k-correction (Bloom et al. 2001) factor in the rest-frame 1 − 10 4 keV band, and derive the isotropic energy E γ,iso as, the isotropic luminosity at the break (L X ) can be derived as: where F X,b is the flux at the break time (t b ), which was obtained from the light curve fit. The luminosity at break time (L X,iso ), and the isotropic release energy of the plateau phase (E X,iso ); together with the time intervals (t s , t e ) for the light curve fitting, the temporal indices of the plateau phase (α 1 ), and the normal decay or steeper decay phase (α 2 ), the break time (t b ), the X-ray fluence for the plateau phase (S X ), which is the integrated observed flux during the time interval (t s , t b ); the photon indices 14 for the shallow decay/plateau phase (Γ 1 ) and a normal or steeper decay phase (Γ 2 ) are obtained from the 13 https://gcn.gsfc.nasa.gov 14 The photon indices of the plateau phase were taken from the table published in Dainotti et al. (2010) Table also gives the spectral indices of the plateau phase and of the following segment, both obtained from either the Swift archive or from the closure relation.

Jet opening angle
In this section, we study the beaming correction with three methods and compare their results. The beaming correction is defined as It is used to correct the γ-ray and kinetic energies. We obtain the geometrically-corrected gamma-ray The later quantity is used for collimated outflows from a BH central engine (See §3. 2) The jet opening angle can be derived from the observational data (Rhoads 1999;Sari et al. 1999;Frail et al. 2001).
here θ j is the jet break time, E K,iso is the isotropic kinetic energy, and n is the medium density. If the light curve shows a jet break, we measure θ j using Equation (16). In our sample, 18 GRBs show a jet break. For the rest of the GRBs without a jet break, we use the observed time of the last data point as jet break time (t j ) to estimate θ j , and get an upper limit.
Another method to estimate the jet opening angle is by empirical relations (Nemmen 2012), θ j ≈ 5.0/Γ 0 , where Γ 0 is the initial Lorentz factor, also derived by the empirical relation as, here E γ,iso,52 is the isotropic energy, in units 10 52 erg. Fig.8 shows the distributions of beaming factors f b . We compare the distributions of these two methods. These distributions can be fitted with Gaussian functions, we have logf b = −1.68 ± 0.61 from using the late data point method, and logf b = −2.14 ± 0.47 from the empirical relation method. We find that the f b values of the first method is greater than the second method, and f b of θ j = 5 o lies between the two distributions. In Nemmen (2012), they adopted a coefficient 5 for the estimation. In turn, if it is initially defined by θ j ≡ a j /Γ 0 , we then can estimate a j . Fig.8b shows the distribution of a j , with the best Gaussian fitting a j = 5.11 ± 7.21. This result is consistent with the finding in Nemmen (2012). The results include θ j , f b , a j , θ j , f b are presented in Table 4.
We finally adopted the following criterion to calculate beaming factor: if a burst has an observed light curve jet-break, we then use it to estimate θ j ; otherwise, we adopt θ j = 5 o .

γ-ray Radiative efficiency
Another point we want to investigate is whether the γ-ray radiative efficiency is different between the samples. The radiative efficiency is defined as (Lloyd-Ronning & Zhang 2004), Since the kinetic energy (E K,iso ) increases with time during the plateau phase but stay the same during the normal decay phase, so the radiation efficiency η γ depends on which E K epoch is adopted.
Here we focus on the radiation efficiency at end of the shallow decay phase, η γ (t b ). The physical meaning for this is the efficiency of converting the spin-down energy of a BH/magnetar to γ− ray radiation. The X-ray efficiency can also be defined by a BH/magnetar to convert its spin-down energy to radiation, i.e.
In Fig.9, we compare E γ,iso and E K,iso (Fig.9a), L X,iso and L K,iso (Fig.9b) for different samples. The same distribution is also manifested in Figs. 10a and 10b. The typical value is log η γ (t b )= -0.53±0.56 for η γ and is log η x (t b )=-0.54±0.58 for η x . The radiative efficiencies in the γ and X-ray energy bands, in general, present the same distribution.

Central Engine Properties
In this section, we derive the properties of the central engine based on our observations, assuming two models.

Black Hole Central Engine
One of the leading models for the GRB central engine is a stellar mass BH surrounded by a hyper-accreting disk (e.g., Popham et al. 1999;Narayan et al. 2001;Frank et al. 2002;Kato et al. 2008;Abramowicz & Fragile 2013;Lei et al. 2013;Blaes 2014;Yuan & Narayan 2014), with typical accretion rate 0.01-1 M s −1 . There are two main energy reservoirs that provide the jet power namely the accretion energy in the disk that is carried by neutrinos and anti-neutrinos, which annihilate and power a bipolar outflow (e.g., Di Matteo et al. 2002;Kohri & Mineshige 2002;Gu et al. 2006;Chen & Beloborodov 2007;Janiuk et al. 2007;Liu et al. 2007;Lei et al. 2008Lei et al. , 2009 and the spin energy of the BH, which can be tapped through magnetic fields via the Blandford & Znajek (1977)(hereafter BZ) mechanism (see also, Lee et al. 2000;Li 2000;Lei et al. 2013).
As suggested by Lei et al. (2013), the jet may be dominated by the BZ power especially at late time We then connect the kinetic energy of the jet to the spin energy of BH through the relation where η is the efficiency of converting the spin energy to kinetic energy of jet with the BZ mechanism and f b is beaming factor (see §2.5 for a detailed discussion). For a maximally rotating BH, about 31% of the rotational energy can be taken out to power a GRB (Lee et al. 2000). This efficiency is not sensitive to the BH spin a • . Therefore we take η 0.3 in the following. From Eq.(20) and Eq.(22), one can constrain the BH spin parameter once we know the BH mass M • , jet opening angle θ j and the kinetic energy of jet E K,iso . Following Kochanek (2015), the typical value for BH mass in a long GRB is 6M . In this paper, we adopt M • = 3M for short GRBs 16 and 6M for long GRBs.
For a • > 0.1, the bimodule distribution of the index q can be explained by the spin evolution of the black hole. The black hole evolves with time during a GRB, since it would be spun-up by accretion while spun-down by the Blandford-Znajek mechanism (Lei et al. 2017). Due to and the different initial parameter set for different GRBs, the characteristics evolution vary from burst to burst. Therefore, we have q > 0 for the spin-down evolution and q < 0 for the spin-up case.
For a • < 0.1, the BH spin does not necessary need to be very low. Instead, the value we obtained based on Eq.(22) should be the change of BH spin ∆a • during a GRB. The very small value of ∆a • means that the BH spin does not evolve. This is possible when the BH spin comes to a critical value where the spin-up by accretion is comparable to the spin-down by the BZ mechanism, or the accretion rate is quite low. A clear evidence comes from the distribution of q when a • < 0.1, for which q is close to zero.
The correlation of a • − q for our different samples is presented in Fig.11. We find that the Gold sample cluster into the region which a • > 0.1. The results are consistent with the prediction of the black hole central engine model. Here we note that since the Bronze sample bursts could still have a BH central engine, we also calculate a • for them. Interestingly, we found all Bronze sample bursts have a • < 0.1.
The BH parameters (a • ), and duration of the prompt emission T 90 , the isotropic energy release in gamma-ray band E γ,iso , the isotropic kinetic energy E K,iso , electron spectral index p, luminosity inject index q and a redshift z (searched from the published papers or GCN Circulars Archive) are reported in Table 2 for the Gold and the Silver samples.
For a magnetized BH-accretion disk system, there is another magnetic mechanism: the magnetic coupling effect between the BH and the disk through closed magnetic field lines (van Putten 1999;Li & Paczyński 2000;Wang et al. 2002;Lei et al. 2009;Janiuk & Yuan 2010). Similar to the Blandford-Znajek mechanism, the magnetic coupling effect also extracts rotational energy from the spinning BH. Only if the BH spin is initially small, the magnetic coupling would act as an additional spin-up process. A similar discussion on this aspect was made by Dai & Liu (2012) within the context of the magnetar central engine model. In more general cases, the magnetic coupling effect would not significantly affect the BH spin evolution (Lei et al. 2009).

Magnetar Central Engine
A rapidly spinning magnetars model, which invokes a rapidly spinning, strongly magnetized neutron star is called a "magnetar". There are two critical magnetar parameters for a magnetar undergoing bipolar spindown: The initial spin period P 0 and the surface polar cap magnetic field B p , which can be derived using the characteristic luminosity L 0 and the spindown time scale τ . We derive the magnetar parameters (P 0 and B P ) using Eqs. (23)-(27), assuming isotropic emission from the magnetar.
For a magnetar, the characteristic spin down luminosity is: where R is the stellar radius 17 . In Eg.
(1), t is a measure of time in the rest frame, τ is the initial spin-down timescale which is defined by: whereṖ is the rate of period increase due to magnetic dipole radiation, I 45 is the stellar moment of inertia in units of 10 45 g cm −2 (Datta & Kapoor 1988;Weber & Glendenning 1992). The initial spin-down timescale, τ , in general, can be identified as the observed break time, t b , 17 We assume the typical value for R 6 =1, M = M 1.4 and I 45 ∼ 2. and the characteristic spin-down luminosity L 0 , in general, should include two terms: where L X,iso is the X-ray isotropic luminosity due to the internal dissipation of the magnetar wind, which is the observed X-ray luminosity of the internal plateau and measured at the observed break time (one can only derive an upper limit for external plateaus), and is the kinetic isotropic luminosity which is injected into the blast wave during the energy injection phase. It depends on the isotropic kinetic energy E K,iso after the injection phase is over and the observed break time t b . The isotropic kinetic energy E K,iso can be derived from afterglow modeling, as discussed above.
Since it is not possible to completely rule out the magnetar progenitor for the Silver sample, we investigate the magnetar parameters for both the Silver and Bronze samples.
The derived values of P 0 and B P are presented in Table 3 and Fig.12. We found that all bursts in the Bronze sample fall into the reasonable range in the (P 0 − B P ) plot (Fig. 12). B P with a typical value ∼ 10 15 G and ranges within [10 14 G,10 16 G], no case have B P below 10 14 G and only one case is greater than 10 16 G. We find most derived P 0 are close to 1 ms, this is consistent with the break-up limit of a neutron star, which is about 0.96 ms (Lattimer & Prakash 2004) 18 .
However, we also found that approximately half of the Silver sample bursts have B P < 10 14 G. Moreover, most of the bursts have P 0 is far away below to 1 ms. The result implies that the Silver sample bursts are more likely to have a BH central engine.
Using the derived value of P 0 , the total rotation energy of the millisecond magnetar within this scenario is, where I is the moment of inertia, Ω = 2Π/P 0 is the initial angular frequency of the neutron star, M 1.4 = M/1.4M , and E rot is the maximal total energy available and can be compared to the observed and inferred energies (see §3.3).
In conclusion, the properties of the bursts in the Bronze sample are consistent with the predictions of magnetar central engine model.

Total Energy and Correlations
In Fig.13, we compare the inferred E γ,iso +E K,iso with the total rotation energy E rot of the BH (Silver; §3.1) and the millisecond magnetar (Bronze; §3.2). It is found that the GRBs are generally above and close above the E rot = E γ,iso + E K,iso line. This is consistent with the hypothesis, namely, all the emission energy ultimately comes from the spin energy of the BH or magnetar.
In Fig.14, we show the T 90 (rest frame) as a function of various energies for both the BH and magnetar candidate samples. All isotropic energies for most cases are above 2×10 52 erg while all energetic for most cases with a beam correction are below 2×10 52 erg.
In Fig.15, we compare the histograms of the isotropic energies (E γ,iso , E K,iso , E X,iso ) of the BH and magnetar candidate samples. The typical values of all kinds of energy (E γ,iso , E K,iso , E X,iso ) for the BH candidate samples (Gold+Silver) are systemically larger than the magnetar sample (Bronze). Fig.16 displays the T 90 (rest frame) as a function of redshift for the samples. Two interesting findings: first, 4 out of 5 short GRBs are in the magnetar candidate sample (Bronze); second, the redshifts of all 5 short bursts below 1. Fig.17 shows the correlations of L X,iso − E γ,iso and L X,iso − t b /(1+z) for the entire sample (BH candidate Gold+Silver samples and magnetar candidate Bronze sample). We find that a higher isotropic γ-ray energy implies a higher isotropic X-ray break luminosity and the bursts of the Gold sample cluster into the region of high isotropic γray energy and isotropic X-ray break luminosity. The best fit to the correlation is log L X,iso,49 = (−1.33 ± 0.10) + (1.06 ± 0.09) log E γ,iso,52 , with a Spearman correlation coefficient 0.78 and a chance probability p < 10 −4 . For the L X,iso -t b /(1+z) correlations, the best fit is log L X,iso,49 = (2.74 ± 0.40) + (−1.09 ± 0.09) log t b /(1+z), with a Spearman correlation coefficient 0.79 and a chance prob-ability p < 10 −4 . The anti-correlation between the luminosity and the rest-frame duration of the plateau phase are consistent with the previous findings (Dainotti et al. 2008(Dainotti et al. , 2010(Dainotti et al. , 2013(Dainotti et al. , 2015(Dainotti et al. , 2017Dall'Osso et al. 2011;Rowlinson et al. 2013Rowlinson et al. , 2014Rowlinson et al. , 2017. In Fig.18, we show the scatter plots of η X compare η γ (t b ), E γ,iso , E K,iso (t b ) and, E rot (Bronze). One interesting finding is the η X -E K,iso (t b ) correlation present an anti-correlation for the BH candidate samples (Gold+Silver) while present an positive correlation for the magnetar candidate sample (Bronze).

Comparison of Statistical Character
Since most of the bursts in our sample have an external plateau, this indicates that our GRB samples comply with the external shock models well, and they offer an excellent sample to study external shock model parameters. Here we investigate whether there are noticeable differences between the BH candidate (Gold+Silver) and magnetar candidate (Bronze) bursts for the observation properties and model parameters.
In addition, we also investigate the external plateau and the internal plateau in our bursts. Fig.19 displays the distributions of temporal indices α in different sub-samples and different temporal segments. They are all well fitted with Gaussian distributions for each sample. For α x,1 distributions, we have α x,1 =0.55±0.30 for the BH candidate samples, and α x,1 =0.32±0.48 for the magnetar candidate sample. The same analysis based on the type of plateau gives α x,1 =0.60±0.47 for internal plateau sample, α x,1 =0.50±0.31 for external plateau ones (Fig.19a). For α x,2 distributions, we have logα x,2 =0.18±0.11 for the BH candidate samples, and logα x,2 =0.07±0.13 for the magnetar candidate. For the internal/external plateau samples, one has logα x,2 =0.46±0.33 for the internal plateau sample, and logα x,2 =0.17±0.11 for the external plateau sample, respectively (Fig.19b).
The temporal indices (both α 1 and α 2 ) for the BH candidates samples (Gold+Silver), in general, steeper than the magnetar sample (Bronze); For the internal and external example, the pre-break temporal indices α 1 are consistent with each other, but the post-break temporal indices α 2 for internal sample is much steeper than the external sample, consistent with an internal origin. The degrees of the break ∆α for the BH candidates samples (Gold+Silver) is significantly less than the Silver and the magnetar sample (Bronze). For Γ x,1 distributions, we have Γ x,1 = 1.91 ± 0.20 for the BH candidate samples, and Γ x,1 = 1.99 ± 0.23 for the magnetar candidate sample (Fig.20a). For Γ x,2 distributions, one has Γ x,2 = 1.94 ± 0.20 for the BH candidate samples, and Γ x,2 = 1.88 ± 0.25 for the magnetar candidate sample (Fig.20b). For ∆Γ x distributions, one has ∆Γ x = 0.03 ± 0.19 for the BH candidate samples, and log∆Γ x = −0.08 ± 0.24 for the magnetar candidate sample (Fig.20c).
The spectral indices (both Γ 1 and Γ 2 ), and the degrees of the break ∆Γ for different sub-samples (the BH candidate samples and the magnetar candidate sample), in general, present the same distributions. For the internal and external samples, a change is not found in the external sample while a significant change in the internal sample is found.
The results indicate that there is no spectral evolution for the sub-samples and the external plateau sample, while internal plateau sample has a significant spectral change between a plateau phase and the followed steeper decay phase. This is consistent with the internal plateau having an internal physical origin. Fig.21 shows the distributions of the observed break times (t b ).

Break time t b
A Gaussian fit to the distributions gives log(t b /s)=(4.05±0.68), log(t b /s)=(4.38±0.63) for the BH candidate samples and the magnetar candidate sample respectively. Separating the "internal" plateau and "external" plateau samples, we have log(t b /s)=(4.22±0.54) for the "internal" plateau sample, log(t b /s)=(4.07±0.72) for the "external" plateau sample. So both external and internal are consistent for t b .
The typical break time for the BH candidates samples (Gold+Silver) are statistically earlier than that of the magnetar candidates sample (Bronze). The similar results can be found between the internal (earlier) and external (later) samples. The result indicate that the end energy injection for the "external plateau" is in general on average earlier than "internal plateau".

Electron spectral index p
We derive the electron spectral indices p, derived using the temporal indices, The p distributions for our sub-samples: p = 2.63 ± 0.63 for the BH candidate samples, and p = 2.49 ± 0.15 for the magnetar candidate sample (Fig.22). The global p value ranges within [1.89, 3.70], and it has a Gaussian distribution with a center value p = 2.69 ± 0.05, except 11 GRBs are internal plateau", all the rest have p > 2.0, except one case (GRB 100302A) had p < 2.0 (Fig.22).
The results are consistent with the typical value of p for relativistic shocks due to 1st-order Fermi acceleration (e.g., Achterberg et al. 2001;Ellison & Double 2002).

Energy injection parameters q
Similarly, the luminosity injection index q (see Eqs. 1) can be determined from the temporal index α and the spectral index β for different afterglow models, here α and β are measured in the shallow decay segment preceeding the break.
The q distributions for our sub-samples are: q = 0.24 ± 0.50 for the BH candidate samples, and q = 0.38 ± 0.35 for the magnetar candidate sample (Fig.23). Also, for the internal plateau, we have q = 0.39 ± 0.35, and external plateau, we have q = 0.22 ± 0.35.
The different sub-samples present different central values, the typical value for the BH candidate samples (Gold+Silver) is slightly less than the magnetar candidate sample (Bronze). The typical q value for the internal sample is greater than the external sample. Totally, the central q value for sub-samples is around 0.3, which are consistent with the predictions of an energy injection model (Zhang et al. 2006).

Discussions and Conclusions
In the analysis above, several assumptions have been made which, if relaxed, makes the magnetar energy budget even more strained. For instance, the radiative efficiency is unknown but could indeed be low (e.g. Beniamini et al. 2017). Therefore the limiting, radiative energy for a magnetar burst could be only a fraction of the maximum allowed rotational energy. Moreover, the energy estimated from the X-ray afterglows could underestimate the true energy, if they are strongly suppressed by Compton cooling (Beniamini et al. 2015). On the other hand, the total energy from a magnetar GRB should come from the rotational energy of the magnetar itself. Except for the kinetic energy E K,iso (t b ), which measured at the end of the plateau phase that is considered in this paper, the energy of the GRB itself (E γ,iso ) including GeV emission should also be considered and should be added into the total energy, and therefore could cause some GRBs in our magnetar sample to exceed the maximum energy allowed by the magnetar. The main assumption made is, however, on the sphericity of the magnetar wind which we discuss in the next section.
Furthermore, the related detailed studies have been done by many authors (e.g., Dall'Osso et al. 2011;Rowlinson et al. 2013Rowlinson et al. , 2014Rowlinson et al. , 2017Rea et al. 2015). For instance, Dall'Osso et al. (2011) considered long GRBs are formed by rapidly spinning, newly born magnetars, and claimed that the high spindown luminosity cause by MRD losses represents a natural mechanism for prolonged energy injection in the external shock in the first hours after the formation of the magnetars. The individual light curves of the cases can be fitted well by their model, and the initial spin period P 0 and magnetic dipole field B p which are derived from their fitting also match well values expected in the magnetar model. The results are consistent with our finding. The simulation results in Rea et al. (2015) suggest that the if a rapidly spinning magnetar powered the GRB X-ray plateau, magnetars are required to have two different progenitors and formation path and different magnetic field formation efficiencies. Rowlinson et al. (2013) show evidence of energy injection for many short GRBs. In addition, they suggest that the remnant of NS-NS star mergers could produce an unstable magnetar rather than collapse directly and immediately to a black hole. The unstable magnetar power an Xray plateau. The results show that nearly half of events forming magnetars could collapse to a BH after a few hundred seconds. Moreover, in Rowlinson et al. (2014), they studied GRB magnetar central engines using the observed plateau luminosity and duration correlation based both long and short GRBs, and they found the magnetar emission most likely is narrowly beamed (< 20 o ) and has a lower efficiency ( 20%) of conversion of rotation energy to observed X-ray plateau emission.

Sphericity of the magnetar wind
One important assumption in the analysis above is whether the emitted energy is isotropic or not, since this affects the derived energetics of the bursts. In the framework of a black hole central engine, the rotational energy can naturally be assumed to be extracted through the jet and is deposited into the X-ray plateau. In this case the a beaming correction needs to be made: On the other hand, for a magnetar central engine, the rotational energy is initially emitted in an isotropic wind. Indeed, for short GRBs, the total energy is typically expected to be emitted quasi-spherical (e.g. Metzger & Piro 2014;Fan & Xu 2006). Here, we have also assumed the energy to be emitted quasi-isotropic even for long GRBs during their X-ray afterglow phase (see, e.g. Mao et al. 2010;Dall'Osso et al. 2011;Mazzali et al. 2014). Then only a small part of the rotational energy of a magnetar is injected into the external shock of the jet. This means that for the magnetar sample, using this assumption, we should not use the beaming correction, so that We note, however, that it is still unclear whether the energy release for long GRBs during their X-ray afterglow phase is spherical or collimated (Metzger et al. 2011) and there needs to be further theoretical and numerical modeling, as well as observational constraints, to understand the late time energy release from a magnetar. Indeed, there are arguments for a collimated outflow from a magnetar central engine. For instance, the progenitor star of long GRBs is expected to cause a collimation of the emitted energy, at least during the prompt phase (Bucciantini et al. 2007(Bucciantini et al. , 2009, and this would change the derived total energetics. However, for very late times, such as during the X-ray plateau phase, which is hundred seconds to a day following the trigger of the GRB, there should be less material surrounding the magnetar. Therefore, the magnetar emission could be expected to come out more isotropically. Similarly, Lü & Zhang (2014) found, by using Swift data of GRB afterglows, that while an isotropic assumption for the emitted energy is reasonable for short GRBs, for most long GRBs this assumption leads to unrealistic values of the derived rotation period of the neutron star. Lü & Zhang (2014) therefore argue that this means that the assumption of isotropic energy release during the afterglow phase is not correct. We note, however, that it could also point towards that the progenitor is not a magnetar.
Oppositely, an indirect argument for quasiisotropic magnetar winds is given in Mazzali et al. (2014) who studied the energetics of GRBs associated with SNe. They noted that the SN energy is narrowly distributed around the characteristic magnetar energy of ∼ 10 52 erg and is always larger than the GRB energy. On the other hand, the GRB energy varies a lot. They therefore suggest that the central engine for the SN associated GRBs is a magnetar and that the SN is driven by a close-to-spherical magnetar wind. A consequence is that the jet does not drive the stellar explosion, indeed, it is not clear how the jet would transfer energy to the star after it has escaped. Mazzali et al. (2014) had to make an assumption, though, that the magnetar wind is not fully isotropic based on the fact that SNe, in general, is observed to have some degree of asphericity (e.g. Maeda et al. 2002;Tanaka et al. 2007). The asphericity is, however, assumed to be much smaller than the highly collimated, jetted flows from BH central engines, with only f b <∼ 0.5−0.2. Moreover, Greiner et al. (2015) argued for a magnetar central engine in GRB 111209A, since the associated SN 2011kl did not show any evidence for Ni production. However, its properties are peculiar compared to other GRB/SN cases and might thus not be representative (Cano et al. 2017). In the case of GRB 111209A a correction for collimation of f b < 1/50 was required for E γ , which would suggest an initial jetted outflow. On the other hand, considering the gamma-ray radiative efficiency found in §2.6 (Eq.18, Fig.9 and Fig.10) and the derived E X,iso the total energy budget could be strongly strained for the simple estimates of a magnetar central engine, in particular if E X is assumed to be isotropic. Indeed, Cano et al. (2016) argues that additional energy sources, such as radioactive heating, have to be considered apart from the magnetar itself for a consistent scenario.
Following these lines of argumentation, the assumption made in this paper of f b = 1 might be too strong. Instead, by employing a factor of f b < 0.5 − 0.2, as argued by Mazzali et al. (2014), the energies become smaller and more bursts are consistent with the magnetar bound. For instance, using f b = 0.5 (0.2), there are 2 (18) bursts in the Silver sample that would be classified as possibly harboring a magnetar central engine and thus belonging to the Bronze sample (see Table 5). We also calculate the beaming correction, f b , that individual bursts need in order to reach the magnetar limit (green line in Fig.8a). These f b values are significantly larger than the f b values obtained from the jet break. On the other hand, it shows that only a slight anisotropy is needed in order for the total energy to get below the magnetar limit. In such cases, the central engine could still be a magnetar.
This points towards the fact that there could be a continuum of degrees of asphericity (or collimation) of the magnetar wind. By assuming a fully isotropic wind we argue in this paper, that most bursts are inconsistent with the magnetar central engine (our Gold and Silver samples). Only a small fraction of bursts, both short GRBs and long GRBs in our Bronze sample, are consistent with a magnetar central engine assumption. Oppositely, Lü & Zhang (2014) assume a high degree of collimation and find a large fraction of all bursts to be consistent with the magnetar central engine. As argued above the degree of asphericity could lie in between the assumption of Lü & Zhang (2014) (high degree of collimations) and this paper (isotropic wind), and could, for instance, be determined by varying mass of the envelope surrounding the magnetar which would have a different impact on the focusing of the wind.

Conclusions
In this paper, we systematically analyzed the Swift/XRT light curves having a plateau phase in 101 GRBs which were detected before May 2017. Assuming an isotropic magnetar wind we draw the following conclusions: • A Gold sample includes 9 out of 101 GRBs.
Both E X,iso and E K,iso exceed the energy budget of a magnetar. Therefore these most likely harbors a BH central engine.
• A Silver sample includes 69 out of 101 GRBs. Estimations of their kinetic energy E K,iso is larger than the energy budget of a magnetar. Likewise, these most likely have a BH central engine.
• A Bronze sample includes 23 out of 101 GRBs. These bursts have energy output lower than the maximal energy budget of magnetars, and as such can be considered candidates for the magnetar central engine.
• We further test the data with the BZ mechanism based on the black hole model, and find that the observations of the 'Gold' and 'Silver' samples are consistent with expectations. We also test the magnetar model for the Silver and Bronze samples, and find that the magnetar surface magnetic filed (B p ) and initial spin period (P 0 ) fall into a reasonable range for the Bronze sample but not for the Silver sample. This implies that the Bronze sample is consistent with the magnetar central engine, but the Silver sample, in general, is inconsistent with the magnetar model.
We also compare the properties between the two central engine candidate samples, and two type of plateau bursts.
• BH candidate (Gold+Silver) & Magnetar candidate (Bronze) We found the observational properties (α 1 , α 2 , ∆α, Γ 1 , Γ 2 , and ∆Γ) and the model parameters (q), in general, do not differ significantly. This indicates it could be insensitive to the central engine or, in the case these parameters were to be sensitive to the central engine, it would argue for a common central engine type. However, one interesting result is that the p value distribution for the Magnetar candidate sample is significantly narrower than the BH candidate sample. A significant difference is also found for t b , z, T 90 , L X,iso , energies (E X,iso , E K,iso , and E γ,iso ), and p, which could imply a different physical reason between the different samples.

• Internal plateau and External plateau
For the internal and the external plateaus, we found α 2 , ∆α, Γ 2 , ∆Γ, t b , and q in general, present different distributions. The results are consistent with the expectation for different physical progenitors for two types of plateaus; while the α 1 and Γ 1 generally show similar distributions.
We conclude that, under the assumption of isotropic energy release from a magnetar, most GRBs should harbor a black hole central engine.
We acknowledge the use of the public data from the Swift data archive and the UK Swift Science Data Center. We thank Damien Bégué, Yu Wang, David Yu, Hüsne Dereli, Jin-Jun Geng, Yun-Feng Liang and Christoffer Lundman for helpful discussions.
We appreciate the valuable comments by the referee and thanks Bing Zhang for useful discussions that greatly improved the paper. This work is supported by    b The X-ray photo indices before and after break time.
c The fluence during the plateau phase, calculated by integrating the fitting light curve from start to break time.
d The luminosity at break time, in units of 10 47 ergs −1 . The X-ray release energy during plateau phases, which is calculated by integrating the luminosity during the time interval of the plateau phase. In units of 10 50 erg.  a E γ,iso is searched by published papers or is calculated using fluence and redshift extrapolated into 1-10 4 keV (rest frame) with a spectral model and a k -correction, in units of 10 52 erg.
b The kinetic energy, which is derived by normal decay phase. The initial kinetic energy (upper limit), which is estimated by kinetic energy at break time (t b ).
c The BH spin parameter a• , in units 10 −1 second.
d We did not calculate the kinetic energy for the bursts which have an internal plateau.      b a j value was estimated by set the ratio of the θ j,0 equal to θ j .
a How many the Silver bursts become to the Bronze bursts after with a beaming correction (E K,iso )?
b How many GRB-SN bursts from the Silver bursts become to the Bronze bursts after with a beaming correction (E K,iso )?
37 10 -4 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3 10 4 E γ, iso 10 52 (erg) The E γ,iso − E K,iso scattered plot for all the GRBs (excluded internal plateaus). Slanted dashed lines mark the constant γ − ray efficiency (η γ ) lines. E K,iso is calculated at t b ; (b) The L X,iso − L K,iso scattered plot for the BH/magnetar samples (excluded internal plateaus). The constant X-ray efficiency η X lines are over plotted.  -The inferred magnetar parameters, initial spin period P 0 vs. surface polar cap magnetic field strength B p , derived from the Silver and Bronze samples. The vertical solid line is the breakup spin-period for a neutron star (Lattimer & Prakash 2004   . q values are derived from in a shallow decay segment. The same symbol as Fig.19 but for the luminosity injection index q distributions.