{'fraction_non_alphanumeric': 0.041013141161509115, 'fraction_numerical': 0.03743641373463332, 'mean_word_length': 4.054908263023973, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The Second Plateau in X-ray Afterglow Providing Additional Evidence for Rapidly Spinning Magnetars as the GRB Central Engine

Evidence for the central engine of gamma-ray bursts (GRBs) has been collected in the Neil Gehrels \emph{Swift} data. For instance, some GRBs show an internal X-ray plateau followed by very steep decay, which is difficult to be interpreted within the framework of a black hole (BH) central engine, but are consistent within a rapidly spinning magnetar engine picture. The very steep decay at the end of the plateau suggests a sudden cessation of the central engine, which is explained as the collapse of a supra-massive magnetar into a black hole when it spins down. Here we propose that some additional evidence, such as a second X-ray plateau feature would show up, if the fall-back accretion could activate the newborn BH and sufficient energy could be transferred from the newborn BH to the GRB blast wave. With a systematic data analysis for all long GRBs, we find three candidates in \emph{Swift} sample, i.e., GRBs 070802, 090111, and 120213A, whose X-ray afterglow light curves contain two plateaus, with the first one being an internal plateau. We find that in a fairly loose and reasonable parameter space, the second X-ray plateau data for all 3 GRBs could be well interpreted with our proposed model. Future observations are likely to discover more similar events, which could offer more information of the properties of the magnetar as well as the newborn BH.


INTRODUCTION
Gamma-ray bursts (GRBs) have been extensively explored since their discovery more than 50 years ago, but the nature of the GRBs central engine remains a mystery. In the literature, two main kinds of central engine have been well discussed: hyper-accreting black holes (BH) or rapidly spinning magnetars (Zhang 2018, for a review). It has long been proposed that some interesting signatures showing in some GRB's X-ray afterglow could help us to determine their central engines (Dai & Lu 1998;Rees, & Mészáros 1998;Zhang & Mészáros 2001;Zhang et al. 2006;Nousek et al. 2006). For instance, systematic analysis for the Swift GRB X-ray afterglow shows that bursts with X-ray plateau features likely have rapidly spinning magnetars as their central engines (Liang et al. 2007;Zhao et al. 2019;Tang et al. 2019). In particular, when the X-ray plateau is followed by a steep decay with temporal decay index α 3, hereafter called "internal X-ray plateaus", the sharp decay at the end of the plateau is difficult to be interpreted within the framework of a BH central engine, but are consistent within a magnetar engine picture, where the abrupt decay is understood as the collapse of a supra-massive magnetar into a BH after the magnetar spins down (Troja et al. 2007;Lyons et al. 2010;Rowlinson et al. 2010Rowlinson et al. , 2013Lü & Zhang 2014;Lü et al. 2015;Gao et al. 2016a;De Pasquale et al. 2016;Zhang et al. 2016).
Recently, Chen et al. (2017) proposed that if the sudden drop after the internal plateau indeed indicates the collapse of a supra-massive NS into a BH, signatures from this newborn BH should be expected. For instance, for long GRBs if a fraction of the envelope material fall back and activate the accretion onto the newborn BH (Kumar et al. 2008a,b;Wu et al. 2013;Gao et al. 2016b), the hyper-accreting BH system can launch a relativistic jet via Blandford-Znajek (BZ) mechanism (Blandford & Znajek 1977). Supposing a fraction ξ of the jet energy would undergo internal dissipation, some detectable signals, such as a X-ray bump following the internal plateau is expected (Chen et al. 2017). Searching through Swift -XRT data archive, the authors found a particular case, GRB 070110, which showed a small X-ray bump following its internal plateau, and successfully interpreted its multi-band data with their model. But it is worth noting that for GRB 070110, Troja et al. (2007) found the optical data indicate a decaying light curve feature sitting underneath the X-ray plateau without showing the rapid drop for the "internal plateau", which means unlike the X-ray band, the optical emission trend could be hardly disturbed by the late central engine activity.
Note that most internal plateaus are not found to be followed by the X-ray bump, it may be because that the late time fallback processes are intrinsically weak or the energy fraction for the internal dissipation (ξ) is relatively small. For the latter situation, a food fraction 1 − ξ of the BZ jet energy would continuously inject into the GRB blast wave, which may generate a second plateau in X-ray afterglow, if the injected energy is comparable or even larger than the blast wave kinetic energy. We thus propose that GRBs with two X-ray plateaus (the first one is internal plateau) may provide further support to the magnetar central engine model.
In this work, we first study how the energy injection from the fallback accretion onto the newborn BH would alter the GRB afterglow emission and analyze the influence of model parameters on the theoretical lightcurves. In section 3, we systematically search for GRBs with two X-ray plateaus from the Swift-XRT sample. Considering that long GRBs, which are likely related to the core-collapse of massive stars (Woosley 1993;MacFadyen & Woosley 1999), are more easier to have enough envelop material to provide late time fall-back accretion than short GRBs, that have been proposed to originate from the merger of two neutron stars (Eichler et al. 1989;Narayan et al. 1992) or the merger of a neutron star and a black hole (Paczynski 1991), here we focus on long GRB samples 1 . We find three candidates, GRB 070802, GRB 090111 and GRB 120213A, whose second plateaus could be well fitted by the model. The conclusion and implications of our results are discussed in Section 4. Throughout the paper, the convention Q = 10 n Q n is adopted in c.g.s. units.

MODEL DESCRIPTION
Rapidly spinning magnetars have long been proposed as one of the candidate of GRB central engines. In this scenario, a collimated jet could be launched by invoking (1) hyper-accretion onto the NS (Zhang & Dai 2009Bernardini et al. 2013); (2) magnetic bubbles from a differentially millisecond proto-NS (Dai et al. 2006); (3) or from a protomagnetar wind (Metzger et al. 2011). The internal dissipation of the jet could power the prompt gamma ray emission of a GRB and the interaction between the jet and the ambient medium could produce a strong external shock that gives rise to bright broadband afterglow emission (Gao et al. 2013, for a review). After launching the jet, the magnetar would also eject a near-isotropic Poynting-flux-dominated outflow, the internal dissipation of which could power a bright X-ray emission, whose temporal profile would follow the spin-down profile of the magnetar, i.e., t 0 at early stage and decay as t −2 after the magnetar spinning down (Zhang & Mészáros 2001). If this emission component is brighter than the external shock afterglow emission, an X-ray plateau would show up. Sometimes the central engine magnetar might be a supra-massive NS, which would collapse to a BH when a good fraction of its rotational energy is lost and the centrifugal support can no longer support gravity. In this case, the X-ray plateau emission would suddenly stop and follows by a sharp decay at the end of the plateau, due to the abrupt cessation of the magnetars central engine. This explains the internal X-ray plateau discovered in both long and short GRBs (Troja et al. 2007;Lyons et al. 2010;Rowlinson et al. 2010Rowlinson et al. , 2013Lü & Zhang 2014;Lü et al. 2015). After the collapse of the magnetar, if a fraction of the envelope material (especially for long GRBs) fall back and activate the accretion onto the newborn BH, the rotational energy of the BH could be extracted via BZ mechanism, and a good fraction of the energy would eventually inject into the afterglow blast wave. If the injected energy is comparable or even larger than the blast wave kinetic energy, the broadband afterglow lightcurve could be significantly altered, for instance, a second plateau in the X-ray afterglow may emerge. Figure 1 presents a physical picture for several emission components at different temporal stages. Here we focus on study how the energy injection from the fallback accretion onto the newborn BH would alter the GRB afterglow emission and analyze the influence of model parameters on the theoretical lightcurves. Model illustration: the "internal plateau" is powered by the spin-down power from a supra-massive magnetar and the steep decay marks the collapse of the magnetar into a BH when it spins down. In addition, the second plateau is caused by energy transfer from new born BH through BZ mechanism to GRB blast wave.

The Fallback Accretion onto the Newborn BH
Assuming the fallback accretion could trigger the energy extraction from the newborn BH via BZ mechanism, in this case, the BZ power from a BH with mass M • and angular momentum J • could be estimated as (Lee et al. 2000;Li 2000;Wang et al. 2002;Lei et al. 2005Lei et al. , 2013Lei et al. , 2017McKinney 2005;Lei & Zhang 2011;Chen et al. 2017;Liu et al. 2017;Lloyd-Ronning et al. 2018) where m • = M • /M is the dimensionless BH mass and a = J • c/(GM 2 • ) is the dimensionless spin parameter of the BH. Here, F (a) = [(1 + q 2 )/q 2 ][(q + 1/q) arctan q − 1], q = a/(1 + √ 1 − a 2 ), and B •,15 is the magnetic-field strength threading the BH horizon in units of 10 15 G. The evolution of the BZ jet power depends on m • , a, and B • .
The evolution of the BH spin and mass governed the competition between spin up by accretion and spin down by the BZ mechanism. The evolution equations of the BH mass M • and the BH spin a are given by Wang et al. (2002) and whereṀ is the BH accretion rate, and T BZ is BZ magnetic torque (Li 2000;Lei & Zhang 2011;Lei et al. 2017) where E ms and L ms are the specific energy and the specific angular momentum at innermost radius r ms of the disk (Novikov & Thorne 1973): Here, R ms = r ms /r g and r g = GM • /c 2 . The innermost stable radius of disk is (Bardeen et al. 1972) where Z l ≡ 1 + (1 − a 2 ) 1/3 (1 + a) 1/3 + (1 − a) 1/3 and Z 2 ≡ (3a 2 + Z 2 1 ) 1/2 for 0 a 1. The strength of the magnetic field is the major uncertainty in estimating the BZ jet power. By assuming the ram pressure of the innermost part of the disk balances the magnetic pressure on the BH horizon (Moderski et al. (1997)), one can estimate the magnetic field strength B • , where r • is the radius of the BH horizon. With this assumption, the BZ jet power could be written as a function mass accretion rate and BH spin, i.e. (Wu et al. 2013;Chen et al. 2017) whereṁ =Ṁ /(M s −1 ) is the dimensionless BH accretion rate. The accretion rate of BH can be estimated by adopting a simple model described in Kumar et al. (2008a) where the viscous timescale τ vis ∼ 1/αΩ K , here α is the standard dimensionless viscosity parameter, Ω K is the Kepler angular velocity of accretion disk. The mass of the disk M d evolves with time, it increases as a result of fall-back from the envelope and decrease as a result of accretion. Thus, one has (Kumar et al. 2008a;Lei et al. 2017) Combining Eqs (10) and (11), one can obtain the accretion rate onto the BH (Kumar et al. 2008a;Lei et al. 2017), The evolution of the fall-back accretion rate is described with a broken-power-law as (Chevalier 1989;MacFadyen et al. 2001;Zhang et al. 2008;Dai & Liu 2012) where t 0 is the start time of the fall-back accretion in the local frame, t p is peak time of fallback and andṀ p is the peak fall-back rate.
For the rapid accretion case, τ vis t, the BH accretion rate would follow the fall-back rate, i.e.,Ṁ =Ṁ fb . For a large value of the viscosity timescale τ vis , the BH accretion rate would be flat until t > τ vis , and then starts to decline with time, see Figure 7 in Lei et al. (2017).

Energy injection into the GRB afterglow blast wave
The energy flow from the BZ process would continuously inject into the external shock, and cause a significant raise to the Lorentz factor of the blast wave Γ, which may produce the second plateau following the steep decay. Huang et al. (2000) proposed a generic dynamical model to describe the dynamical evolution of GRB outflow, which has been widely applied for modeling the afterglow light curve. Based on their model and taking the energy injection into account, the evolution equation of the outflow's bulk Lorentz factor can be written as (Liu & Chen 2014)  where M sw is the swept-up mass by shock, M ej is the initial mass of the GRB outflow. The blast wave energy continuously increase with time, due to the continuously injected energy from the BZ process into the blast wave. The initial kinetic energy of GRB outflow can be estimated as E 0 = Γ 0 M ej c 2 , and the injected energy from the BZ process can be calculated as 2 E inj = (1 − ξ)L BZ dt. The total energy in the blast wave could be expressed as We introduce a parameter η ≡ E inj /E 0 to denote the ratio between the injected energy and the initial energy in blast wave. For η > 1, the injected energy is dominated in the total energy of the blast wave, otherwise the initial kinetic energy of GRB outflow is dominated.
In order to obtain the shock dynamical evolution, three additional differential equations are required (Huang et al. 2000). The evolution of the radius of shock R, the swept-up mass M sw , and the opening angle of the jet θ are described by Huang et al. (2000) where n is the number density of the unshocked ISM and c s is the sound speed, where the adiabatic indexγ = (4Γ + 1)/(3Γ).
For an electron with an energy γ e m e c 2 in the co-moving frame of the shock, the observed frequency from sychrotron emission is (Rybicki, & Lightman 1979) where m e is the electron mass, q e is the electron charge, the bulk Lorentz factor Γ is introduced by transferring the shock co-moving to the observer rest frame and B is the magnetic field strength in the shock co-moving frame. In GRB problems, one usually assumes the magnetic field density (B 2 /8π) is a fraction of the internal energy of post-shocked medium with a shock equipartition parameter B . Therefore, the co-moving magnetic field is in form of (Sari et al. 1998) The energy distribution of shock accelerated electron is usually assumed to be a power law, with N (γ e )dγ e ∝ γ −p e dγ e . The minimum Lorentz factor of the electron, γ m could be obtained by the conservation laws of particle number and energy. Assuming a constant fraction e of the post-shock internal energy goes into the electrons, one has (Sari et al. 1998) Another critical electron Lorentz factor is the cooling Lorentz factor γ c . When γ e > γ c , electrons would loss most of their energies by synchrotron radiation, otherwise the cooling caused by the radiation can be ignored. If synchrotron radiation is dominated, the cooling Lorentz factor is given by (Sari et al. 1998) where t refers to time in the rest frame of observer, σ T is Thomson cross-section. As shown in Eq (20), electrons with different Lorentz factors γ e have different radiation frequencies ν(γ e ). Two characteristic frequencies, ν m = ν(γ m ) and ν c = ν(γ c ) would determine the synchrotron spectrum. The evolution of these frequencies with time can be derived from shock dynamics.
For a given dynamical time t, if γ c > γ m , it means only a small fraction of electrons could be cooled, This is called slow cooling regime. In this case, for an observational frequency ν, the synchrotron spectrum is described a broken power law characterized by ν m and ν c as follows (Sari et al. 1998) On the other hand, for γ c < γ m , all the accelerated electrons could be cooled in the dynamical timescale t. This is named fast cooling regime. The radiation spectrum of the shock is (Sari et al. 1998) where F ν,max is the observed peak flux at luminosity distance D L from the source, which can be estimated as (Sari et al. (1998)) where N e,tot = 4πnR 3 /3 is total number of swept-up electrons in the post-shock medium.
With the dynamics and radiation equations described above, we can calculate the time-dependent spectra of the blast wave emission and the light curves for any observational frequencies. For the purpose of this work, here we focus on the X-ray light curve within the Swift-XRT energy band (0.3 − 10 keV). Note that the dominated radiation mechanism for X-ray afterglow emission is synchrotron radiation (Mészáros & Rees 1997;Sari et al. 1998), therefore, inverse Compton mechanism is ignored in our calculation.

Influence of the Model Parameters on the X-ray Light Curve
In this subsection, we show the numerical results of our model. Here, we explore the influence of the model parameters on the theoretical X-ray light curves. There are several free parameters in our model, which can be divided into two categories. The first category is related to the BZ process and the newborn BH, including the initial mass of the newborn BH M •,0 , the initial BH spin a 0 , the viscosity timescale of disk τ vis , peak time of fallback t p and the peak fallback rateṀ p . The second category is associated with the external shock, including initial kinetic energy of GRB outflow E 0 , the initial bulk Lorentz factor of GRB outflow Γ 0 , the initial opening angle of the jet θ 0 , the equipartition parameters for the magnetic field and electrons B and e , and the electron distribution index p.
For the first category, considering that the newborn BH is produced by the collapse of a supermassive NS, the initial mass of the BH, M •,0 should be close to the maximum mass of NS, which should be between 2 to 3 solar mass (Lattimer 2012). Within this range, the exact value of M •,0 hardly affect the BZ power. Here we adopt M •,0 = 2.2M as a fiducial value. The influence of a 0 , τ vis and t p are shown in Figure 2. We find that: 1) a larger value of a 0 leads to a more luminous and longer duration plateau; 2) a longer viscosity timescale of disk τ vis results in a longer duration, but lower luminous plateau, which is understandable, since according to Eq (12), a larger value of τ vis corresponds to a slower and weaker BH accretion; 3) a larger value of t p leads to lower luminosity of lightcurves.
For the second category, the initial opening angle of the jet hardly affect the final results. Here we adopt θ 0 = 0.1 as a fiducial value. The influence of Γ 0 , n, B , e and p are shown in Figure 2. We find that: 1) when n is larger, the blast wave collects more material in a higher density ISM, which leads to more luminous lightcurve; 2) since synchrotron radiation intensity increases with the magnetic energy density and the electrons energy in the external shock, the increase of B and e would brighten the X-ray flux; 3) the value of Γ 0 mainly affects the early behavior of the light curve, but not the late behavior when energy injection happens; 4) the value of p mainly affects the decline slope of the light curve.
In order to better reflect the injected energy E inj and the initial kinetic energy of the blast wave E 0 , we test the influence of parameter η ≡ E inj /E 0 instead ofṀ p and E 0 . As expected, the larger of η value, the more significant of the light curve re-brightening (see Figure 2).

Data Reduction And Sample Selection
For the purpose of this work, we systemically search sources consisting two X-ray plateaus, where the first one should be an "internal plateau" followed by a steep decay. The XRT light curves data were downloaded from the Swift/XRT team website 3 (Evans et al. , 2009, and processed with HEASOFT v6.12. There were 1291 GRBs detected by Swift/XRT between 2004 February and 2017 July, with 625 GRBs have well-sampled XRT light curves, which including at least 6 data points, excluding the upper limit. XRT light curves for all selected sample are then fitted with multi-segment broken power law function (in logarithmic scale). Here we adopt the multivariate adaptive regression spline (MARS) technique (e.g. Friedman (1991)) to fit the light curves. MARS technique can automatically determine both variable selection and functional form, resulting in an explanatory predictive model. Some previous works have proven that MARS can automatically fit the XRT light curve with multi-segment broken power-law function (results in general consistent with fitting results provided by the XRT GRB online catalog (Evans et al. , 2009)), detect and optimize all breaks, and record all break times and power-law indices for each segment (see Zhang et al. (2014) and Zhao et al. (2019) for details). Here we treat the adjacent segments with index difference smaller than 0.3 as one component when calculating the segment time span, in order to avoid the potential over-fitting problem from MARS technique. With the fitting results provide by MARS, we searched for candidates having two shallow decay components, where segments with decay slope shallower than 0.65 and time span in log scale larger than 0.4 dex are defined as the shallow decay components 4 . For the purpose of this work, we add one more criteria for the sample selection, i.e., the decay slop following the first shallow decay component should be 3. Eventually, we find 3 long GRBs (GRB 070802,GRB 090111 and GRB 120213A) meeting all our requirement.
GRB 070802 triggered the BAT at 07:07:25 UT on 2007 August 2. T 90 is 16.4 s± 1.0 s. The time-averaged spectrum from T + 4.9 s to T + 23.2 s is best fitted by a single power law (SPL) function and the power law index of spectrum is Γ γ = 1.79 ± 0.27. The fluence in the 15-150 keV band is S γ = 2.8 ± 0.5 × 10 −7 ergs/cm 2 (Cummings et al. (2007)). The XRT began taking data 138 s after the trigger ). The fitting result of XRT data provided by MARS is shown in Figure 3. It is worth noting that on this case, the fitting curve given by the on-line XRT catalogue is a little different from that given by MARS. Since there are few data points in the late stage, the on-line catalogue take all the data after the first plateau into one segment, which is considered as a normal decay component. According to its own algorithm, however, MARS automatically fits the late data with two segments, and the second segment has a decline slope less than 0.65. In this case, we take GRB 070802 as one of our candidates, but with relatively weak evidence. The X-ray fluence and photon index for the two shallow decay segments are listed in Table  3. The UV/Optical Telescope (UVOT; Roming et al. (2005)) start to collect data 141 s after the trigger. No new source in the UVOT observations at the location of the refined XRT position (Kuin & Immler (2007)). Prochaska et al. (2007) observed the afterglow of GRB 070802 with the ESO VLT + FORS. From the detection of several Fe lines in a 30 minutes spectrum starting on August 2.378 UT, the redshift was measured as z = 2.45. Krühler et al. (2008) presented the optical and near-infrared photometry of the afterglow obtained with the multichannel imager GROND. Unfortunately, the late optical data points are also scarce, and there is no data point around the time when the second X-ray plateau emerges. The late optical-IR data could be basically consistent with a single decay segment, so that the optical emission trend might be different with the X-ray band, just like what is found for the case of GRB 070110 (Troja et al. 2007). On the other hand, as shown in the next section, the late optical-IR data could also be well fitted by our proposed model simultaneously with the X-ray data.
The BAT triggered and located GRB 090111 at 23:58:21 UT on 2009 January 11. T 90 is 24.8± 2.7 s. The timeaveraged spectrum from T − 2.9 s to T + 25.6 s is best fitted by a SPL function. The power law index of the time-averaged spectrum is Γ γ = 2.37 ± 0.17. The fluence in the 15-150 keV band is S γ = 6.2 ± 0.6 × 10 −7 ergs/cm 2 (Stamatikos et al. (2009)). XRT observations started at 76.6s after the trigger (Hoversten & Sakamoto 2009). The fitting result of XRT data provided by MARS is shown in Figure 3. The X-ray fluence and photon index for the two shallow decay segments are listed in Table 3. Note that GRB090111 has data more difficult to interpret mainly due to the orbital gap around thousands of seconds, where the data could also be interpreted as a flare followed by a decay (this would also be consistent with the variation in hardness ratio seen in the XRT data).The UVOT to collect data starting 86 s after the trigger. No source was detected by the UVOT at the X-ray afterglow position (Hoversten & Sakamoto (2009)). No prompt ground-based observation was reported, probably due to the vicinity (46 • ) to the Sun (Margutti et al. (2009)).
GRB 120213A triggered the BAT at 00:27:19 UT on 2012 February 12. T 90 is 48.9± 12 s. The time-averaged spectrum from T − 6.31 s to T + 74.46 s is best fitted by a SPL function. The power law index of the time-averaged spectrum is Γ γ = 2.37 ± 0.09. The fluence in the 15-150 keV band is S γ = 1.9 ± 0.1 × 10 −6 ergs/cm 2 (Baumgartner et al. (2012)). The XRT began collect data 54 s after the trigger (Oates & Sakamoto 2012). The fitting result of XRT data provided by MARS is shown in Figure 3. The X-ray fluence and photon index for the two shallow decay segments are listed in Table 3. For this case, there is a fairly clear evidence for a second shallow decay component, but unfortunately the data do not exist to know when the shallow decay stops and at what slop the even later time emission decays. The UVOT to collect data starting 58 s after the trigger. No optical afterglow consistent with the XRT position and no prompt ground-based observation was reported (Oates & Sakamoto (2012)).

Model Application to Selected GRBs
In this section, we apply the model described in Section 2 to interpret the X-ray light curve data of GRBs 070802, 090111, and 120213A. Since the first X-ray plateau data could be well explained with magnetar spin-down power (Troja et al. 2007;Lyons et al. 2010;Rowlinson et al. 2010Rowlinson et al. , 2013Lü & Zhang 2014;Lü et al. 2015), here we focus on fitting the second X-ray plateau data by consider that the second plateau is produced by energy injection from the new born BH driving BZ power. In order to minimize the χ 2 of the fitting, a Markov Chain Monte Carlo (MCMC) method is adopted. In our MCMC fitting, emcee code (Foreman-Mackey et al. 2013) is used with a walkers number 160 and 10 4 burn-in iterations in the ensemble. Considering that the total number of observational data points available are not enough to constrain all the model parameters. In order to reduce the number of the free parameters in our fitting, we fix several parameters at their typical values. For instance, we set E 0 = 10 52 ergs, Γ 0 = 100, n = 1 cm −3 , θ 0 = 0.1, e = 0.1, B = 10 −4 , p = 2.5, and we only take the initial BH spin a 0 , the viscosity timescale of disk τ vis , the peak time of the fallback t p , and the ratio between the injected energy and initial kinetic energy η as free parameters. We set the allowed ranges for the four free parameters in our fitting as: a 0 ≡ [0, 1], log 10 τ vis ≡ [t 0 , t 1 ], log 10 t p ≡ [t 0 , t 1 ], log 10 η ≡ [−3, 3]. Here, t 0 is the ending timescale of internal plateau, which can also be used as the start time of the fall-back accretion and t 1 is the time of last observational data point, the values of them for each burst are shown in Figure 3. Figure 4 shows our fitting results for 3 selected GRBs, where the upper panel shows the fitting light curves and the lower panel shows the corresponding corner plot of the posterior probability distribution for the fitting. We can see that, the second plateau for all 3 GRBs could be well fitted with our proposed model. GRB 070802 has a redshift measurement z = 2.45, and the fitting results at 1σ confidence level are are a 0 = 0.68 +0.15 −0.24 , log 10 τ vis = 4.50 +0.44 −0.36 s , log 10 η = 1.95 +0.54 −0.38 and log 10 t p = 4.12 +0.14 −0.08 s. It is interesting to note that the late optical-IR afterglow data of GRB 070802 could also be well fitted by our proposed model simultaneously with the X-ray data (in Figure 4, we show the data and fitting result for Ks band as an example). For GRB 090111 and GRB 120213A, due to the lack of redshift measurement, here we adopt z = 1 in our analysis. In this case, for GRB 090111, the model parameters at 1σ confidence level are a 0 = 0.7 +0.14 −0.18 , log 10 τ vis = 4.09 +0.08 −0.08 , log 10 η = 1.10 +0.48 −0.39 and log 10 t p = 2.79 +0.19 −0.10 ; for GRB 120213A, the results are a 0 = 0.69 +0.14 −0.19 , log 10 τ vis = 5.43 +0.18 −0.23 , log 10 η = 1.63 +0.44 −0.36 and log 10 t p = 3.95 +0.26 −0.17 . From the fitting results, we find that the constraints on model parameters are relatively loose, mainly due to the lack of enough high quality observation data. Even in this case, some general conclusions could still be made: for all 3 GRBs, 1) the fallback accretion model could easily explain the second X-ray plateau data with fairly loose parameter requirements. Note that in the fitting, we have fixed several model parameters, which means the parameter constraints could become even looser if these parameters were also released. 2) The constraints on η are relatively tight. η was constrained to the order of 10-100, inferring that the BZ power much be 10 or 100 times larger than the initial GRB blast wave kinetic energy, in order to produce the second X-ray plateau feature. 3) Although the allowed parameter space are wide, a 0 tends to have a large value, the distribution peaks are larger than 0.6, which is expected since larger a 0 would easily give larger BZ power. 4) The distributions of the viscosity timescale of disk τ vis are also wide. The distribution peaks of τ vis are relatively large, inferring that the fall-back accretion all falls into the slow accretion regime. 5) The distributions for the peak time of the fallback t p tend to peak around the ending time of internal plateau, which are around 10 3 − 10 4 s. Taken this as the start time of the fall-back accretion, the minimum radius around which matter starts to fall back could be estimated as r fb ∼ 2.8 × 10 11 (M • /2.2M ) 1/3 (t 0 /10 4 s) 2/3 , which is consistent with the typical radius of a Wolf-Rayet star.
The fitting mass fallback rateṀ fb for these three GRBs reaches the peak value around 10 −4 M s −1 at the time about 10 3 -10 4 s. It is interesting to check whether such a mass fallback rate at that time could be supplied by the progenitor envelope. Here we estimate the mass supply rate from the envelope with the presupernova structure models (e.g., Suwa & Ioka 2011;Woosley & Heger 2012;Matsumoto et al. 2015;Liu et al. 2018), i.e., in whichρ = 3M r /(4πr 3 ) is the average density within radius r, M r is the mass coordinate of a shell, t ff = 3π/(32Gρ) = π 2 r 3 /(8GM r ) denotes the free-fall timescale. By taking some representative progenitor density profiles with different metallicities and masses from Liu et al. (2018) and the references therein, we reproduce the mass supply rate changing with time for those progenitor models. In calculation, we set the time when the central accumulated mass reaches M r = M 0 = 2.2M (our fiducial value of the mass of a new born magnetar) as the zero time reference point, i.e., we take t = t ff (r) − t ff (r 0 ), M r = M 0 + r r0 4πr 2 ρdr, here r 0 is the radial coordinate where the enclosed mass is M 0 . As shown in Figure 5, we find that our fitting resulted mass fallback rates are compatible with the theoretical mass supply rate of some low metallicity massive progenitor stars such as those ones with (Z 10 −1 , M 80M ), (Z 10 −2 , M 40M ), (Z 10 −4 , M 20M ), and so on. While the solor metallicity stars might not be so possible to play as the progenitors for our sample.
Based on the fitting results, the magnetic field strength of the new-born BH (B • ) for these three GRBs reaches the peak value around 10 13−14 G at the time about 10 3 -10 4 s. According to the dipole spin-down model, one can make estimation for the surface magnetic field B p and the initial spin period P 0 of the rapidly spinning magnetar with the first plateau data for all 3 GRBs in our sample (Rowlinson et al. 2013;Lü et al. 2015). Here we adopt the constant values of the moment of inertia I = 1.5 × 10 45 g cm 2 and radius R=10 km for a typical neutron star. For GRB 070802, the plateau luminosity and break time are L b = 1.4 × 10 47 erg s −1 , t b = 1 × 10 4 s, respectively, one can  thus derive 5 P 0 < 4.64 × 10 −3 s and B p < 2.59 × 10 15 G. For GRB 090111, the plateau luminosity and break time are L b = 2.9 × 10 47 erg s −1 , t b = 427 s, respectively, one can derive P 0 < 1.55 × 10 −2 s and B p < 4.19 × 10 16 G. For GRB 120213A, the plateau luminosity and break time are L b = 2.2 × 10 47 erg s −1 , t b = 4.57 × 10 3 s, respectively, one can derive P 0 < 5.44 × 10 −3 s and B p < 4.49 × 10 15 G. We find that the magnetic field strength of the new-born BH required to power BZ jets is comparable or slightly lower than the magnetic field strength of the magnetar, which is understandable if we consider that the magnetic flux should be roughly conserved when the magnetar collapse into the BH, and some magnetic energy might dissipate due to the interaction between the magnetosphere and the fall back flow (Lloyd-Ronning et al. 2019).

DISCUSSION AND CONCLUSION
In the 15 years of Swift's operation, it has brought us a lot of observations of GRB X-ray afterglow, which provides valuable information for understanding the GRB central engine. One particular example is the discovery of internal X-ray plateau, (a plateau followed by a very steep decay phase), which is commonly taken as the smoking gun evidence of a rapidly spinning magnetar as the central engine. The very steep decay at the end of the plateau suggests a sudden cessation of the central engine, which is explained as the collapse of a supra-massive magnetar into a black hole when it spins down. If this interpretation is correct, the fall-back accretion from the envelope of progenitor star into the newborn BH could generate some detectable signatures.
Here we propose that the energy extracted from the newborn BH could be continuously injected into the GRB blast wave. In this scenario, We find that with appropriate parameters for the fall back accretion and new born BH, it is possible to produce a second plateau following the steep decay phase of the internal plateau. With a systematical search through the Swift -XRT sample, we find three interesting long GRBs, i.e., GRBs 070802, 090111, and 120213A, whose X-ray afterglow light curves contain internal X-ray plateau followed by a second plateau. Here we focus on fitting the second X-ray plateau for these 3 GRBs with our proposed model. We find that in a fairly loose and reasonable parameter space, the second X-ray plateau data could be well interpreted with our model. It is worth noting that the quality of current observation data is not good enough to completely eliminate the degeneracy of parameters. In our sample, GRB 090111 has the most observational data points in the period of the second plateau. GRBs 070802 and 120213A only have four data points in the period of the second plateau. Even in this case, some general conclusions could still be reached, for instance, in order to interpret the second plateaus, the initial spin of the new-born BH tends to have a large value (the peak of its posterior probability distribution is larger than 0.6), the later injected energy should be 10 or 100 times larger than the initial GRB blast wave kinetic energy, and the viscosity timescale of disk tends to be large, inferring that the fall-back accretion all falls into the slow accretion regime . The mass fall-back rate reaches the peak value around 10 −4 M s −1 at the time about 10 3 -10 4 s, which is compatible with the late mass supply rate of some low metallicity massive progenitor stars. Based on the fitting results, one can infer the total accreted masses M acc ∼ 0.08 − 0.26M and the fallback radii r fb ∼ a few × 10 11 cm, which is consistent with the typical radius of a Wolf-Rayet star. Combining the X-ray data of the internal plateaus of the three GRBs with the magnetic dipole radiation model of the magnetar, we coarsely estimate the strength of magnetic field of the magnetar before its collapsing to a BH. We find the magnetic field is strong enough to drive a BZ jet.
If our interpretation is correct, the 3 GRBs with two X-ray plateaus provide additional evidence for rapidly spinning magnetar as the GRB central engine. It is worth noticing that most GRBs with internal plateaus do not show a second X-ray plateau, which means for most cases the fallback accretion may be relatively weak, so that the injected energy is smaller than initial kinetic energy of GRB blast wave. Future observations are likely to discover more similar events, which could offer more information of the properties of the magnetar as well as the newborn BH.
In this paper, we ignore the mass into the outflow from disk, which will reduce the accretion rate onto BH horizon during late central engine activity. Usually, the distribution of accretion rate with disk radius is simply described with a power-law model due to the poor knowledge of the disk outflow. The effects of such outflow are thus highly relied on the uncertain power-law index parameter. The existence of disk outflow may also be important to comprehend the baryon loading into GRB jet  and 56 Ni synthesis for associated supernovae (Song & Liu 2019). We hope future general-relativistic magnetohydrodynamic (GRMHD) simulation for a better understanding.